Mailman 3 July 2008 - SfePy

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

Comments added to le.py
by Ryan Krauss Aug. 1, 2008

Aug. 1, 2008

I added some comments to le.py that I think are helpful to other beginners. If they are wrong, I need to know what I am still misunderstanding about the setup file. If they are right and useful, let me know what the best format is for doing that sort of thing. My hg repo is here: http://freehg.sympy.org/u/ryankrauss/sfepy (I had to recreate it for some reason). Ryan

3 2

how to install paraview in ubuntu
by Ondrej Certik July 31, 2008

July 31, 2008

Hi, since couple of you had problems installing paraview in ubuntu, here is how to do it: Add these lines to /etc/apt/sources.list: deb http://cz.archive.ubuntu.com/ubuntu/ intrepid main restricted deb http://cz.archive.ubuntu.com/ubuntu/ intrepid universe (you can change the "cz" to "us" or any other mirror) Then: $ sudo apt-get update $ sudo apt-get install paraview And that's it. No hassle, no setups, no unofficial packages. I tested it on hardy in virtualbox (in Debian) and Ryan tested this on gutsy which he uses. Ondrej

2 1

Problem after HG update
by Ryan Krauss July 31, 2008

July 31, 2008

I am having problems running my script after an HG update this morning: In [1]: run simple.py from_Robert.py sfepy: error: required missing: ['filename_mesh'] sfepy: left over: ['E', '__builtins__', '__file__', 'tractionLoad', '__name__', 'lamb', 'mu', 'fileName_mesh', 'v', '__doc__'] --------------------------------------------------------------------------- <type 'exceptions.ValueError'> Traceback (most recent call last) /home/ryan/siue/Research/SfePy/first_problem/simple.py in <module>() 100 101 if __name__ == '__main__': --> 102 main() 103 104 /home/ryan/siue/Research/SfePy/first_problem/simple.py in main() 85 other.extend( ['equations'] ) 86 ---> 87 conf = ProblemConf.from_file( filename_in, required, other ) 88 ## print conf 89 ## pause() /home/ryan/svn/sfepy/sfepy/base/conf.py in from_file(filename, required, other) 155 obj.__dict__.update( define_dict ) 156 --> 157 other_missing = obj.validate( required = required, other = other ) 158 for name in other_missing: 159 setattr( obj, name, None ) /home/ryan/svn/sfepy/sfepy/base/conf.py in validate(self, required, other) 218 219 if err: --> 220 raise ValueError 221 222 return other_missing <type 'exceptions.ValueError'>: WARNING: Failure executing file: <simple.py> In [2]: %debug > /home/ryan/svn/sfepy/sfepy/base/conf.py(220)validate() 219 if err: --> 220 raise ValueError 221 The files I am trying to use are attached. Thanks, Ryan

3 12

lxml dependency not documented (I think)
by Ryan Krauss July 31, 2008

July 31, 2008

I don't think the dependency on lxml is documented. I had some issues building the docs using the gen script. Installing lxml, and the dev pacakges for libxml and libxslt solved the problem. Ryan

2 3

Beginner's errors with runTests.py
by rpmu...＠gmail.com July 24, 2008

July 24, 2008

Hi everyone! I'm a quantum chemist who would like to learn about FEM analysis, and Ondrej pointed me to this package and this group. I'm running on Mac OS X, and I've built numpy 1.1.0 and scipy-0.6.0 from source. I installed the sfepy-release-00.41.03 that was available from mac.softpedia. When I try to run ./runTests.py --debug, I get the following: $ ./runTests.py --debug <<< directory: tests, test files: 18 <<< tests/test_base.py sfe: warning: other missing: ['fileName_mesh', 'field_[0-9]+|fields', 'ebc_[0-9]+|ebcs', 'fe', 'equations', 'region_[0-9]+|regions', 'variable_[0-9]+|variables', 'material_[0-9]+|materials', 'integral_[0-9]+|integrals', 'solver_[0-9]+|solvers', 'functions', 'modules', 'epbc_[0-9]+|epbcs', 'lcbc_[0-9]+|lcbcs', 'nbc_[0-9]+| nbcs', 'options'] sfe: warning: left over: ['TestCommon'] >>> test instance prepared (2 test(s)) +++ test_structAdd: ok +++ test_structIAdd: ok >>> all passed in 0.00 s <<< tests/test_elasticity_small_strain.py sfe: warning: other missing: ['equations', 'functions', 'modules', 'epbc_[0-9]+|epbcs', 'lcbc_[0-9]+|lcbcs', 'nbc_[0-9]+|nbcs', 'options'] sfe: warning: left over: ['allYourBases', 'fileName_meshes', 'getPars', 'equations_general', 'equations_iso', 'TestCommon'] >>> test instance prepared (2 test(s)) [('test_get_solution', <bound method Test.test_get_solution of Test>), ('test_linear_terms', <bound method Test.test_linear_terms of Test>)] >>> <type 'exceptions.AttributeError'> Traceback (most recent call last): File "./runTests.py", line 222, in <module> main() File "./runTests.py", line 217, in main op.walk( options.testDir, wrapRunTests( options ), stats ) File "/Library/Frameworks/Python.framework/Versions/2.5/lib/ python2.5/posixpath.py", line 290, in walk func(arg, top, names) File "./runTests.py", line 148, in runTests nFail, nTotal, testTime = runTest( confName, options ) File "./runTests.py", line 111, in runTest ok, nFail, nTotal = test.run( options.debug ) File "/Users/rmuller/Python/sfepy/sfepy-release-00.41.03/sfe/base/ testing.py", line 38, in run ret = testMethod() File "tests/test_elasticity_small_strain.py", line 156, in test_get_solution from sfe.solvers.generic import solveStationary File "/Users/rmuller/Python/sfepy/sfepy-release-00.41.03/sfe/solvers/ __init__.py", line 3, in <module> from ls import * File "/Users/rmuller/Python/sfepy/sfepy-release-00.41.03/sfe/solvers/ ls.py", line 6, in <module> import scipy.linsolve.umfpack as um AttributeError: 'module' object has no attribute 'umfpack' Trying to dig a little deeper, I see that in my version of scipy I can successfully import scipy.linsolve.umfpack: >>> import scipy.linsolve.umfpack but when I try to install umfpack *as* something, I run into problems: >>> import scipy.linsolve.umfpack as um Traceback (most recent call last): File "<stdin>", line 1, in <module> AttributeError: 'module' object has no attribute 'umfpack' Any suggestsions? Or is this a bug that should be taken up with the scipy-users list?

3 9

Re: [sage-support] Re: symbolic vector calculus
by Ondrej Certik July 21, 2008

July 21, 2008

On Mon, Jul 21, 2008 at 10:20 PM, William Stein <wst...(a)gmail.com> wrote: > > On Mon, Jul 21, 2008 at 10:17 PM, Ondrej Certik <ond...(a)certik.cz> wrote: >> >> On Mon, Jul 21, 2008 at 10:12 PM, Tim Lahey <tim....(a)gmail.com> wrote: >>> >>> >>> On Jul 21, 2008, at 3:41 PM, Gary Furnish wrote: >>> >>>> >>>> On Mon, Jul 21, 2008 at 10:00 AM, William Stein <wst...(a)gmail.com> >>>> wrote: >>>>> On Mon, Jul 21, 2008 at 6:28 PM, <jdo...(a)gmail.com> wrote: >>>>>> >>>>>> Hello all, >>>>>> >>>>>> I have a simple question about the capabilities of Sage that I have >>>>>> not been able to resolve by looking at the documentation. I often >>>>>> find myself manipulating somewhat complex functions that take vector >>>>>> arguments. I then need to derive gradients, hessians, etc. I >>>>>> need to >>>>>> do this with out knowing the dimensions of the vectors. So for >>>>>> example, what I would like to do is something like the following. >>>>>> First, I would define the function f(x)=.5 * x' * A * x + b, perhaps >>>>>> something like: >>>>>> >>>>>> A = matrix(); >>>>>> x = vector(); >>>>>> b = vector(); >>>>>> f = function( x' * A * x + b); >>>>>> >>>>>> Then, I would like Sage to do the calculus for me, something like, >>>>>> say: >>>>>> >>>>>> f.gradient() >>>>>> sage) 2*A*x >>>>>> f.hessian() >>>>>> sage) A >>>>>> >>>>>> (Of course, I wouldn't bother using a computer algebra system for >>>>>> such >>>>>> a simple function, but you get the idea.) My question is: can >>>>>> Sage do >>>>>> things like this? I would like to avoid specifying the dimensions >>>>>> of >>>>>> x, or giving the entries of A when I define the function. >>>>> >>>> >>>> My system does not currently support this, and I have no plans to >>>> implement this, but it could be done pretty easily if there was a >>>> dimensionless vectorspace and matrixspace object, but I'm not >>>> completely convinced I see the point. If one wanted to do this they >>>> could just choose an arbitrary dimension, and as long as they don't do >>>> anything that depends on the dimension it should still give the same >>>> answer (in my system, not the current one). >>> >>> The purpose for this is computations like this come up in Finite Element >>> Analysis when you are deriving the system of equations (and probably >>> elsewhere). >>> The reason for not specifying the dimensions is that you don't know the >>> length of the vectors in advance since it is dependent upon the number >>> of >>> elements you choose. They have a definite structure, so once the size is >>> specified, you can then fill in the vectors, but you don't specify the >>> size >>> in advance. The most you would do is say that A is n by n and x is n >>> by 1. >>> >>> Maxima and Maple can't do these calculations (I've tried), but >>> Mathematica >>> can. I had to go through a lot of effort to figure out an alternative >>> to use >>> in Maple for deriving Finite Element equations, although I haven't >>> been able >>> to convert it into Sage. >> >> >> Do you have some FEM code in Python? >> >> We develope sfepy (CCing sfepy-devel): >> >> http://code.google.com/p/sfepy/ >> >> and we'd like to hook symbolic capabilities to it, so that you can do >> this kind of things symbolically. We just started to use SymPy for >> checking the numerical solutions symbolically and my secret plan is to >> be able to just write an equation in SymPy or Sage, optionally specify >> some boundary conditions and get it solved using FEM. > > I just want to quickly mention in case the originally poster doesn't > know that (1) sympy is part of Sage, and (2) Ondrej is the lead > sympy developer. Anyway, I think it would be extremely cool if Sage could solve stuff using FEM. And sfepy could perfectly be hooked in, it is pure python + some C stuff (we use swig for historical reasons, but could switch to cython, or just leave the generated C file in there) + scipy + numpy. And umfpack and pysparse. But those are all either small libraries or already in Sage. It's still too early to talk about inclusion, we first need to create a spkg, test it, etc. etc., also I'd wait with this until sfepy matures a bit more, but I think it'd make Sage extremely useful for engineering applications. Imagine downloading sage on any computer, typing "sage", then one simple command to type in the equation and then get it solved automatically and correctly. Wow. Before that one would type one command to download it from optional spkg packages. I created an issue for that: http://code.google.com/p/sfepy/issues/detail?id=50 Ondrej

1 0

Re: [sage-support] Re: symbolic vector calculus
by Ondrej Certik July 21, 2008

July 21, 2008

On Mon, Jul 21, 2008 at 10:12 PM, Tim Lahey <tim....(a)gmail.com> wrote: > > > On Jul 21, 2008, at 3:41 PM, Gary Furnish wrote: > >> >> On Mon, Jul 21, 2008 at 10:00 AM, William Stein <wst...(a)gmail.com> >> wrote: >>> On Mon, Jul 21, 2008 at 6:28 PM, <jdo...(a)gmail.com> wrote: >>>> >>>> Hello all, >>>> >>>> I have a simple question about the capabilities of Sage that I have >>>> not been able to resolve by looking at the documentation. I often >>>> find myself manipulating somewhat complex functions that take vector >>>> arguments. I then need to derive gradients, hessians, etc. I >>>> need to >>>> do this with out knowing the dimensions of the vectors. So for >>>> example, what I would like to do is something like the following. >>>> First, I would define the function f(x)=.5 * x' * A * x + b, perhaps >>>> something like: >>>> >>>> A = matrix(); >>>> x = vector(); >>>> b = vector(); >>>> f = function( x' * A * x + b); >>>> >>>> Then, I would like Sage to do the calculus for me, something like, >>>> say: >>>> >>>> f.gradient() >>>> sage) 2*A*x >>>> f.hessian() >>>> sage) A >>>> >>>> (Of course, I wouldn't bother using a computer algebra system for >>>> such >>>> a simple function, but you get the idea.) My question is: can >>>> Sage do >>>> things like this? I would like to avoid specifying the dimensions >>>> of >>>> x, or giving the entries of A when I define the function. >>> >> >> My system does not currently support this, and I have no plans to >> implement this, but it could be done pretty easily if there was a >> dimensionless vectorspace and matrixspace object, but I'm not >> completely convinced I see the point. If one wanted to do this they >> could just choose an arbitrary dimension, and as long as they don't do >> anything that depends on the dimension it should still give the same >> answer (in my system, not the current one). > > The purpose for this is computations like this come up in Finite Element > Analysis when you are deriving the system of equations (and probably > elsewhere). > The reason for not specifying the dimensions is that you don't know the > length of the vectors in advance since it is dependent upon the number > of > elements you choose. They have a definite structure, so once the size is > specified, you can then fill in the vectors, but you don't specify the > size > in advance. The most you would do is say that A is n by n and x is n > by 1. > > Maxima and Maple can't do these calculations (I've tried), but > Mathematica > can. I had to go through a lot of effort to figure out an alternative > to use > in Maple for deriving Finite Element equations, although I haven't > been able > to convert it into Sage. Do you have some FEM code in Python? We develope sfepy (CCing sfepy-devel): http://code.google.com/p/sfepy/ and we'd like to hook symbolic capabilities to it, so that you can do this kind of things symbolically. We just started to use SymPy for checking the numerical solutions symbolically and my secret plan is to be able to just write an equation in SymPy or Sage, optionally specify some boundary conditions and get it solved using FEM. Do you have some pointers to the manipulations in Mathematica? Let's implement the same in Sage using the same or similar syntax. I don't think it is difficult. Ondrej

1 0