Mailman 3 August 2008 - SfePy

large deformations
by Robert Cimrman Sept. 1, 2008

Sept. 1, 2008

Hi all, and Ryan in particular, I have started implementing large deformations into the code. Nothing is finished yet, but I have already some proof-of-concept. I will commit relatively soon (i.e. after my holidays) a patch adding the hyperelastic neo-Hookean material. Ultimately, I would like to have in SfePy all the functionality I had in its predecessor (a matlab code called mafest1), but it is a long term goal, as it was more a one-purpose tool than a general FE system. All that the old code can do is nicely summarized in my thesis, that you can find at http://ui505p06-mbs.ntc.zcu.cz/sfe/RobertCimrman at the bottom (thesis.pdf). Ryan, try to look in particular at Section 2.1.2 (TL formulation), Chapter 3 (FE discretization in TL context) and Chapter 6 (FE discretization in a linear problem context), as it may provide you some insight. Appendix A might be also interesting, as the usual constitutive (stress-strain) relations are summarized there, e.g. the neo-Hookean one. r.

2 3

make sfepy usable a'la mathematica using sympy
by Ondrej Certik Sept. 1, 2008

Sept. 1, 2008

Hi, I'd like to be able to do these things with sfepy+sympy together: http://reference.wolfram.com/mathematica/tutorial/NDSolvePDE.html I.e. really just specify the equation in sympy and sfepy should be clever enough to solve it correctly, or say why it cannot solve it. Ondrej

2 5

next release
by Robert Cimrman Aug. 27, 2008

Aug. 27, 2008

Hi all, it's time to make another release (this week or start of the next one), so if you have time and will, try the latest hg version, please, and report any problems. Mainly the new hyperelastic terms surely need some testing (input/hyperelastic.py), as well as solutions in complex domain (input/complex.py). r.

1 0

Compilation/configuration problems
by Michael Aug. 24, 2008

Aug. 24, 2008

Hello everybody, while compiling sfepy I found that makefile was not able to find a correct include directory for the numphy library. I had to change it by hand and it worked. I set: PYTHON_INCL := -I/usr/include/python$(PYVER) -I/usr/lib64/python2.4/ site-packages/numpy/core/include/ My Linux : uname -s -r -p -i -m Linux 2.6.18-92.1.6.el5 x86_64 x86_64 x86_64 The numphy library I compiled from the distribution numpy-1.1.1.tar.gz. I just followed their instruction: sudo python setup.py install. Michael.

2 1

holidays etc.
by Robert Cimrman Aug. 13, 2008

Aug. 13, 2008

Hi all, I will be hiking in Alps from 14. to 24. August, so keep posting your questions/bugreports/patches, I will process the unresolved ones as fast as possible after. This week is unfortunately an IT infrastructure maintenance week at our university, so I may have problems getting on-line even before Thursday. It is already very slow. best wishes, r.

4 3

Installing on a Mac + manual
by Dominique Aug. 12, 2008

Aug. 12, 2008

Hi all, I just downloaded the latest devel version of sfepy for installation on an Intel Mac. I wanted to let you know that in order for sfepy to compile properly, I had to make the following changes: - In the top directory Makefile: add the following to OPTFLAGS -fno-strict-aliasing -fno-common -dynamic - In sfepy/fem/extmods/Makefile.inc and sfepy/terms/extmods/ Makefile.inc, replace occurences of '-shared' with $(SHARED), where I defined SHARED = -dynamiclib -undefined dynamic_lookup at the top of the Makefile.inc (there surely is a nicer way!) I thought this might be useful to others. Next, upon running the tests, I obtained: [dpo@pod:sfepy-fe15c3e45e20]$ ./runTests.py --filter-more <<< directory: tests, test files: 23 +++ test_vector_matrix: ok +++ test_struct_add: ok +++ test_struct_i_add: ok +++ test_get_solution: ok +++ test_linear_terms: ok +++ test_input: ok +++ test_input: ok +++ test_input: ok nodes: 100% | #####################################################################| Time: 00:00:00 elements: 100% | ##################################################################| Time: 00:00:00 +++ test_input: ok +++ test_input: ok +++ test_input: ok +++ test_input: ok +++ test_input: ok +++ test_recursive_dict_hdf5: ok +++ test_sparse_matrix_hdf5: ok --- test_linear_rigid_body_bc: failed! --- test_linear_rigid_body_bc: failed! +++ test_solvers: ok nodes: 100% | #####################################################################| Time: 00:00:00 elements: 100% | ##################################################################| Time: 00:00:00 +++ test_read_meshes: ok +++ test_compare_same_meshes: ok +++ test_read_dimension: ok --- test_msm_laplace: failed! +++ test_msm_symbolic_diffusion: ok +++ test_msm_symbolic_laplace: ok +++ test_parse_equations: ok +++ test_parse_regions: ok +++ test_pbc: ok +++ test_pbc: ok +++ test_problem_creation: ok +++ test_consistency_d_dw: ok 23 test file(s) executed in 27.13 s, 3 failure(s) of 30 test(s) How can I check what went wrong with the three failures? I see no log files in the 'output_tests' directory. Moreover, I can't seem to do anything with the vtk files generated by the tests, neither with that generated with 'simply.py input/ poisson.py'. What am I missing? I have the latest Enthought Python Distribution and Enthought Tool Suite. Mayavi2 cannot seem to display the vtk scenes and neither can Paraview. Any idea what's wrong? Note that in order for the tests to run, I had to bypass the Scipy distributed as part of ETS since it is outdated (!). I had a build of 0.7.0 lying around, with UMFPACK, and that seemed to work with sfepy. Finally, yes I know, I really should read the manual. But I can't generate it! While pexpect is a breeze to install, lxml is a real pain on the Mac. Could anyone post the sfepy manual on the Wiki or on the newsgroup page? Many thanks, Dominique

3 6

Paraview question: wireframe labeling points
by Ryan Krauss Aug. 11, 2008

Aug. 11, 2008

I have poked around in google for a while and even tried mayavi2 and haven't come up with a good answer to this slightly off topic question. I am trying to create a tutorial for myself (and hopefully others) on simple usage of SfePy for an axially loaded bar. I created a simple mesh of two bricks connected in series (using blockgen.py) and now I would like to show my mesh and label the points (to fully illustrate the vtk file). I can't find a simple way to generate a wireframe and add point #'s so that these lines of the vrk file CELLS 2 18 8 0 6 9 3 1 7 10 4 8 1 7 10 4 2 8 11 5 make sense (i.e. I would like to generate a wireframe with point #'s 0-11 labeled). Is there an easy way to do this using paraview (or someother program)? Thanks, Ryan

2 1

trying to understand linear elasticity
by Ryan Krauss Aug. 5, 2008

Aug. 5, 2008

So, I dug into sfepy/sfepy/terms/termsLinElasticity.py(61)__call__() a bit to try and understand how linear elasticity is being handled: def __call__( self, diff_var = None, chunk_size = None, **kwargs ): Pdb().set_trace() material, virtual, state = self.get_args( **kwargs ) ap, vg = virtual.get_approximation( self.get_current_group(), 'Volume' ) shape, mode = self.get_shape( diff_var, chunk_size, ap ) fargs = self.build_c_fun_args( material, state, ap, vg ) for out, chunk in self.char_fun( chunk_size, shape ): status = self.function( out, *fargs + (chunk, mode) ) yield out, chunk, status I think I understand the first line: material, virtual, state = self.get_args( **kwargs ). It seems to extract the material, including its Lame parameters and all that, along with the state. I am not entirely sure what state means: ipdb> print state Variable:u kind: unknown _order: 0 name: u family: field dual_var_name: v dpn: 3 dof_conns: Struct:dof_conns indx: slice(0, 1062, None) dof_name: u dtype: <type 'numpy.float64'> field: Field:3_displacement step: None eq_map: Struct flags: set([0, 10]) i_dof_max: 1062 key: variable_1 current_ap: None n_dof: 1062 dofs: ['u.0', 'u.1', 'u.2'] data: deque([array([ 0., 0., 0., ..., 0., 0., 0.])]) history: None Some of those terms make sense. But what are deque and flags? What is the indx used for? Is that to get this node out of the mesh? But it seems like the real action happens when self.function is called. self.function refers to sfepy/sfepy/terms/extmods/terms.py(166)dw_lin_elastic_iso() which just calls a C function. I am scared that implementing a nonlinear material model will require C programming. I haven't done much C lately. Ryan

2 6

Help expressing 1d eigenproblems in sfepy
by rpmu...＠gmail.com Aug. 5, 2008

Aug. 5, 2008

Here's yet another plea for help on a "first problem". I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here: http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like: - create a mesh - work out a bunch of integrals over the mesh elements - then a miracle occurs - write down your final solution Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well. My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer. Rick

5 17

understanding my first problem
by Ryan Krauss Aug. 1, 2008

Aug. 1, 2008

So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing. My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done. Referring to the attached file, which gets passed to simple.py as an input, I have the following questions: The file starts with field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} } variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', } All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)? >From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else? I assume 'Omega' is a domain that includes the whole problem, but that I can also name domains and some how use different domains for different purposes. Is this correct? What does the '3_4_P1' mean? Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is ultimately being sought, reported, and visualized? What are the units on 'u'? ====================== The next section I have questions about is this: equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", } material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. } What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related. It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct? It seems like this function: def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt ) print val return {'val' : val} defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it? I think that is enough for now. Thanks for your ongoing help and patience. Ryan

4 20