FEM_BASIS
Finite Element Basis Functions
on an M-dimensional Simplex

FEM_BASIS is a FORTRAN77 library which can evaluate basis functions associated with an M-dimensional simplex (a 1D interval, a 2D triangle, a 3D tetrahedron, and the higher-dimensional generalizations).

A complete polynomial space is generated, up to the user-specified degree. Thus, for instance, in 2D, simply by specifying the degree, the user can request a space of:

        degree 0,  1 function,  1;
        degree 1,  3 functions, 1, x, y;
        degree 2,  6 functions, 1, x, y, x^2, xy, y^2;
        degree 3, 10 functions, 1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3

and so on. Corresponding families of basis functions in higher dimensions are generated by specifying the degree.

The feature of this library is to generate a Lagrange form of the basis set of the given degree, and to evaluate those basis functions at any point.

The Triangular Case:

We will discuss the triangular case in detail; it is general enough to show all the details of the procedure, while still being simple enough to describe in pictures and tables.

Given the maximum degree D for the polynomial basis defined on a reference triangle, we have ( ( D + 1 ) * ( D + 2 ) ) / 2 monomials of degree at most D. Normally, the reference triangle is supplied with barycentric coordinates (xsi1,xsi2,xsi3) which are nonnegative and add to 1. For our purposes, we will use barycentric coordinates that have been multiplied by D. This means that all the coordinates we are interested in will be integers, which we will now identify by (I,J,K) = D * (xsi1,xsi2,xsi3). We can then define a grid of integer coordinates so that 0 <= I, J, K <= D and I+J+K = D, with (I,J,K) corresponding to

the basis point (X,Y)(I,J,K) = ( I/D, J/D );
the basis monomial P(I,J,K)(X,Y) = X^I Y^J.

For example, with D = 2, we have simply:

        F
        |\
        C-E
        |\|\
        A-B-D

with

         I J K    X    Y    P(I,J,K)(X,Y) 
 
      A (0 0 2) (0.0, 0.0)  1
      B (1 0 1) (0.5, 0.0)  x
      C (0 1 1) (0.0, 0.5)  y
      D (2 0 0) (1.0, 0.0)  x^2
      E (1 1 0) (0.5, 0.5)  x y
      F (0 2 0) (0.0, 1.0)  y^2

Note that in this arrangement, there is an obvious relationship between a particular choice of (I,J,K) and the corresponding point (X,Y), as well as a relationship between (I,J,K) and the polynomial P(I,J,K)(X,Y). However, there is no obvious connection between the point (X,Y) and the polynomial P(I,J,K)(X,Y). This will change when we move to the Lagrange basis.

Instead of the monomials P(I,J,K)(X,Y), we want a set of polynomials L(I,J,K)(X,Y) which span the same space, but have the Lagrange property, namely L(I,J,K)(X,Y) is 1 if (X,Y) is equal to (X,Y)(I,J,K), and 0 if (X,Y) is equal to any other of the basis points. One effect of this change is to tie each basis point (X,Y) to a particular basis polynomial L(I,J,K)(X,Y), namely, the polynomial which is 1 there.

This is easily arranged. Given an index (I,J,K), we compute a polynomial that has:

I factors of the form (X-0) * (X-1/D) * ... * (X-(I-1)/D);
J factors of the form (Y-0) * (Y-1/D) * ... * (Y-(J-1)/D);
K factors of the form (X+Y-(D-0)/D) * (X+Y-(D-1)/D) * ... * (X+Y-(D-(K-1))/D).

This polynomial is the product of I+J+K linear factors, in other words, it is a polynomial of degree D. This polynomial is 0 at all basis points except (X,Y)(I,J,K). If we divide this polynomial by its value at the basis point, we arrive at the desired Lagrange polynomial L(I,J,K)(X,Y).

This method is applicable for a basis set of any polynomial degree; moreover, it is a straightforward matter to extend it to higher dimensions.

The Triangular Prism:

A (right) triangular prism is a 3D shape which can be formed by drawing a triangle in the XY plane and then "lifting" the triangle along a line in the Z direction a fixed distance.

A finite element basis family can be defined for this shape, using the product of basis functions in the XY triangle and basis functions for the Z line. The basis functions are then defined by scaled barycentric coordinates I(1), I(2), I(3) for the triangle, and an independent set of scaled barycentric coordinates J(1) and J(2) for the line. A function is provided to evaluate basis functions for such an element.

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

Languages:

FEM_BASIS is available in a C version and a C++ version and a FORTRAN77 version and a FORTRAN90 version and a MATLAB version and and a Python version.

Related Data and Programs:

FEM2D_PACK, a FORTRAN77 library which contains utilities for finite element calculations.

FEM3D_PACK, a FORTRAN90 library which contains utilities for 3D finite element calculations.

Reference:

John Burkardt,
The Finite Element Basis for Simplices in Arbitrary Dimensions,
http://people.sc.fsu.edu/~jburkardt/presentations/fem_basis_2013_fsu.pdf.
Gilbert Strang, George Fix,
An Analysis of the Finite Element Method,
Cambridge, 1973,
ISBN: 096140888X,
LC: TA335.S77.
Olgierd Zienkiewicz,
The Finite Element Method,
Sixth Edition,
Butterworth-Heinemann, 2005,
ISBN: 0750663200,
TA640.2.Z54.

Source Code:

fem_basis.f, the source code.
fem_basis.sh, BASH commands to compile the source code.

Examples and Tests:

fem_basis_prb.f, a sample calling program.
fem_basis_prb.sh, BASH commands to compile and run the sample program.
fem_basis_prb_output.txt, the output file.

List of Routines:

FEM_BASIS_1D evaluates an arbitrary 1D basis function.
FEM_BASIS_2D evaluates an arbitrary triangular basis function.
FEM_BASIS_3D evaluates an arbitrary tetrahedral basis function.
FEM_BASIS_MD evaluates an arbitrary M-dimensional basis function.
FEM_BASIS_PRISM_TRIANGLE evaluates a triangular prism basis function.
R8_FRACTION uses real arithmetic on an integer ratio.
TIMESTAMP prints the current YMDHMS date as a time stamp.

You can go up one level to the FORTRAN77 source codes.

Last revised on 26 June 2013.

FEM_BASIS Finite Element Basis Functions on an M-dimensional Simplex