Re: trying to understand linear elasticity

1 Aug 2008


      Thanks for the explanations.  Prototyping in Python makes a lot of
sense and is far less intimidating .
Ryan
On Fri, Aug 1, 2008 at 4:50 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
...
Ryan Krauss wrote:
...
So, I dug into sfepy/sfepy/terms/termsLinElasticity.py(61)__call__() a
bit to try and understand how linear elasticity is being handled:
Let me go line-by-line:
...
def __call__( self, diff_var = None, chunk_size = None, **kwargs ):
Pdb().set_trace()
material, virtual, state = self.get_args( **kwargs )
here we extract the arguments passed to the term.
...
    ap, vg = virtual.get_approximation( self.get_current_group(), 'Volume' )
'ap' stores various FE-related data corresponding to the given variable,
here 'virtual'. state.get_approximation(...) would return the same, as
both the test (virtual) and unknown (state) variables use the same
approximation here.
'vg' is a VolumeGeometry instance, it holds derivatives of shape (base)
functions transformed from a reference element to the real ones. It
holds also the necessary jacobians and weights required for performing a
numerical quadrature on the domain of the term.
...
    shape, mode = self.get_shape( diff_var, chunk_size, ap )
'shape' is a shape of element contributions (rezidual if mode == 0,
matrix if mode == 1). It has always the following structure:
(n_el, n_qp, n_row, n_col)
n_el - number of elements to assemble
n_qp - number of quadrature points in each element (= 1 here, as the C
function already integrates the data over elements)
n_row, n_col - actual dimensions
rezidual mode assembles D * u, matrix mode just D
...
    fargs = self.build_c_fun_args( material, state, ap, vg )
this just builds arguments passed to the C function. subclasses may use
different function, they just override build_c_fun_args.
...
    for out, chunk in self.char_fun( chunk_size, shape ):
        status = self.function( out, *fargs + (chunk, mode) )
        yield out, chunk, status
all actual computation is done here, the function is implemented in C
for speed, see sfepy/terms/extmods/termsLinearElasticity.c, function
dw_lin_elastic_iso()
C is not required though, it is an implementation detail - you can well
prototype your code in Python.
...
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?
indx - if the problem is described by more variables, then indx of each
variable point to the global state vector to data of corresponding to
the variable.
data - in time-dependent problems we usually need to remember some
previous time-steps - they are stored in a deque, which is just a list
which you can prepend and rotate.
You do not need to care much about this, just use u = state() to get to
the actual displacements.
...
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.
As I said, you can well start with just Python, speed can be obtained later.
r.
...