The Mathematica Journal
Feature Articles
New Products
New Publications
News Bulletins
Write Us
About the Journal
Staff and Contributors
Back Issues
Download this Issue

Let's Do Some Element Integrals

Most functions of elements can be defined in either the physical (global) element or a reference (local) element. The change of basis transformation, the Jacobian matrix, provides the mapping between derivatives of the physical and reference element. Similarly, basis functions can be defined in physical or the reference space. I will consider basis functions in the reference space only.

Most of the entities of importance to an element are element matrix functions, which are matrix/vector-valued functions defined on the element. They can be expressed as being defined either on the reference element E or on the physical element e. I will use the following notation: the variable ranges over the reference element E, and the variable x ranges over the physical element e. The correspondence between them is denoted by [Graphics:../Images/index_gr_4.gif], its inverse by [Graphics:../Images/index_gr_5.gif], and the Jacobian of the map [Graphics:../Images/index_gr_6.gif] is written [Graphics:../Images/index_gr_7.gif]. In fact, [Graphics:../Images/index_gr_8.gif] is a matrix for each , and so [Graphics:../Images/index_gr_9.gif] is one of the matrix-valued functions of the element. In simple terms this says that, if you can express something on a nicely shaped simple element, then the Jacobian will help you translate into the very ugly-looking complex physical counterpart.

Imagine the physical variable x that is defined on the physical space. The basis functions of the element e are defined as [Graphics:../Images/index_gr_10.gif], where m is the number of nodes of the element, the map [Graphics:../Images/index_gr_11.gif] is written as:


For example, [Graphics:../Images/index_gr_13.gif] is a temperature field defined as [Graphics:../Images/index_gr_14.gif] (just imagine that!), and the nodal temperatures [Graphics:../Images/index_gr_15.gif] are [Graphics:../Images/index_gr_16.gif]. Using one-dimensional element ([Graphics:../Images/index_gr_17.gif], [Graphics:../Images/index_gr_18.gif]), the mapping [Graphics:../Images/index_gr_19.gif] is obtained as [Graphics:../Images/index_gr_20.gif]. In the following example with [Graphics:../Images/index_gr_21.gif], [Graphics:../Images/index_gr_22.gif], [Graphics:../Images/index_gr_23.gif], and [Graphics:../Images/index_gr_24.gif], the temperature [Graphics:../Images/index_gr_25.gif] varies linearly in the element and can be plotted as:


Let's define this element and its nodes using the following framework.


These statements create two nodes, node1 and node2, with the coordinates {12.5} and {30.5}, and IDs 1 and 2 respectively. The element is defined as:


This statement creates a TwoNodeLine element whose topology is defined by the ordered list of nodes {node1, node2}. The TwoNodeLine element is a subclass of the FemElement class. The next two statements create the interpolated temperature map [Graphics:../Images/index_gr_29.gif] and plot it in the reference element for the points [Graphics:../Images/index_gr_30.gif].

T = interpolate[element, {197.75, 1025.75}]
Plot[T, {r,-1,1}];


The Jacobian operator J of the element relates the physical coordinate derivatives to the reference coordinate derivatives:


Using the same one-dimensional element above, the Jacobian of the element can be computed by:

J = jacobian[element]

It is also possible to do all of these operations symbolically. Using the same problem as before, but with symbolic variables for coordinates [Graphics:../Images/index_gr_35.gif] and temperatures ([Graphics:../Images/index_gr_36.gif] and [Graphics:../Images/index_gr_37.gif]) at the nodes 1 and 2 respectively, the next example is obtained.

node1   = new[ FemNode, 1, {a}];
node2   = new[ FemNode, 2, {b}];
element = new[ TwoNodeLine, {node1,node2}];
T = interpolate[element, {Ta, Tb}]
J = jacobian[element]

A more comprehensive list of the element matrix operations for the FemElement class follows.

1. jacobian. The reference Jacobian, [Graphics:../Images/index_gr_40.gif]. This function returns an element matrix of size [Graphics:../Images/index_gr_41.gif], where M is the dimension of the reference space and N is the dimension of the physical space. For example,


    J = jacobian[element]
    JInv = inverseJacobian[element]
represents a Jacobian with a three-dimensional reference space [Graphics:../Images/index_gr_43.gif] and three-dimensional physical space [Graphics:../Images/index_gr_44.gif].
2. shapeFunctions. The nodal scalar basis matrix [Graphics:../Images/index_gr_45.gif]. This function returns an [Graphics:../Images/index_gr_46.gif] matrix, where M is the number of nodal basis functions.


    R = shapeFunctions[element]
3. gradient. The gradient of [Graphics:../Images/index_gr_48.gif] in the reference space is an [Graphics:../Images/index_gr_49.gif] matrix, where N is the dimension of the reference element and M is the number of nodal basis functions.

    GradR = gradient[element]
4. physicalGradient. The gradient of [Graphics:../Images/index_gr_51.gif] in the physical space is an [Graphics:../Images/index_gr_52.gif] matrix, where N is the dimension of the physical element and M is the number of nodal basis functions.


    GradR = physicalGradient[element]
5. interpolate. This returns the map [Graphics:../Images/index_gr_54.gif] for an element. It requires a nodal function as an argument. The nodal function is either a discrete list of scalars [Graphics:../Images/index_gr_55.gif], where M is the number of nodal basis functions, or a pure function.
interpolate[element, {a,b}]
interpolate[element, Sin[r]]
6. physicalIntegrate. In general, a system of continuum mechanics equations that form the basis of the physical phenomenon are described. These equations are then discretized using one of the standard finite element procedures such as the Galerkin formulation, least squares method, and so on.

The resulting equations, also called the governing FEM equations, are generally provided in the form of inner products in terms of other functions of the element, such as the basis functions, gradients, and so on, state variables in the problem, and material properties. The most common inner product that is utilized in FEM is the element integration.

A typical example is the Laplacian operator [Graphics:../Images/index_gr_58.gif], where is a scalar field. FEM discretization of the Laplacian operator can be written as the following inner product:


where [Graphics:../Images/index_gr_60.gif] is the finite dimensional space test function space.

With some manipulation and using the standard Galerkin FEM discretization the Laplacian operator is written as the integral

[Graphics:../Images/index_gr_61.gif] 0
[Graphics:../Images/index_gr_62.gif] [Graphics:../Images/index_gr_63.gif]

where [Graphics:../Images/index_gr_64.gif] is a scalar value defined at each node


where is the error introduced by the approximation.

The standard Galerkin FEM suggests

[Graphics:../Images/index_gr_66.gif] 0
[Graphics:../Images/index_gr_67.gif] 0

Applying the divergence theorem

This can be evaluated by the following integral


where [Graphics:../Images/index_gr_71.gif] is the subdomain of the element in the physical space.

In the reference space the integral can be evaluated as:


physicalIntegrate is a function provided by the package. It does the element integration in the reference space. The function physicalIntegrate internally performs the [Graphics:../Images/index_gr_73.gif] part of the integral. Therefore, the only part of the expression that the programmer must provide is [Graphics:../Images/index_gr_74.gif]. physicalIntegrate accepts this expression as an argument. The expression can be any function of the element. For example, if the expression for [Graphics:../Images/index_gr_75.gif] is

LaplacianIntegrand[e_] := 
  Module[ {GradR } ,
    GradR = physicalGradient[e];
    ( Transpose[GradR ] . GradR)

then the integral [Graphics:../Images/index_gr_76.gif] is

Laplacian[ e_ ] := physicalIntegrate[e, LaplacianIntegrand];

The use of the Laplacian operator, Laplacian[e_], can be demonstrated by an example with a two-dimensional element. In the following example a new element class called the LaplacianElement is created as a subclass of the FourNodeQuad element. FourNodeQuad is a subclass of the FemElement.

Class[ LaplacianElement, FourNodeQuad,

Let's define an instance of the Laplacian element by the following setup.


n1 = new[ FemNode, 1, {0,0} ];
n2 = new[ FemNode, 2, {x,0} ];
n3 = new[ FemNode, 3, {x,y} ];
n4 = new[ FemNode, 4, {0,y} ];
element = new[LaplacianElement, {n1,n2,n3,n4}];
lap = Laplacian[element];

If x and y were set to [Graphics:../Images/index_gr_78.gif] as in the following example,

        Show[Graphics3D[display[element]  /. {x->1, y->1}]];


then the Laplacian operator for the element is:

        MatrixForm[lap /. {x->1, y->1}]

Converted by Mathematica     

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