fem_basis, a FORTRAN90 code 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 code is to generate a Lagrange form of the basis set of the given degree, and to evaluate those basis functions at any point.
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
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:
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.
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.
The computer code and data files described and made available on this web page are distributed under the MIT license
fem_basis is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version and a Python version.
FEM2D_PACK, a FORTRAN90 code which contains utilities for finite element calculations.
FEM3D_PACK, a FORTRAN90 code which contains utilities for 3D finite element calculations.