Mailman 3 April 2015 - SfePy

Re: Torque
by Robert Cimrman 09 Mar '18

09 Mar '18

In examples/large_deformation/hyperelastic.py a rotation by displacements is applied. By using a similar function the vectors defining the force couples could be defined for dw_surface_ltr (IMHO). Does it make sense? r. ----- Reply message ----- From: "Andre Smit" <freev...(a)gmail.com> To: <sfepy...(a)googlegroups.com> Subject: Torque Date: Sat, Dec 18, 2010 05:10 What is the best way to apply a torque load to a model? -- Andre -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...(a)googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.

4 7

new interactive example with time stepping
by Robert Cimrman 13 Mar '17

13 Mar '17

Hi, I have just updated the time stepping solvers in sfepy for interactive use, as demonstrated in the new example [1]. For basic use, ignore the probing code - the time stepper can be used as simply as: tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 100.0, 'n_step' : 11}, problem=problem) tss.init_time() for step, time, state in tss(): pass r. [1] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html

2 8

1D line elements in SfePy
by Nimish 22 Jan '16

22 Jan '16

I am currrently looking for FEM packages to help me solve a system of beams and columns, basically a collection of 1D bernoulli/timoshenko line elements. I started reading SfePy docs and i am getting the idea that doing the above is not really possible here, am i right? Are only 2D area elements permitted in SfePy? Or is there any direct support for solving 1D line elements too.. Cheers Nimish

4 9

Elasticity with periodic boundary condition
by Ronghai Wu 03 Jun '15

03 Jun '15

Dear Robert Cimrman, I am trying to implement elastic inclusion problem. It is a 2D mesh with an inclusion(and therefore eigenstrain) at mesh centre. I developed my codes(as attached) based on following two previous issues. https://groups.google.com/forum/#!searchin/sfepy-devel/Ronghai/sfepy-devel/… https://groups.google.com/forum/#!searchin/sfepy-devel/elastic$20periodic$2… But I have troubles with the boundary conditions. I would like to try two kinds of boundary conditions: (1) Top and bottom boundaries: uniaxial loading(or surface traction). Left and Right boundaries: periodic. I used "dw_tl_surface_traction" but some errors come out. (2) Periodic at all boundaries. I think there must be a easy way for this, but failed to find out in SfePy documentation and mailing list. Could you please give me some help. Thanks. Regards Ronghai

2 19

Probing occasionally returns nan for no apparent reason
by Jonatan Lehtonen 30 Apr '15

30 Apr '15

Hi, I'm using SfePy to solve a simple Laplace problem in a unit square with Dirichlet boundary conditions (well, what I'm actually doing is a bit more involved, but this is just a minimal example). I'm having some problems with the fact that probing with a PointsProbe occasionally returns nan instead of a proper value. I have attached a mesh file which I created by using gmsh with the attached .geo file and then manually converting it to the desired format for SfePy (that is, I removed the edges and the z-coordinates of the nodes). The following code, with the attached mesh file, should be able to reproduce this issue: import numpy as np from sfepy.discrete import (FieldVariable, Integral, Equation, Equations, Problem) from sfepy.discrete.fem import (Mesh, FEDomain, Field) from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.postprocess.viewer import Viewer mesh = Mesh.from_file('rectangle2D.mesh') domain = FEDomain('domain', mesh) omega = domain.create_region('Omega', 'all') # Groups 1 and 3 correspond to the bottom and top of the square, respectively gamma0 = domain.create_region('Gamma 0', 'vertices of group 1', 'facet') gamma1 = domain.create_region('Gamma 1', 'vertices of group 3', 'facet') field = Field.from_args('fu', np.float64, 'scalar', omega, approx_order=2) u = FieldVariable('u', 'unknown', field) v = FieldVariable('v', 'test', field, primary_var_name='u') integral = Integral('i', order=2) term = Term.new('dw_laplace(v, u)', integral, omega, v=v, u=u) eq = Equation('Laplace equation', term) eqs = Equations([eq]) ebc0 = EssentialBC('Dirichlet 0', gamma0, {'u.all': 0.0}) ebc1 = EssentialBC('Dirichlet 1', gamma1, {'u.all': 1.0}) pb = Problem('Laplace problem', equations=eqs) pb.time_update(ebcs=Conditions([ebc0, ebc1])) vec = pb.solve() pb.save_state('rectangle2D.vtk', vec) view = Viewer('rectangle2D.vtk') view() from sfepy.discrete.probes import PointsProbe probe = PointsProbe(np.random.rand(100000,2)) pars, res = probe.probe(pb.get_variables()['u']) print(np.any(np.isnan(res))) print(probe.points[np.nonzero(np.isnan(res))[0]]) The code above solves the Laplace problem and creates a PointsProbe with 100000 random points. It then prints True if any of the probed values equal nan, and prints the coordinates of the corresponding points. When I run this code, it usually turns up around 5-10 points which produce a nan value. This is an example of the output on my system: sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (rectangle2D.mesh)... sfepy: ...done in 0.03 s sfepy: updating variables... sfepy: ...done sfepy: setting up dof connectivities... sfepy: ...done in 0.00 s sfepy: matrix shape: (13343, 13343) sfepy: assembling matrix graph... sfepy: ...done in 0.01 s sfepy: matrix structural nonzeros: 151301 (8.50e-04% fill) sfepy: updating materials... sfepy: ...done in 0.00 s sfepy: nls: iter: 0, residual: 1.819854e+01 (rel: 1.000000e+00) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.07 [s] sfepy: matrix: 0.01 [s] sfepy: nls: iter: 1, residual: 4.729175e-14 (rel: 2.598656e-15) sfepy: point scalars u at [-0.5 -0.5 0. ] sfepy: range: 0.00e+00 1.00e+00 l2 norm range: 0.00e+00 1.00e+00 sfepy: iconn: 0.016510 s sfepy: kdtree: 0.000420 s sfepy: evaluating in 100000 points... sfepy: reference field: 3436 vertices sfepy: kdtree query: 0.291867 s sfepy: ref. coordinates: 0.466726 s sfepy: interpolation: 0.135751 s sfepy: ...done True [[ 0.48132644 0.32787828] [ 0.1449193 0.69671521] [ 0.48135807 0.32673518] [ 0.77758596 0.35547628] [ 0.64769394 0.79028833] [ 0.66355676 0.78836725] [ 0.4799345 0.32667195] [ 0.48147468 0.3281688 ]] After some further examination, I've found that the problem occurs in some small triangular regions of the domain. To illustrate this problem, I made a (rather ugly, but possibly helpful) scatter plot of the values in a grid around the first point of the list above, using this (extremely messy) code: temp = probe.points[np.nonzero(np.isnan(res))[0][0]] X, Y = np.mgrid[temp[0]-0.003:temp[0]+0.003:40j, temp[1]-0.003:temp[1]+0.003:40j] gridprobe = PointsProbe(np.array(np.array([X.flatten(), Y.flatten()]).transpose(), order='C')) plot_data = gridprobe.probe(pb.get_variables()['u'])[1].flatten() from pylab import * figure() scatter(X.flatten()[~np.isnan(plot_data)], Y.flatten()[~np.isnan(plot_data)]) closest_points = vec.variables['u'].field.get_coor(np.argsort(np.linalg.norm(vec.variables['u'].field.get_coor() - temp, 2, axis=1))[:4]) scatter(closest_points[:,0], closest_points[:,1], color='r') scatter([temp[0]], [temp[1]], color='g') show() The resulting plot is attached. The triangular area inside the blue grid is where probing results in nans. The green point ('temp' above) is the first nan point which was found earlier (i.e. the point (0.48132644 0.32787828)), and the red points are the four mesh nodes closest to the green point. I included the mesh points because my first assumption was that the nan area would correspond to a mesh triangle, but that does not seem to be the case. Now the obvious question at this point is: am I doing something wrong, or is this a bug? Regards, Jonatan P.S. I was originally using evaluate_at() instead of a PointProbe, and that function returns 0.0 instead of nan. I'm not entirely sure if this is intended, but at least I found it extremely confusing.

2 12

major update in master - no element groups
by Robert Cimrman 17 Apr '15

17 Apr '15

Hi, I have just merged the branch 'no-ig' (aka element groups, go away). The element groups served no particular purpose for some time already, because the fields required the same element type in the whole mesh anyway. As you can see in the merge message below, I was mostly deleting and simplifying code: 63 files changed, 1761 insertions(+), 3931 deletions(-) The functionality stays, while there should be a lower memory footprint and the code is slightly faster (about 10% when running tests on my computers). Your code will probably need an update, if you access material data or physical quadrature points directly. If you never needed to know about element groups, you probably do not need to do anything. There might be also some left-overs in the code/docs, so let me know about any problems or inconsistencies. r.

1 0

Region-specific node values of a variable
by Mark Davis 15 Apr '15

15 Apr '15

I'm just starting to use this package, and I'm enjoying it so far. I'm working on a problem where two adjacent regions share an interface, and I would like to list variable values at the actual nodes along the interface in each region. Is there an easy way to do this? Regards, Mark

2 3

Node variable values along a surface boundary
by Mark Davis 14 Apr '15

14 Apr '15

Hello. I'm new to this package and am trying to find my way around. It seems to be great work with lots of capability. Is there an easy way to extract node values from a FieldVariable in a specific region? I have a shared internal surface between two regions, and I want to extract variables at the surface between the two regions. "evaluate_at" with interface coors works for one region, but not the other. If the boundary on one side is a Neumann boundary, is it possible that variable data is not saved on those boundary nodes? Any help is appreciated. Thanks, Mark

2 1

Why it removed biot in the new version sfepy?
by Kid Guo 07 Apr '15

07 Apr '15

I want to create a model by biot,and i can't look at the link of biot/biot.py.What should I do? Thank you advanced Kid Guo

2 1