Mailman 3 March 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

sfepy installation error
by Kid Guo 31 Mar '15

31 Mar '15

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

2 1

Looking to contribute to Sfepy
by Alex Eftimiades 29 Mar '15

29 Mar '15

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

2 7

Recommended solvers for Laplace problem?
by Jonatan Lehtonen 25 Mar '15

25 Mar '15

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.

2 10

Windows 7 installation problem
by berat.s...＠gmail.com 13 Mar '15

13 Mar '15

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

3 5

import error
by Sumesh K.C. 03 Mar '15

03 Mar '15

>>> 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.

2 1