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