The Mathematica Journal
Departments
Feature Articles
Columns
New Products
New Publications
Classifieds
Calendar
News Bulletins
Mailbox
Letters
Write Us
About the Journal
Staff and Contributors
Submissions
Subscriptions
Advertising
Back Issues
Home
Download this Issue

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

[Graphics:../Images/index_gr_102.gif] [Graphics:../Images/index_gr_103.gif]
[Graphics:../Images/index_gr_104.gif] [Graphics:../Images/index_gr_105.gif]
[Graphics:../Images/index_gr_106.gif] [Graphics:../Images/index_gr_107.gif]
[Graphics:../Images/index_gr_108.gif] [Graphics:../Images/index_gr_109.gif]

The Package

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

BeginPackage["Fem`","Classes`",
"NumericalMath`GaussianQuadrature`","LinearAlgebra`MatrixManipulation`"]

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])& }
        {gradient, Module[{i,grad},
            grad = Table[(D[ # ,ref[[i]] ]& /@  functions[self])
            ,{i,1,dimension[self]} ];
            Return[Simplify[grad]] ]& }
        {physicalGradient,
                    (inverseJacobian[self].gradient[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 [Graphics:../Images/index_gr_110.gif] 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,
                     load,bcondensed,x,answer},
                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} ]];
                answer     = Table[0,{i,1,n*dof}];
                (answer[[#]]=known[[ Position[fixeddof,#][[1,1]] ]])& 
                                    /@ fixeddof;
                load = bcondensed -
                    ((Kremoved.(List/@known))[[freedof ]]);
                x = Flatten[LinearSolve[Kcondensed, load]];
                (answer[[#]] = x[[ Position[freedof,#][[1,1]] ]] )& 
                    /@ freedof;
                answer]]}

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]]& },
        {gradientSpatial, Module[{i,grad},
            grad = Table[(D[ # ,ref[[i]] ]& /@  functions[self])
            ,{i,1,dimension[self]-1} ];
            Return[Simplify[grad]] ]& } ,
        {gradientTime, Module[{i,grad},
            grad = {D[ # ,thau ]& /@  functions[self]};
            Return[Simplify[grad]] ]& } 

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 ;
        gaussquad = GaussianQuadratureWeights[npoint,-1,1,50];
        gausspointsbrick  = Partition[Flatten[Table[
            {gaussquad[[i,1]], gaussquad[[j,1]], gaussquad[[k,1]]},
        {i,1,npoint}, {j,1,npoint}, {k,1,npoint}]],3];
        gaussweightsbrick = Flatten[Table[
            gaussquad[[i,2]] gaussquad[[j,2]] gaussquad[[k,2]],
        {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]