### The Gory Details

It is time to look under the hood. The full source code is available in the file Fem.m. Interested users may consult the `::usage` messages for a detailed description of the arguments and implementation of the operations. It is necessary to introduce some of the basic terminology and the syntax of object oriented programming in Mathematica before the description of the package itself.

Object oriented programming (OOP) combines the data and functional aspects of a program into one integrated view. An OOP is described as interactions of objects in a collective manner. An object encapsulates the representation (data/instance variables) with the operations that provide its behavior. A class is the meta-description of an object. It is a template from which the instances of objects are created. Inheritance is an important aspect of OOP and is a mechanism for code reuse. New classes can be derived from existing classes. A subclass inherits all the allowed operations and representation of its superclass. It can modify and override some of the behavior of its superclass and extend its representation.

Mathematica is not an OOP language. Maeder [4] has provided a framework for Smalltalk-like [5] OOP in Mathematica, which makes the Fem package possible. It is the best of both worlds: OOP and Mathematica all in one. OOP in Mathematica is provided by the package named Classes. In Maeder's framework, a new class can be described as:

Class[class,superclass,variables,methods]

where

#### The Package

The Fem package requires the standard Mathematica packages Classes, Gaussian Quadrature, and MatrixManipulation.

BeginPackage["Fem`","Classes`",

Many FEM computations are done in the reference map of an element. A set of symbols is reserved to represent the dimensions of the reference space. The notation is generalized for arbitrary 1D, 2D, 3D, and 4D space and time coordinates. `r`, `s`, and `t` are reserved for the Cartesian reference coordinate frames corresponding to the space elements in 1D (`r`), 2D (`r,s`), and 3D (`r,s,t`). The time dimension is referred to as `thau`.

```Unprotect[ref,r,s,t,thau];
Clear[r,s,t,thau];
ref={r,s,t,thau};
Protect[ref,r,s,t,thau];
```
##### FemNode Class

The class `FemNode` represents the node in an FEM problem. The `FemNode` is a simple class that represents a global degree of freedom and a position in space.

Class[FemNode, Object,{pos, localId}, ....

The instantiation method `new` requires two arguments, `i` and `location`, which are the integer node ID and the list that represents the coordinate of the node in any dimension.

{new,
Function[{i,location},
new[super];
setCoordinate[self,location];
setId[self,i];
self]}

Here is the protocol that provides the get-set operations for the `FemNode` class.

{coordinate,    pos&}
{id,        localId&}
{setId,        Function[ {i}, localId = i] }
{setCoordinate,    Function[ {p}, pos = p] }

The protocol to support graphics operations includes two methods. `asPoint` returns the coordinate of the node as a 3D point. `displayNode` returns a graphics list that contains a `Point` and `Text` object corresponding to the location of the node and its ID.

{asPoint,
Module[ {pstn,i},
If[ Length[pos] < 3,
pstn = Table[0,{i,1,3}];
Do[ pstn[[i]] = pos[[i]], {i,1,Length[pos]}];
Return[pstn],
Return[pos]] ]&}
{displayNode,
{Point[asPoint[self]],
Text[FontForm[InputForm[id[self]],
{"Times-Roman",14}],asPoint[self]]}& }

##### FemElement Class

This is the abstract class that defines the behavior of an element. There should not be any instances of this class. The subclasses of `FemElement` define the behavior of the specific elements. The subclass of the `FemElement` must implement at least four methods.

1. functions. The method that returns the basis functions for an element in a list. Basis functions must be functions of (`r`, `s`, `t`, and `thau`).
2. gausspoints. The list of Gauss points for the element. This list is used for the Gaussian quadrature. The list contains a list of coordinates for the Gauss points. Each coordinate is a list of n numbers, where n is the dimension of the reference space for this element.
3. gaussweights. The list of Gauss weights corresponding to the Gauss points. This is a list of numbers. gausspoints and gaussweights must be present; otherwise no element function that performs `physicalIntegrate` will yield valid results.
4. dimension. The dimension of the reference space for the element.
```        Class[FemElement, Object, {nodz}, ....
```

The instantiation method `new` requires one argument `n` which is the list that contains the nodes of the element. The order of the nodes is important because it represents the topology of the element.

```        {new,
Function[ {n},
new[super];
setnodes[self,n];
self]}
```

Get-set operations for the `FemElement` class are:

```        {dimension,    0&}(*Default-Subclass must implement*)
{setnodes,    Function[{c},nodz = c]}
{nodes,        nodz&}
{numberOfNodes,    Length[functions[self]]&}
{getlocations,     (coordinate[#]& /@ nodz )&}
```

Methods for basis functions are:

```        {functions,  {1} & } (*Default-Subclass must implement*)
{shapeFunctions,  (List[#]&  /@ functions[self])& }
grad = Table[(D[ # ,ref[[i]] ]& /@  functions[self])
,{i,1,dimension[self]} ];
{interpolate, Function[{f}, Plus@@(f*functions[self])]}
```

Methods that compute the Jacobian, its inverse, and its determinant for the element are:

```        {jacobian, Module[{rdim,pdim,pts,xyzt,jac,i},
rdim = dimension[self];(*Dim of reference space*)
pts  = Transpose[getlocations[self]];
pdim = Length[pts] ; (*Dimension of physical space*)
xyz  = Table[ Plus@@(functions[self]*pts[[i]]) ,
{i,1,pdim}];
jac  = Table[ D[#, ref[[i]] ]& /@ xyzt ,
{i,1,rdim} ];
Return[jac] ] }
{inverseJacobian,    ( Inverse[jacobian[self]] )& }
{detJacobian,        ( Det[jacobian[self]] )& }
```

The following methods implement Gaussian quadrature for the element. Integration over the reference element is done by evaluating the element function at each Gauss point and returning the sum of products of these matrices with the weight and determinant of the Jacobian at those Gauss points.

```        {gausspoints,    {}& }(*Default-Subclass must implement*)
{gaussweights,    {}& }(*Default-Subclass must implement*)
{physicalIntegrate,
Function[ {integ}, Module[ {f,i,detJ,gw,gp,lcs,num} ,
f = integ[self];
gw = gaussweights[self];
num = Length[gw];
gp = gausspoints[self];
lcs = Take[{ref, dimension[self] ];
detJ = detJacobian[self];
Sum[    gw[[i]] (detJ /. Thread[lcs->gp[[i]] ])
( f /. Thread[lcs->gp[[i]] ] ),
{i,1,num}] ]] }
```

Graphics support is provided by the `display` function, which returns a list of polygons that represent the element and its nodes.

```        {display,  (Module[{l,i,polys,p},
l =  (List[displayNode[#]]& /@ nodes[self] );
polys = asPolygon[self];
p = Table[ Polygon @ (asPoint[#]& /@ polys[[i]] ))
,{i,1,Length[polys] } ];
Join[l, p] ])& }
```
##### Fem Class

The class that represents the global FEM problem and the mesh is called the `Fem` class. It holds the list of nodes and elements in the problem.

```        Class[ Fem, Object, { nodz, elmts }, .....
```

The instantiation method `new` requires two arguments, `n` and `e`, which are the list of nodes and elements respectively.

```        {new,
Function[
{n,e},
new[super];
setnodes[self,n];
setelements[self,e];
self]}
```

The protocol that provides the get-set operations for the `Fem` class is:

```        {elements,    elmts&            }
{nodes,        nodz&                }
{setnodes,    Function[{c},nodz = c]}    }
{setelements,    Function[{e},elmts = e]    }
```

The mesh of the problem can be visualized. The method `display` returns a `Graphics3D` object that can be displayed.

```        {display,  (Graphics3D[display[#]]& /@ elmts) & }
```

The element matrices and vectors describe governing equations that correspond to the local degrees of freedom in the elements. The assembly protocols provide methods to join the local element's dof matrix to the global matrix. Block assembly assumes that each degree of freedom at a node of the element is assembled as one block. The global matrix is organized similarly. The assembly function requires an argument that is an element function that returns a matrix.

```        {assembleblock,
Function[ {elementFunction , dofNode},
Module[ {nd,n,m,em,ng,gi,gj,li,lj,mg,d,i,j,k,nelz},
n       = Length[nodz]*dofNode;
nelz    = Length[elmts];
gm      = ZeroMatrix[n]; (* Global Matrix *)
{ng,mg} = Dimensions[gm]/dofNode;
Do[ (* Loop over the elements *)
em    = elementFunction [elmts[[k]]] ;
nd    = nodes[elmts[[k]]];
{n,m} = Dimensions[em]/dofNode;
Do[ (* Assemble local matrix em into gm *)
gi =  id[nd[[i]]] + ((d-1)*ng);
gj =  id[nd[[j]]] + ((d-1)*mg);
li =  i + ((d-1)*n);
lj =  j + ((d-1)*m);
gm[[ gi,gj]] +=  em[[li,lj]];
,{d,1,dofNode},   {i,1,n},   {j,1,m}] ,
{k,1,nelz}];
gm
]]}
```

`linearSolve` solves the linear system of equations and returns the answer x as a vector. In addition to solving the linear system of equations, the function can also impose the boundary conditions defined by the arguments `fixed` and `fixedvalues`, which correspond to the fixed degrees of freedom in the problem.

```        {linearSolve,
Function[ {K,b,fixed, fixedvalues},
Module[ {freedof,fixeddof, known,
dof,n, Kremoved,Kcondensed,
dof      = Length[fixed];
n        = Dimensions[K][[1]]/dof;
fixeddof = Flatten[Table[((i-1)n + #)& /@ fixed[[i]],
{i,1,dof}]];
freedof    = Complement[Range[1,dof*n],fixeddof];
known      = Flatten[fixedvalues] ;
Kremoved   = Transpose[  Transpose[K][[fixeddof]]   ];
Kcondensed = K[[freedof,freedof ]];
bcondensed = b[[freedof,{1} ]];
/@ fixeddof;
((Kremoved.(List/@known))[[freedof ]]);
(answer[[#]] = x[[ Position[freedof,#][[1,1]] ]] )&
/@ freedof;
```

#### Implementing Subclasses of FemElement

The `FemElement` class is an abstract class. It defines some behavior without defining implementation details of a real FEM element, (e.g., `methods`, `dimension`, `function`, `gausspoints`, `asPolygon`). The minimum set of methods that must be implemented by a concrete subclass (a class that can have instances) of the `FemElement` are: (1) functions, (2) gausspoints, (3) gaussweights, and (4) dimension. In addition to these methods, if the method `asPolygon` is implemented, the element can be displayed.

Additional abstract classes are also provided for convenience. These are: `SpatialElement`, `TimeElement`, and `SpaceTimeElement`. `SpatialElement` provides some additional structure for classifying elements that represent a spatial subdomain. `TimeElement` is a subdomain of time. The `SpaceTimeElement` class describes a subdomain in space-time. `SpaceTimeElement` and `TimeElements` override the methods that compute Gauss, its inverse, and so on. They also provide additional methods for gradients of the basis functions with respect to space and time coordinates.

##### SpatialElement

This is a shadow subclass of `FemElement`.

```    Class[SpatialElement, FemElement,{}, {} ];
```
##### TimeElement

The Jacobian of the `TimeElement` has the derivatives of the basis functions with respect to the time dimension only.

```    Class[TimeElement, FemElement, {}, {
{jacobian, Module[{dim,pts},
dim = dimension[self];
pts = Transpose[getlocations[self]];
jac = ZeroMatrix[dim];
jac[[dim,dim]] = (D[#, thau]& /@
(Plus@@(functions[self]*pts[[dim]])));
Return[jac] ]& }} ]
```
##### SpaceTimeElement

The Jacobian of the `SpaceTimeElement` has the derivatives of the basis functions with respect to space and time dimensions. `SpaceTimeElement` also provides methods for spatial and time derivatives of basis functions.

```    Class[SpaceTimeElement, FemElement, {}, {
{inverseJacobian, Module[{invjac,i,jac,coltime,dim,sdim},
jac     = jacobian[self];
dim     = dimension[self];
sdim    = Range[1,dim-1];
invjac  = Inverse[jac[[sdim , sdim ]] ];
coltime = -invjac . jac[[ dim ,sdim ]];
invjac  = AppendRows[invjac, Transpose[ {coltime}]];
invjac  = AppendColumns[invjac,ZeroMatrix[1,dim] ];
invjac[[dim,dim]] = 1/jac[[dim,dim]];
Simplify[invjac]]&},
{detJacobian, Module[{jac,dim,sdim},
jac = jacobian[self];
dim = dimension[self];
sdim = Range[1,dim-1];
(jac[[dim,dim]] Det[jac[[ sdim , sdim ]]])]&} ,
{detJacobianSpatial, Module[{jac,sdim},
jac = jacobian[self];
sdim =Range[1,dimension[self]-1];
Det[jac[[sdim,sdim ]]]]&} ,
{shapeFunctionsSpatial,Module[{fcs,numnodz},
numnodz = numberOfNodes[self];
fcs = functions[self];
Return[Simplify[
Partition[fcs,1][[ Range[1,Length[fcs]/2 ] ]] ]
/. thau->0]]& },
grad = Table[(D[ # ,ref[[i]] ]& /@  functions[self])
,{i,1,dimension[self]-1} ];
grad = {D[ # ,thau ]& /@  functions[self]};
```

#### Concrete Subclasses of the FemElement

A typical class hierarchy for the `FemElement` was provided earlier. The description of only one of these classes is included here for completeness. It demonstrates the stages of implementing a subclass of the `FemElement`. The 3D trilinear brick element `EightNodeBrick` with a two point Gaussian quadrature rule can be implemented as follows.

```        npoint = 2 ;
gausspointsbrick  = Partition[Flatten[Table[
{i,1,npoint}, {j,1,npoint}, {k,1,npoint}]],3];
gaussweightsbrick = Flatten[Table[
{i,1,npoint},   {j,1,npoint}, {k,1,npoint}]];
Class[ EightNodeBrick, SpatialElement,   {},{
{dimension, 3& },
{functions, {1/8(1+r)(1+s)(1+t) ,
1/8(1-r)(1+s)(1+t) ,
1/8(1-r)(1-s)(1+t) ,
1/8(1+r)(1-s)(1+t) ,
1/8(1+r)(1+s)(1-t) ,
1/8(1-r)(1+s)(1-t) ,
1/8(1-r)(1-s)(1-t) ,
1/8(1+r)(1-s)(1-t) }&},
{gausspoints, gausspointsbrick&},
{gaussweights, gaussweightsbrick&},
{asPolygon, Module[ {faces,i,j},
faces = { {6,7,3,2},
{7,8,4,3},
{7,6,5,8},
{8,5,1,4},
{5,6,2,1},
{4,1,2,3}};
Table[ (nodes[self])[[ faces[[i,j]] ]] ,
{i,1,6}, {j,1,4}]]&}
}];
```

Converted by Mathematica

[Article Index] [Prev Page][Next Page]