Generic compiled linear and non-linear elements

Contents Functions Index

PDF

7.16 Generic compiled linear and non-linear elements

7.16.1 Principles

To improve the ease of development of new elements, OpenFEM now supports a new category of generic element functions. Matrix assembly, stress and load assembly calls for these elements are fully standardized to allow optimization and generation of new element without recompilation. All the element specific information stored in the EltConst data structure.

Second generation volume elements are based on this principle and can be used as examples. These elements also serve as the current basis for non-linear operations.

The adopted logic is to develop families of elements with different toplogies. To implement a family, one needs

shape functions and integration rules. These are independent of the problem posed and grouped systematically in integrules.
topology, formatting, display, test, ... information for each element. This is the content of the element function (see hexa8, tetra4, ...) .
a procedure to build the constit vectors from material data. This is nominally common to all elements of a given family and is used in integinfo element call. For example p_solid('BuildConstit').
a procedure to determine constants based on current element information. This is nominally common to all elements of a given family and is used in groupinit phase (see fe_mk). The GroupInit call is expected to generate an EltConst data structure, that will be stored in the last column of Case.GroupInfo. For example hexa8 constants which calls p_solid('ConstSolid').
a procedure to build the element matrices, right hand sides, etc. based on existing information. This is compiled in of_mk MatrixIntegration and StressObserve commands. For testing/development purposes is expected that for sdtdef('diag',12) an .m file implementation in elem0.m is called instead of the compiled version.

The following sections detail the principle for linear and non-linear elements.

7.16.2 What is done in the element function

Most of the work in defining a generic element is done in the element property function (for initializations) and the compile of_mk function. You do still need to define the commands

integinfo to specify what material property function will be called to build integ, constit and elmap. For example, in hexa8, the code for this command command is
```
 if comstr(Cam,'integinfo')
  %constit integ,elmap                       ID,pl,il
  [out,out1,out2]= ...
   p_solid('buildconstit',[varargin{1};24;8],varargin{2},varargin{3});
```
input arguments passed from fe_mknl are ID a unique pair of MatId and ProId in the current element group. pl and il the material and element property fields in the model. Expected outputs are constit, integ and elmap, see Case.GroupInfo. Volume elements hexa8, q4p, ... are topology holders. They call p_solid BuildConstit which in turn calls as another property function as coded in the type (column two of il coded with fe_mat('p_fun','SI',1)). When another property funciton is called, it is expected that constit(1:2)=[-1 TypeM] to allow propagation of type information to parts of the code that will not analyze pl.
constants to specify what element property function will be called to initialize EltConst data structure and possibly set the geometry type information in pointers(4,:). For example, in hexa8, the code for this command command is
```
...
elseif comstr(Cam,'constants') 
  integ=varargin{2};constit=varargin{3};
  if nargin>3; [out,idim]=p_solid('const','hexa8',integ,constit);
  else; p_solid('constsolid','hexa8',[1 1 24 8],[]);return;
  end
  out1=varargin{1};out1(4,:)=idim; % Tell of_mk('MatrixInt') this is IDIM 
...
```
input arguments passed from fe_mknl are pointers,integ,constit the output arguments are EltConst and a modified pointers where row 4 is modified to specify a 3D underlying geometry.

If constit(1:2)=[-1 TypeM] p_solid calls the appropriate property function.

For elements that have an internal orientation (shells, beams, etc.) it is expected that orientation maps are built during this command (see beam1t, ...).
standard topology information (commands node, dof, prop, line, patch, face, edge, parent) see section 7.14.2.

hexa8 provides a clean example of what needs to be done here.

7.16.3 What is done in the element property function

*_fcn

Expected commands common to both p_* and m_* functions are the following

Subtype

With no argument returns a cell array of strings associated with each subtype (maximum is 9). With a string input, it returns the numeric value of the subtype. With a numeric input, returns the string value of the subtype. See m_elastic for the reference implementation.

database

Returns a structure with reference materials or properties of this type. Additional strings can be used to give the user more freedom to build properties.

dbval

Mostly the same as database but replaces or appends rows in model.il (for element properties) or model.pl (for material properties).

PropertyUnitType

For each subtype returns the units of each value in the property row.

p_fcn

Commands specific to p_* are associated to the implementation of a particular physical formulation.

BuidConstit

As shown in section 7.14.1 and detailed under fe_mknl the FEM initialization phase needs to resolved constitutive law information from model constants (integinfo call to the element functions) and to fill in integration constants and other initial state information (groupinit and constant calls to the element functions).

Many aspects of a finite element formulation are independent of the supporting topology. Element property functions are thus expected to deal with topology independent aspects of element constant building for a given family of elements.

Thus the element integinfo call usually just transmits arguments to a property function that does most of the work. That means defining the contents of integ and constit columns. For example for an acoustic fluid, constit colums generated by p_solid BuildConstit contain [1/

C²

].

Generic elements (hexa8, q4p, ...) all call p_solid BuildConstit. Depending on the property type coded in column 2 of the current material, p_solid attempts to call the associated m_
ti mat function with a BuildConstit command. If that fails, an attempt to call p_
ti mat is made (this allows to define a new familly of elements trough a single p_fcn p_heat is such an example).

integ nominally contains MatId,ProId,NDofPerElt,NNodePerElt,IntegRuleNumber.

Const

Similarly, element constant generation of elements that support variable integration rules is performed for an element family. For example, p_solid const supports for 3D elastic solids, for 2D elastic solids and 3D acoustic fluid volumes. p_heat supports 2D and 3D element constant building for the heat equation.

Generic elements (hexa8, q4p, ...) all use the call [EltConst,NDNDim] = p_solid('Const',ElemF, integ, constit). User extendibility requires that the user be able to bypass the normal operation of p_solid const. This can be achieved by setting constit(1)=-1 and coding a property type in the second value (for example constit(1)=fe_mat('p_heat','SI',1). The proper function is then called witht the same arguments as p_solid.

7.16.4 Generic multi-physic linear elements

This element family supports a fairly general definition of linear multi-physic elements whose element integration strategy is fully described by an EltConst data structure. hexa8 and p_solid serve as a prototype element function. Element matrix and load computations are implemented in the of_mk.c MatrixIntegration command with StrategyType=1, stress computations in the of_mk.c StressObserve command.

EltConst=hexa8('constants',[],[1 1 24 8],[]);
integrules('texstrain',EltConst)
EltConst=integrules('stressrule',EltConst);
integrules('texstress',EltConst)

Elements of this family are standard element functions (see section 7.14) and the element functions must thus return node, prop, dof, line, patch, edge, face, and parent values. The specificity is that all information needed to integrate the element is stored in an EltConst data structure that is initialized during the fe_mknl GroupInit phase.

For DOF definitions, the family uses an internal DOF sort where each field is given at all nodes sequentially 1x 2x ... 8x 1y ... 8 y ... while the more classical sort by node 1x 1y ... 2x ... is still used for external access (internal and external DOF sorting are discussed in section 7.14.2).

Each linear element matrix type is represented in the form of a sum over a set of integration points

k^(e) =

ji,jj

[{B_ji} D_{ji jk}(w(jw)) {B_jj}^T ]J(w(jw)) W((jw)) (7.3)

where the jacobian of the transformation from physical xyz to element rst coordinates is stored in EltConst.jdet(jw) and the weighting associated with the integration rule is stored in EltConst.w(jw,4).

The relation between the Case.GroupInfo constit columns and the D_ij constitutive law matrix is defined by the cell array EltConst.ConstitTopology entries. For example, the strain energy of a acoustic pressure formulation (p_solid ConstFluid) is given by

The integration rule for a given element is thus characterized by the strain observation matrix B_ji(r,s,t) which relates a given strain component

_ji and the nodal displacements. The generic linear element family assumes that the generalized strain components are linear functions of the shape functions and their derivatives in euclidian coordinates (xyz rather than rst).

The first step of the element matrix evaluation is the evaluation of the EltConst.NDN matrix whose first Nw columns store shape functions, Nw next their derivatives with respect to x, then y and z for 3D elements

[NDN]_{Nnode × Nw (Ndims + 1)} = [[N(r,s,t)][

¶ N

¶ x

] [

¶ N

¶ y

] [

¶ N

¶ z

] ] (7.4)

To improve speed the EltConst.NDN and associated EltConst.jdet fields are preallocated and reused for the assembly of element groups.

For each strain vector type, one defines an int32 matrix

EltConst.StrainDefinition{jType} with each row describing row, NDNBloc, DOF, NwStart, NwTot giving the strain component number (these can be repeated since a given strain component can combine more than one field), the block column in NDN (block 1 is N, 4 is ¶ N/¶ z), the field number, and the starting integration point associated with this strain component and the number of integration points needed to assemble the matrix. The default for NwStart NwTot is 1, Nw but this formalism allows for differentiation of the integration strategies for various fields. The figure below illustrates this construction for classical mechanical strains.

To help you check the validity of a given rule, you should fill the

EltConst.StrainLabels{jType} and EltConst.DofLabels fields and use the
integrules( 'texstrain', EltConst) command to generate a LATEX printout of the rule you just generated.

The .StrainDefinition and .ConstitTopology information is combined automatically in integrules to generate .MatrixIntegration (integrules MatrixRule command) and .StressRule fields (integrules StressRule command). These tables once filed properly allow an automated integration of the element level matrix and stress computations in OpenFEM.

7.16.5 Generic RHS computations

Right hand side (load) computations can either be performed once (fixed set of loads) through fe_load which deals with multiple loads, or during an iterative process where a single RHS is assembled by fe_mknl into the second column of the state argument dc.def(:,2) along with the matrices when requiring the stiffness with MatDes=1 or MatDes=5 (in the second case, the forces are assumed following if implemented).

There are many classical forms of RHS, one thus lists here forms that are implemented in of_mk.c MatrixIntegration. Computations of these rules, requires that the EltConst.VectMap field by defined. Each row of EltConst.RhsDefinition specifies the procedure to be used for integration.

Two main strategies are supported where the fields needed for the integration of loads are stored either as colums of dc.def (for fields that can defined on DOFs of the model) or as entries in InfoAtNode (Case.GroupInfo{7}). In the later case, each column of InfoAtNode specifies an input field specified at nodes of the model or with a NodePos field at specific nodes. The new strategy is compatible with both continuous and discontinous fields at each node. THIS HAS BEEN REVISED FOR 2007a and is still unstable.

Initialization of InfoAtNode is performed with fe_mknl('Init -gstate') calls.

Currently the only accepted format for rows of EltConst.RhsDefinition is

101(1) InfoAtNode1(2) InStep(3) NDNOff1(4) FDof1(5) NDNCol(6)
NormalComp(7) w1(8) nwStep(9)

Where InfoAtNode1 gives the first row index in storing the field to be integrated in InfoAtNode. InStep gives the index step (3 for a 3 dimensional vector field), NDNOff1 gives the block offset in the NDN matrix (zero for the nominal shape function). FDof1 gives the offset in force DOFs for the current integration. NDNCol. If larger than -1, the normal component NormalComp designs a row number in EltConst.bas, which is used as a weighting coefficient. tt w1 gives the index of the first gauss point to be used (in C order starting at 0). nwStep gives the number of gauss points in the rule being used.

volume forces not proportional to density

ó
õ

W₀
f_v(x) . du(x) = {F_v}_k =

j_w
( {N_k(j_w)} {N_j(j_w)} f_v(x_j) ) J(j_w) W(j_w) (7.5)

are thus described by
```
 opt.RhsDefinition=int32( ...
   [101 0 3 0     0 0 -1    rule+[-1 0];
    101 1 3 0     1 0 -1    rule+[-1 0];
    101 2 3 0     2 0 -1    rule+[-1 0]]);
```
for 3D solids (see p_solid).

Similarly, normal pressure is integrated as 3 volume forces over 3D surface elements with normal component weighting

F_m =

ó
õ

¶ W₀
p(x) n_m(x) . dv(x)

=

j_w
( {N_k(j_w)} {N_j(j_w)} p(x_j) n_m ) J(j_w) W(j_w)

(7.6)
inertia forces (volume forces proportional to density)

F= ó
õ

W₀
(x) f_v(x) . dv(x) (7.7)
stress forces (will be documented later)

7.16.6 Non-linear iterations, what is done in of_mk

Non linear problems are characterized by the need to perform iterations with multiple assemblies of matrices and right hand sides (RHS). To optimize the performance, the nominal strategy for non-linear operations is to

perform an initialization (standard of_mknl init call)
define a deformation data structure dc with two columns giving respectively the current state and the non linear RHS.

At a given iteration, one resets the RHS and performs a single fe_mknl call that returns the current non-linear matrix and replaces the RHS by its current value (note that fe_mknl actually modifies the input argument dc which is not an normal MATLAB behaviour but is needed here for performance)

 % at init allocate DC structure
 dc=struct('DOF',model.DOF,'def',zeros(length(model.DOF),2); 
 % ... some NL iteration mechanism here
 dc.def(:,2)=0; % reset RHS at each iteration
 k=fe_mknl('assemble not',model,Case,dc,5); % assemble K and RHS

Most of the work for generic elements is done within the of_mk MatrixIntegration command that is called by fe_mknl. Each call to the command performs matrix and RHS assembly for a full group of elements. Three strategies are currently implemented

Linear multiphysic elements of arbitrary forms, see section 7.16.4
Elastic3DNL fully anisotropic elastic elements in geometrically non-linear mechanics problems. Element matrix are implemented in the of_mk.c MatrixIntegration command with StrategyType=2 for the linear tangent matrix (MatType=5). Other computations are performed using generic elements (section 7.16.4) (mass MatType=2). This formulation family has been tested for the prediction of vibration responses under static pre-load.

Stress post-processing is implemented using the underlying linear element.
Hyperelastic elements for large transformation problems. These elements have been tested through the RivlinCube example.

Simultaneous element matrix and right hand side computations are implemented in the of_mk.c MatrixIntegration command with StrategyType=3 for the linear tangent matrix (MatType=5). In this case (and only this case!!), the EltConst.NDN matrix is built as follow:
for 1£ jw £ Nw
[NDN]_{(Ndims + 1) × Nnode (Nw)} = [[NDN]^jw] (7.8)
with

[NDN]_{(Ndims + 1) × Nnode}^jw = [

[N(r,s,t)]_jw

[

¶ N

¶ x
]_jw

[

¶ N

¶ y
]_jw

[

¶ N

¶ z
]_jw

] (7.9)
This implementation corresponds to case 31 of NDNSwitch function in of_mk_pre.c. The purpose is to use C-BLAS functions in element matrix and right hand side computations implemented in the same file (function Mecha3DintegH) to improve speed.

Other computations are performed using generic elements (section 7.16.4) (mass MatType=2). This formulation family has been tested for the RivlinCube test.

Stress post-processing is not implemented.