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.
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
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
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-deve...https://groups.google.com/forum/#!searchin/sfepy-devel/elastic$20periodic...
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
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.
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.
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
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
