Creating new elements (advanced tutorial)

Contents Functions Index

PDF

7.14 Creating new elements (advanced tutorial)

In this section one describes the developments needed to integrate a new element function into OpenFEM. First, general information about OpenFEM work is given. Then the writing of a new element function is described. And at last, conventions which must be respected are given.

7.14.1 Conventions

Element functions and C functionality

In OpenFEM, elements are defined by element functions. Element functions provide different pieces of information like geometry, degrees of freedom, model matrices, ...

OpenFEM functions like the preprocessor femesh, the model assembler fe_mk or the post-processor feplot call element functions for data about elements.

For example, in the assembly step, fe_mk analyzes all the groups of elements. For each group, fe_mk gets its element type (bar1, hexa8, ...) and then calls the associated element function.
First of all, fe_mk calls the element function to know what is the rigth call form to compute the elementary matrices (eCall=elem0('matcall') or eCall=elem0('call'), see section 7.14.2 for details). eCall is a string. Generally, eCall is a call to the element function. Then for each element, fe_mk executes eCall in order to compute the elementary matrices.

This automated work asks for a likeness of the element functions, in particular for the calls and the outputs of these functions. Next section gives information about element function writing.

Standard names in assembly routines

cEGI	vector of element property row indices of the current element group (without the group header)
constit	real (double) valued constitutive information. The constit for each group is stored in Case.GroupInfo{jGroup,4};.
def.def	vector of deformation at DOFs. This is used for non-linear, stress or energy computation calls that need displacement information.
EGID	Element Group Identifier of the current element group (different from jGroup if an EGID is declared).
elt	model description matrix. The element property row of the current element is given by elt(cEGI(jElt),:) which should appear in the calling format eCall of your element function.
ElemF	name of element function or name of superelement
ElemP	parent name (used by femesh in particular to allow property inheritance)
gstate	real (double) valued element state information.
integ	int32 valued constitutive information.
jElt	number of the current element in cEGI
jGroup	number of the current element group (order in the element matrix). [EGroup,nGroup]=getegroup(elt); finds the number of groups and group start indices.
nodeE	nodes of the current element.
NNode	node identification reindexing vector. NNode(ID) gives the row index (in the node matrix) of the nodes with identification numbers ID. You may use this to extract nodes in the node matrix using something like node(NNode(elt(cEGI(jElt),[1 2])),:) which will extract the two nodes with numbers given in columns 1 and 2 of the current element row (an error occurs if one of those nodes is not in node). This can be built using NNode=sparse(node(:,1),1,1:size(node,1).
pointers	one column per element in the current group gives.

Case.GroupInfo cell array

The meaning of the columns is as follows

DofPos Pointers Integ Constit gstate ElMap InfoAtNode EltConst

DofPos	int32 matrix whose columns give the DOF positions in the full matrix of the associated elements. Numbering is C style (starting at 0) and -1 is used to indicate a fixed DOF.
pointers	int32 matrix whose columns describe information each element of the group. Pointers has one column per element giving [OutSize1 OutSize2 u3 NdNRule MatDes IntegOffset ConstitOffset StateOffset] Outsize1 size of element matrix (for elements issued from MODULEF), zero otherwise. MatDes type of desired output. See the fe_mk MatType section for a current list. IntegOffset gives the starting index (first element is 0) of integer options for the current element in integ. ConstitOffset gives the starting index (first element is 0) of real options for the current element in constit.
integ	int32 matrix storing integer values used to describe the element formulation of the group. Meaning depends on element family and should be documented in the element property function (p_solid BuildConstit for example). The nominal content of an integ column (as return by the element integinfo call) is MatId,ProId,NDofPerElt,NNodePerElt,IntegRuleType where integrules(ElemP,IntegRuleType) is supposed to return the approriate integration rule.
constit	double matrix storing integer values used to describe the element formulation of the group. Meaning depends on element family and should be documented in the element property function (p_solid BuildConstit for example).
gstate	double matrix whose columns describe the internal state of each element of the group. By default, column content is stress at integration points (Nstrain× Nw values). Users are of course free to add any appropriate value for their own elements, a typical application is the storage of internal variables. For an example of gstate initialization see fe_stres thermal.
ElMap	int32 element map matrix used to distinguish between internal and external element DOF numbering (for example : hexa8 uses all x DOF, then all y ... as internal numbering while the external numbering is done using all DOFs at node 1, then node 2, ...). The element matrix in extrenal sort is given by k_ext=ke(ElMap). EltConst.VectMap gives similar reordering information for vectors (loads, ...).

InfoAtNode	double matrix whose rows describe information at element nodes (as many columns as nodes in the model). An other possible format is a structure with .NodePos (int32) with as many columns as elements in the group giving column positions in a .data field.
EltConst	struct used to store element formulation information (integration rule, constitutive matrix topology, etc.) Details on this data structure are given in section 7.14.1.

Element constant data structure

The EltConst data structure is used in most newer generation elements implemented in of_mk.c. It contains geometric and integration rule properties.

.N	nw × Nnode shape functions at integration points
.Nr	nw × Nnode derivative of shape function with respect to the first reference coordinate r
.Ns	nw × Nnode derivative of shape function with respect to the second reference coordinate s
.Nt	nw × Nnode derivative of shape function with respect to the second reference coordinate t
.NDN	Nshape × nw (1+Ndim) memory allocation to store the shape functions and their derivatives with respect to physical coordinates [N N,x N,y N,z]. of_mk currently supports the following geometry rules 3 3D volume, 2 2D volume, 23 3D surface, 13 3D line (see integrules BuildNDN for calling formats). Cylindrical and spherical coordinates are not currently supported. In the case of rule 31 (hyperelastic elements), the storage scheme is modifified to be (1+Ndim) × Nshape × nw which perserves data locality better.
.jdet	Nw memory allocation to store the determinant of the jacobian matrix at integration points.
.bas	9× Nw memory allocation to store local material basis. This is in particular used for 3D surface rules where components 6:9 of each column give the normal.
.Nw	number of integration points (equal to size(EltConst.N,1))
.Nnode	number of nodes (equal to size(EltConst.N,2)=size(EltConst.NDN,1))
.xi	Nnode × 3 reference vertex coordinates
.VectMap	index vector giving DOF positions in external sort. This is needed for RHS computations.

7.14.2 Basic commands of an element function

The first step to create a new element is to write a new element function.

In Matlab version, a typical element function is an .m or .mex file that is in your MATLAB path. In Scilab version, a typical element function is an .sci or mex file that is loaded into Scilab memory (see getf in Scilab on-line help).

The name of the function/file corresponds to the name of the element (thus the element bar1 is implemented through the bar1.m file)

General element information

To build a new element take q4p.m or q4p.sci as an example.

As for all Matlab or Scilab functions, the header is composed of a function syntax declaration and a help section. The following example is written for Matlab. For Scilab version, don't forget to replace % by //. In this example, the name of the created element is elem0.

For element functions the nominal format is

 function [out,out1,out2]=elem0(CAM,varargin);
 %elem0 help section

The element function should then contain a section for standard calls which let other functions know how the element behaves.

 if isstr(CAM) %standard calls with a string command

  [CAM,Cam]=comstr(CAM,1); % remove blanks
  if comstr(Cam,'integinfo')
   % some code needed here
   out= constit; % real parameter describing the constitutive law
   out1=integ;   % integer (int32) parameters for the element
   out2=elmap;   

  elseif comstr(Cam,'matcall')
   out=elem0('call');
   out1=1; % SymFlag
  elseif comstr(Cam,'call');     out = ['AssemblyCall'];
  elseif comstr(Cam,'rhscall');  out = ['RightHandSideCall'];
  elseif  comstr(Cam,'scall');   out = ['StressComputationCall'];
  elseif  comstr(Cam,'node');    out = [NodeIndices]; 
  elseif  comstr(Cam,'prop');    out = [PropertyIndices]; 
  elseif  comstr(Cam,'dof');     out = [ GenericDOF ];
  elseif comstr(Cam,'patch');    
                         out = [ GenericPatchMatrixForPlotting ];
  elseif comstr(Cam,'edge');     out = [ GenericEdgeMatrix ];
  elseif comstr(Cam,'face');     out = [ GenericFaceMatrix ];
  elseif comstr(Cam,'sci_face'); out = [ SciFaceMatrix ];
  elseif comstr(Cam,'parent');   out = ['ParentName'];
  elseif comstr(Cam,'test')
     % typically one will place here a series of basic tests
  end
  return
 end % of standard calls with string command

The expected outputs to these calls are detailed below.

call,matcall

Format string for element matrix computation call. Element functions must be able to give fe_mk the proper format to call them (note that superelements take precedence over element functions with the same name, so avoid calling a superelement beam1, etc.).

matcall is similar to call but used by fe_mknl. Some elements directly call the of_mk mex function thus avoiding significant loss of time in the element function. If your element is not directly supported by fe_mknl use matcall=elem0('call').

The format of the call is left to the user and determined by fe_mk by executing the command eCall=elem0('call'). The default for the string eCall should be (see any of the existing element functions for an example)

  [k1,m1]=elem0(nodeE,elt(cEGI(jElt),:),...
                pointers(:,jElt),integ,constit,elmap);

To define other proper calling formats, you need to use the names of a number of variables that are internal to fe_mk. fe_mk variables used as output arguments of element functions are

k1	element matrix (must always be returned, for opt(1)==0 it should be the stiffness, otherwise it is expected to be the type of matrix given by opt(1))
m1	element mass matrix (optional, returned for opt(1)==0, see below)

[ElemF,opt,ElemP]=feutil('getelemf',elt(EGroup(jGroup),:),jGroup)
returns, for a given header row, the element function name ElemF, options opt, and parent name ElemP.

fe_mk and fe_mknl variables that can be used as input arguments to element function are listed in section 7.14.1.

dof, dofcall

Generic DOF definition vector. For user defined elements, the vector returned by elem0('dof') follows the usual DOF definition vector format (NodeId.DofId or -1.DofId) but is generic in the sense that node numbers indicate positions in the element row (rather than actual node numbers) and -1 replaces the element identifier (if applicable).

For example the bar1 element uses the 3 translations at 2 nodes whose number are given in position 1 and 2 of the element row. The generic DOF definition vector is thus [1.01;1.02;1.03;2.01;2.01;2.03].

A dofcall command may be defined to bypass generic dof calls. In particular, this is used to implement elements where the number of DOFs depends on the element properties. The command should always return out=elem0('dofcall');. The actual DOF building call is performed in p_solid('BuildDof') which will call user p_*.m functions if needed.

Elements may use different DOF sorting for their internal computations.

edge,face,patch,line,sci_face

face is a matrix where each row describes the positions in the element row of nodes of the oriented face of a volume (conventions for the orientation are described under integrules). If some faces have fewer nodes, the last node should be repeated as needed. feutil can consider face sets with orientation conventions from other software.

edge is a matrix where each row describes the node positions of the oriented edge of a volume or a surface. If some edges have fewer nodes, the last node should be repeated as needed.

line (obsolete) is a vector describes the way the element will be displayed in the line mode (wire frame). The vector is generic in the sense that node numbers represent positions in the element row rather than actual node numbers. Zeros can be used to create a discontinuous line. line is now typically generated using information provided by patch.

patch. In MATLAB version, surface representations of elements are based on the use of MATLAB patch objects. Each row of the generic patch matrix gives the indices nodes. These are generic in the sense that node numbers represent positions in the element row rather than actual node numbers.

For example the tetra4 solid element has four nodes in positions 1:4. Its generic patch matrix is [1 2 3;2 3 4;3 4 1;4 1 2]. Note that you should not skip nodes but simply repeat some of them if various faces have different node counts.

sci_face is the equivalent of patch for use in the SCILAB implementation of OpenFEM. The difference between patch and sci_face is that, in SCILAB, a face must be described with 3 or 4 nodes. That means that, for a two nodes element, the last node must be repeated (in generallity, sci_face = [1 2 2];). For a more than 4 nodes per face element, faces must be cut in subfaces. The most important thing is to not create new nodes by the cutting of a face and to use all nodes. For example, 9 nodes quadrilateral can be cut as follows :

Figure 7.1: Lower order patch representation of a 9 node quadrilateral

but a 8 nodes quadrilaterals cannot by cut by this way. It can be cut as follows:

Figure 7.2: Lower order patch representation of a 8 node quadrilateral

integinfo, BuildConstit

[constit,integ,elmap]=elem0('integinfo',[MatId ProId],pl,il,model,Case) is supposed to search pl and il for rows corresponding to MatId and ProId and return a real vector constit describing the element consitutive law and an integer vector integ.

ElMap is used to build the full matrix of an element which initially only gives it lower or upper triangular part. If a structure is return, fe_mknl can do some group wise processing (typically initialization of internal states).

In most elements, one uses [constit,integ,elmap]=p_solid('buildconstit', [varargin{1};Ndof;Nnode],varargin{2:end}) since p_solid passes calls to other element property functions when needed.

node

Vector of indices giving the position of nodes numbers in the element row. In general this vector should be [1:n] where n is the number of nodes used by the element.

prop

Vector of indices giving the position of MatId, ProId and EltId in the element row. In general this vector should be n+[1 2 3] where n is the number of nodes used by the element. If the element does not use any of these identifiers the index value should be zero (but this is poor practice).

parent

Parent element name. If your element is similar to a standard element (beam1, tria3, quad4, hexa8, etc.), declaring a parent allows the inheritance of properties. In particular you will be able to use functions, such as fe_load or parts of femesh, which only recognize standard elements.

rhscall

rhscall is a string that will be evaluated by fe_load when computing right hand side loads (volume and surface loads). Like call or matcall, the format of the call is determined by fe_load by executing the command eCall=elem0('call'). The default for the string eCall should be :

be=elem0(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),...
                        integ,constit,elmap,estate);

The output argument be is the right hand side load. The inputs arguments are the same as those for matcall and call.

Matrix, load and stress computations

The calls with one input are followed by a section on element matrix assembly. For these calls the element function is expected to return an element DOF definition vector idof and an element matrix k. The type of this matrix is given in opt(1). If opt(1)==0, both a stiffness k and a mass matrix m should be returned. See the fe_mk MatType section for a current list.

Take a look at bar1 which is a very simple example of element function.

A typical element assembly section is as follows :

 % elem0 matrix assembly section

 % figure out what the input arguments are
 node=CAM;   elt=varargin{1}; 
 point=varargin{2};  integ=varargin{3};
 constit=varargin{4}; elmap=varargin{5};
 typ=point(5);
   
 % outputs are [k,m] for opt(1)==0
 %             [mat] for other opt(1)
 switch point(5)
 case 0
  [out,out1] = ... % place stiffness in out and mass in out1
 case 1
   out= ...  % compute stiffness
 case 2
   out= ...  % compute mass
 case 100 
   out= ...  % compute right hand side
 case 200 
   out= ...  % compute stress  ...
 otherwise
   error('Not a supported matrix type');
 end

Distributed load computations (surface and volume) are handled by fe_load. Stress computations are handled by fe_stres.

There is currently no automated mechanism to allow users to integrate such computations for their own elements without modifying fe_load and fe_stres, but this will appear later since it is an obvious maintenance requirement.

The mechanism that will be used will be similar to that used for matrix assembly. The element function will be required to provide calling formats when called with elem0('fsurf') for surface loads, elem0('fvol') for volume loads, and
elem0('stress') for stresses. fe_load and fe_stres will then evaluate thes calls for each element.

©1991-2007 by SDTools