### 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 , its inverse by , and the Jacobian of the map is written . In fact, is a matrix for each , and so 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 , where m is the number of nodes of the element, the map is written as:

For example, is a temperature field defined as (just imagine that!), and the nodal temperatures are . Using one-dimensional element (, ), the mapping is obtained as . In the following example with , , , and , the temperature 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 and plot it in the reference element for the points .

``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 and temperatures ( and ) 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, . This function returns an element matrix of size , 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 and three-dimensional physical space .
2. `shapeFunctions`. The nodal scalar basis matrix . This function returns an matrix, where M is the number of nodal basis functions.

```
R = shapeFunctions[element]
```
3. `gradient`. The gradient of in the reference space is an matrix, where N is the dimension of the reference element and M is the number of nodal basis functions.
```

```
4. `physicalGradient`. The gradient of in the physical space is an 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 for an element. It requires a nodal function as an argument. The nodal function is either a discrete list of scalars , 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 , where is a scalar field. FEM discretization of the Laplacian operator can be written as the following inner product:

where 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

 0

where is a scalar value defined at each node

where is the error introduced by the approximation.

The standard Galerkin FEM suggests

 0 0

Applying the divergence theorem

This can be evaluated by the following integral

where 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 part of the integral. Therefore, the only part of the expression that the programmer must provide is . `physicalIntegrate` accepts this expression as an argument. The expression can be any function of the element. For example, if the expression for is

``````LaplacianIntegrand[e_] :=
Module[ {GradR } ,
];``````

then the integral 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 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]