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.