Re: Torque
by Robert Cimrman
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.
1 year, 1 month
1D line elements in SfePy
by Nimish
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
3 years, 3 months
Probing occasionally returns nan for no apparent reason
by Jonatan Lehtonen
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.
3 years, 11 months
sfepy installation error
by Kid Guo
Hi,I am trying to install sfepy from
http://sfepy.org/doc-devel/installation.html#table-of-contents,and when I
got to this step(step 6) in Steps to Get a Working *SfePy* on Windows Using
Python(x,y) <http://sfepy.org/doc-devel/installation.html#table-of-contents>,I
got an error,like this
Administrator@DELL /c/sfepy
$ python setup.py build_ext --inplace --compiler=mingw32
Traceback (most recent call last):
File "setup.py", line 111, in <module>
package_check('tables', INFO.PYTABLES_MIN_VERSION)
File "c:\sfepy\build_helpers.py", line 311, in package_check
raise RuntimeError(msgs['missing'] % pkg_name)
RuntimeError: Cannot import package "tables" - is it installed?
谷歌(国内)翻译结果谷歌(国外)翻译结果有道翻译结果百度翻译结果网页翻译
Then I reinstall but it also existed.What happend?What should I do?
Could you please give me some help. Thanks.
Regards
Kid guo
4 years
Looking to contribute to Sfepy
by Alex Eftimiades
I was working on a discrete exterior calculus framework about a year ago.
When it came time to start benchmarking it, I was disappointed with my
options for mesh generation. I wanted to test my framework on different
geometries, topologies, and dimensions in a consistent way. Since I could
not find a good way of doing this, I wrote my own
<https://github.com/aeftimia/pyfeec/blob/master/pyfeec/grid_utils.pyx> in
Cython. I apologize for the poor documentation. It turns out I ended up
reinventing "Kuhn triangulation" as documented in section 2.2.1 and 2.2.2
in this thesis
<https://dspace.library.cornell.edu/bitstream/1813/6140/1/92-1322.pdf>. I
am hoping to get more involved collaborating with open source projects as
opposed to working on my own projects. I noticed Sfepy distinctly noted
lacking BSD licensed mesh generation code in a 2008 power point
<http://www.ondrejcertik.com/media/euroscipy2008.pdf> but seemed to have
made some progress since then as documented here
<http://sfepy.org/doc-devel/src/sfepy/mesh/mesh_generators.html>.
My implementation allows you to specify an arbitrary number of subdivisions
along an arbitrary number of axes, and contains two algorithms for
generating a corresponding simplicial mesh. Pictures of the result of each
algorithm for a 2D discretization are attached. The process looks like this:
1. Create a dict that maps tuples of integers (number of "steps" along each
axis) to a global index of the corresponding vertex.
>> shape = [nx, ny, ....]
>> g = grid_indices(shape)
2. Connect the vertexes according to one of two algorithms for generating a
uniform simplicial mesh (the difference in the resulting meshes is readily
seen in the attached pictures). This is stored as an NxD dimensional numpy
array of unsigned integers for a mesh with N D-dimensional simplexes. Each
integer in the array is the global index of a vertex in the mesh.
>> simplices = grid(g)
3. Customize topology of mesh by stitching pairs of vertexes together. This
is accomplished by replacing each occurrence of an index in the above array
with the index of the vertex you want to "stitch" it to. So in the case of
a one dimensional chain of 10 vertexes, I could create a "ring"-like
topology by replacing each occurrence of the integer "9" (the rightmost
vertex) with a "0" (the leftmost vertex) in my 9x1 dimensional numpy array
of integers (as described in step 2). Since periodic boundary conditions
are very common, I wrote a routine specifically for toroidal topologies
that are periodic along arbitrary combinations of axes. The routine returns
a dict that maps vertex indexes to other vertex indexes such that the key
index is to be replaced with the value index. So in this example, it would
return {9: 0} to indicate vertex 9 is to be replaced with vertex 0.
>> stitches = pbc_stitches(g, shape, [0,3]) #imposing periodic boundary
conditions along the first and fourth axes.
4. I found it most useful to keep the original "unstitched" simplexes too
so geometry can be specified separately by embedding the vertexes within a
specific coordinate system. For example, I might wish to impose periodic
boundary conditions along the second axis in a 3 dimensional system and
leave the system in cartesian coordinates. I might alternatively choose to
embed the points in a cylindrical coordinate system and have the angular
coordinate periodic.
I don't know if you would find this useful. If not, is there anything else
I can do to contribute to Sfepy? I read the how to contribute guide
<http://sfepy.org/doc-devel/developer_guide.html#how-to-contribute>, but
your development section did say to let you know on your mailing list that
I was interested in contributing to Sfepy. I am generally looking to gain
experience contributing to open source projects. Since I spend a lot of
time writing Python programs that numerically solve PDEs (typically
eigenvalue problems at the moment), I thought sfepy might be a good place
to start.
Thanks,
Alex Eftimiades
4 years
Recommended solvers for Laplace problem?
by Jonatan Lehtonen
Hello,
I am trying to solve a 2D Laplace problem with simple Dirichlet boundary
conditions (u = 0 on one part of the boundary, u = 1 on another). I
mentioned this same problem in my previous post
<https://groups.google.com/forum/#!topic/sfepy-devel/8xgzvEc6Noo>; the only
difference is that I want to be able to solve it in a much more general
domain than a unit square. The code I'm currently using looks like this:
mesh = mesh_file_reader(filepath, bc_indices)
domain = FEDomain('domain', mesh)
omega = domain.create_region('Omega', 'all')
gamma0 = domain.create_region('Gamma 0', ' +v '.join('vertices of group
{}'.format(n) for n in bc_indices[0]),
'facet')
gamma1 = domain.create_region('Gamma 1', ' +v '.join('vertices of group
{}'.format(n) for n in bc_indices[1]),
'facet')
field = Field.from_args('fu', nm.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()
While it's not relevant to the problem at hand, here's two notes about the
code above which may seem odd at first glance. Firstly, mesh_file_reader is
a custom function which simply reads the necessary data from a .mesh-file
generated by gmsh and calls Mesh.from_data() with the appropriate
arguments. Secondly, bc_indices contains two tuples of indices which define
parts of the boundary; this is just to work around some shortcomings with
the gmsh command line interface, but essentially domain.create_region() is
called with something like "vertices of group 1 +v vertices of group 2".
Anyway, my goal is to get the highest possible accuracy out of this. The
only real limitations here are memory and speed. So basically my question
is, given a fixed amount of memory and time, what would be the best choice
of solver methods for this problem? Should I just go with the defaults?
Also, is there some recommendation for settings such as the mesh order?
Memory is my biggest concern; for time consumption something which takes an
order of magnitude longer than the default direct solver would still be
acceptable.
To clarify, I fully understand that there's no silver bullet that magically
solves all my problems, I'm just looking for some pointers about what I
should try, as my own experience in FEM solvers is quite limited. Thus far,
I've successfully solved the problem with the defaults, and I've tried all
the available SciPy iterative solvers to no avail. WIth the iterative
solvers, convergence seems much too slow; this is likely because I did not
supply a preconditioner, but I'm not entirely sure how that should be done
in SfePy, or even what would be an appropriate preconditioner in the first
place.
I would greatly appreciate any help with this problem. I can probably work
with what I already have, but any performance improvements would certainly
be most welcome.
Regards,
Jonatan
PS. An additional somewhat unrelated question. My goal is to solve several
hundred of these problems in one execution; each time, the geometry of the
mesh changes slightly, but everything else stays the same. Once I start
solving a new problem, I don't need the old one anymore, so is there a way
of deleting all the SfePy data related to the old problem? That is, I would
like to be able to run the code I posted above several hundred times in a
row (with different file paths for the mesh_file_reader), without running
into memory leak issues. Using the same mesh each time, running the code I
posted above results in a peak memory consumption of about 600MB (as
reported by the Windows Task Manager) the first time I run it. When I rerun
it, the peak memory usage jumps up to 900MB, and after each subsequent
repetition the peak memory usage increases by 30MB until eventually I reach
1.3 GB and UMFPACK fails due to running out of memory. As each rerun causes
all the variables above to be reassigned, the old ones should be garbage
collected, so the only reason I can think of for this happening is that
SfePy stores the data somewhere internally. It would be great if there was
a way of essentially telling SfePy to start from a clean slate, without
having to do trickery like calling SfePy in a temporary subprocess or
something like that.
PPS. On the subject of things getting left behind from previous problem
solutions, I've noticed that the Probe cache isn't reset when a new problem
is solved. This is not unexpected, I'd just like to check whether my
workaround for this is correct. What I've done is include the following
code after the one I posted above:
from sfepy.discrete.probes import Probe
Probe.cache = pb.domain.get_evaluate_cache(cache=None, share_geometry=True)
Is this the correct way of going about this? The goal is to ensure that
subsequent calls to PointsProbe() use the correct cache without needing to
explicitly give them a cache object.
4 years
Windows 7 installation problem
by berat.s...＠gmail.com
Hi!
I sent two posts regarding to this issue on Sunday, but I cannot see any
of my posts in the list. Therefore, I decided to send a new post in the
hope that it will be visible and has been received.
If somebody is organizing the posts under the hood, please send me an email
that my post has been received.
I followed the installation instructions for Windows 7 and I started
installing Python(x,y), version: 2.7.9.0 . As recommended, I made a "FULL"
installation and the installation went smoothly. Python(x,y) is now up and
running.
Then I started the next step: the installation of "pyparsing" when I ran
the installer, I got the following error message: "Python version 2.6
required, which was not found in the registry".
The "pyparsing" I downloaded has the following version signature:
pyparsing-2.0.3.win32-py2.6. I have spent sometime to investigate the
issue. Unfortunately, there is no information as to what the problem might
be.
What I assume is that there is probably a compatibility issue between 64
bit Python(x,y) and 32 bit pyparsing, but I am not sure.
Now the question is : what do I do next?
I should also add that these pages are unstructured, meaning the posts here
are sorted with respect to the date they were posted and not grouped with
respect to their content. If these posts are categorized, it makes it easy
for people to find the information they need.
I spent sometime reading the posts here to find information about, if
somebody encountered previously similar problem that I have now, before I
wrote this post, but going through all the posts here was and still is a
daunting experience.
I thank you very much for any help.
Regards,
Berat
4 years, 1 month
import error
by Sumesh K.C.
>>> from mayavi import mlab
Traceback (most recent call last):
File "<pyshell#4>", line 1, in <module>
from mayavi import mlab
File "C:\Python27\ArcGIS10.1\lib\site-packages\mayavi\mlab.py", line 27,
in <module>
from mayavi.tools.camera import view, roll, yaw, pitch, move
File "C:\Python27\ArcGIS10.1\lib\site-packages\mayavi\tools\camera.py",
line 25, in <module>
from engine_manager import get_engine
File
"C:\Python27\ArcGIS10.1\lib\site-packages\mayavi\tools\engine_manager.py",
line 9, in <module>
from traits.api import HasTraits, Instance
File "C:\Python27\ArcGIS10.1\lib\site-packages\traits\api.py", line 26,
in <module>
from .trait_base import Uninitialized, Undefined, Missing, Self,
python_version
File "C:\Python27\ArcGIS10.1\lib\site-packages\traits\trait_base.py",
line 283, in <module>
from . import ctraits
ImportError: DLL load failed: The specified procedure could not be found.
4 years, 1 month
