fe_mknl, fe_mk

fe_mknl, fe_mk#

Purpose

Assembly of finite element model matrices.

Syntax

 [m,k,mdof] = fe_mknl(model);
 [Case,model.DOF]=fe_mknl('init',model); 
 mat=fe_mknl('assemble',model,Case,def,MatType);

Description

The exact procedure used for assembly often needs to be optimized in detail to avoid repetition of unnecessary steps. SDT typically calls an internal procedure implemented in fe_caseg Assemble and detailed in section 4.10.7. This documentation is meant for low level calls.

fe_mknl (and the obsolete fe_mk) take models and return assembled matrices and/or right hand side vectors.

Input arguments are

model a model data structure describing nodes, elements, material properties, element properties, and possibly a case.
case data structure describing loads, boundary conditions, etc. This may be stored in the model and be retrieved automatically using fe_case(model,'GetCase').
def a data structure describing the current state of the model for model/residual assembly using fe_mknl. def is expected to use model DOFs. If Case DOFs are used, they are reexpanded to model DOFs using def=struct('def',Case.T*def.def,'DOF',model.DOF). This is currently used for geometrically non-linear matrices.
MatType or Opt describing the desired output, appropriate handling of linear constraints, etc.

Output formats are

model with the additional field model.K containing the matrices. The corresponding types are stored in model.Opt(2,:). The model.DOF field is properly filled.
[m,k,mdof] returning both mass and stiffness when Opt(1)==0
[Mat,mdof] returning a matrix with type specified in Opt(1), see MatType below.

mdof is the DOF definition vector describing the DOFs of output matrices.

When fixed boundary conditions or linear constraints are considered, mdof is equal to the set of master or independent degrees of freedom Case.DOF which can also be obtained with fe_case(model,'gettdof'). Additional unused DOFs can then be eliminated unless Opt(2) is set to 1 to prevent that elimination. To prevent constraint elimination in fe_mknl use Assemble NoT.

In some cases, you may want to assemble the matrices but not go through the constraint elimination phase. This is done by setting Opt(2) to 2. mdof is then equal to model.DOF.

This is illustrated in the example below

% Low level assembly call with or without constraint resolution
 model =femesh('testubeam');
 model.DOF=[];% an non empty model.DOF would eliminate all other DOFs
 model =fe_case(model,'fixdof','Base','z==0');
 model = fe_mk(model,'Options',[0 2]); 
 [k,mdof] = fe_mk(model,'options',[0 0]); 
 fprintf('With constraints %i DOFs\n',size(k,1)); 
 fprintf('Without          %i DOFs',size(model.K{1},1));
 Case=fe_case(model,'gett');
 isequal(Case.DOF,mdof) % mdof is the same as Case.DOF

For other information on constraint handling see section 7.14.

Assembly is decomposed in two phases. The initialization prepares everything that will stay constant during a non-linear run. The assembly call performs other operations.

Init#

The fe_mknl Init phase initializes the Case.T (basis of vectors verifying linear constraints see section 7.14, resolution calls fe_case GetT, Case.GroupInfo fields (detailed below) and Case.MatGraph (preallocated sparse matrix associated with the model topology for optimized (re)assembly). Case.GroupInfo is a cell array with rows giving information about each element group in the model (see section 7.15.3 for details).

Command options are the following

NoCon Case = fe_mknl('initNoCon', model) can be used to initialize the case structure without building the matrix connectivity (sparse matrix with preallocation of all possible non zero values).
Keep can be used to prevent changing the model.DOF DOF list. This is typically used for submodel assembly.
-NodePos saves the NodePos node position index matrix for a given group in its EltConst entry.
-gstate is used force initialization of group stress entries.
new will force a reset of Case.T.

The initialization phase is decomposed into the following steps

Generation of a complete list of DOFs using the feutil('getdof',model) call.
get the material and element property tables in a robust manner (since some data can be replicated between the pl,il fields and the mat,pro stack entries. Generate node positions in a global reference frame.
For each element group, build the GroupInfo data (DOF positions).
For each element group, determine the unique pairs of [MatId ProId] values in the current group of elements and build a separate integ and constit for each pair. One then has the constitutive parameters for each type of element in the current group. pointers rows 6 and 7 give for each element the location of relevant information in the integ and constit tables.
This is typically done using an [integ,constit,ElMap]=ElemF('integinfo') command, which in most cases is really being passed directly to a p_fun('BuildConstit') command.
ElMap can be a structure with fields beginning by RunOpt_, Case_ and eval which allows execution of specific callbacks at this stage.
For each element group, perform other initializations as defined by evaluating the callback string obtained using elem('GroupInit'). For example, initialize integration rule data structures EltConst, define local bases or normal maps in InfoAtNode, allocate memory for internal state variables in gstate, ...
If requested (call without NoCon), preallocate a sparse matrix to store the assembled model. This topology assumes non zero values at all components of element matrices so that it is identical for all possible matrices and constant during non-linear iterations.

Assemble [ , NoT]#

The second phase, assembly, is optimized for speed and multiple runs (in non-linear sequences it is repeated as long as the element connectivity information does not change). In fe_mk the second phase is optimized for robustness. The following example illustrates the interest of multiple phase assembly

% Low level assembly calls
 model =femesh('test hexa8 divide 100 10 10');
 % traditional FE_MK assembly
 tic;[m1,k1,mdof] = fe_mk(model);toc
% Multi-step approach for NL operation
tic;[Case,model.DOF]=fe_mknl('init',model);toc
tic;
m=fe_mknl('assemble',model,Case,2);
k=fe_mknl('assemble',model,Case,1);
toc

MatType: matrix identifiers#

Matrix types (sometimes also noted mattyp or MatType in the documentation) are numeric indications of what needs to be computed during assembly. Currently defined types for OpenFEM are

0 mass and stiffness assembly. 1 stiffness, 2 mass, 3 viscous damping, 4 hysteretic damping
5 tangent stiffness in geometric non-linear mechanics (assumes a static state given in the call. In SDT calls (see section 4.10.7), the case entry 'curve','StaticState' is used to store the static state.
3 viscous damping. Uses info,Rayleigh case entries if defined, see example in section 5.3.2.
4 hysteretic damping. Weighs the stiffness matrices associated with each material with the associated loss factors. These are identified by the key word Eta in PropertyUnitType commands.
7 gyroscopic coupling in the body fixed frame, 70 gyroscopic coupling in the global frame. 8 centrifugal softening.
9 is reserved for non-symmetric stiffness coupling (fluid structure, contact/friction, ...);
20 to assemble a lumped mass instead of a consistent mass although using common integration rules at Gauss points.
100 volume load, 101 pressure load, 102 inertia load, 103 initial stress load. Note that some load types are only supported with the mat_og element family;
200 stress at node, 201 stress at element center, 202 stress at Gauss point
251 energy associated with matrix type 1 (stiffness), 252 energy associated with matrix type 2 (mass), ...
300 compute initial stress field associated with an initial deformation. This value is set in Case.GroupInfo{jGroup,5} directly (be careful with the fact that such direct modification INPUTS is not a MATLAB standard feature). 301 compute the stresses induced by a thermal field. For pre-stressed beams, 300 modifies InfoAtNode=Case.GroupInfo{jGroup,7}.
-1, -1.1 submodel selected by parameter, see section 4.10.7.
-2, -2.1 specific assembly of superelements with label split, see section 4.10.7.

NodePos#

NodePos=fe_mknl('NodePos',NNode,elt,cEGI,ElemF) is used to build the node position index matrix for a given group. ElemF can be omitted. NNode can be replaced by node.

nd#

nd=fe_mknl('nd',DOF); is used to build and optimized object to get indices of DOF in a DOF list.

OrientMap#

This command is used to build the InfoAtNode entry. The 'Info','EltOrient' field is a possible stack entry containing appropriate information before step 5 of the init command. The preferred mechanism is to define an material map associated to an element property as illustrated in section 7.13.

of_mk#

of_mk is the mex file supporting assembly operations. You can set the number of threads used with of_mk('setomppro',8).

obsolete#

Syntax

 model      = fe_mk(model,'Options',Opt)
 [m,k,mdof] = fe_mk( ... ,[0       OtherOptions])
 [mat,mdof] = fe_mk( ... ,[MatType OtherOptions])

fe_mk options are given by calls of the form fe_mk(model,'Options',Opt) or the obsolete fe_mk(node,elt,pl,il,[],adof,opt).

opt(1)	MatType see above
opt(2)	if active DOFs are specified using model.DOF (or the obsolete call with adof), DOFs in model.DOF but not used by the model (either linked to no element or with a zero on the matrix or both the mass and stiffness diagonals) are eliminated unless opt(2) is set to 1 (but case constraints are then still considered) or 2 (all constraints are ignored).
opt(3)	Assembly method (0 default, 1 symmetric mass and stiffness (OBSOLETE), 2 disk (to be preferred for large problems)). The disk assembly method creates temporary files using the sdtdef tempname command. This minimizes memory usage so that it should be preferred for very large models.
opt(4)	0 (default) nothing done for less than 1000 DOF method 1 otherwise. 1 DOF numbering optimized using current ofact SymRenumber method. Since new solvers renumber at factorization time this option is no longer interesting.

[m,k,mdof]=fe_mk(node,elt,pl,il) returns mass and stiffness matrices when given nodes, elements, material properties, element properties rather than the corresponding model data structure.

[mat,mdof]=fe_mk(node,elt,pl,il,[],adof,opt) lets you specify DOFs to be retained with adof (same as defining a case entry with {'KeepDof', 'Retained', adof}).

These formats are kept for backward compatibility but they do not allow support of local coordinate systems, handling of boundary conditions through cases, ...

Notes

fe_mk no longer supports complex matrix assembly in order to allow a number of speed optimization steps. You are thus expected to assemble the real and imaginary parts successively.