User manual¶
Note
This page is work in progress and needs substantial editing.
FIAT (FInite element Automatic Tabulator) is a Python package for defining and evaluating a wide range of different finite element basis functions for numerical partial differential equations. It is intended to make ``difficult’’ elements such as high-order Brezzi-Douglas-Marini [BDM] elements usable by providing abstractions so that they may be implemented succinctly and hence treated as a black box. FIAT is intended for use at two different levels. For one, it is designed to provide a standard API for finite element bases so that programmers may use whatever elements they need in their code. At a lower level, it provides necessary infrastructure to rapidly deploy new kinds of finite elements without expensive symbolic computation or tedious algebraic manipulation. It is my goal that a large number of people use FIAT without ever knowing it. Thanks to several ongoing projects such as Sundance [Sundance], FFC [FFC], and PETSc [PETSc], it is becoming possible to to define finite element methods using mathematical notation in some high-level or domain-specific language. The primary shortcoming of these projects is their lack of support for general elements. It is one thing to ``provide hooks’’ for general elements, but absent a tool such as FIAT, these hooks remain mainly empty. As these projects mature, I hope to expose users of the finite element method to the exotic world of potentially high-degree finite element on unstructured grids using the best elements in \(H^1\), \(H(\mathrm{div})\), and \(H(\mathrm{curl})\).
In this brief (and still developing) guide, I will first present the high-level API for users who wish to instantiate a finite element on a reference domain and evaluate its basis functions and derivatives at some quadrature points. Then, I will explain some of the underlying infrastructure so as to demonstrate how to add new elements.
Using FIAT: A tutorial with Lagrange elements¶
Importing FIAT¶
FIAT is organized as a package in Python, consisting of several modules. In order to get some of the packages, we use the line
from FIAT import Lagrange, quadrature, shapes
This loads several modules for the Lagrange elements, quadrature rules, and the simplicial element shapes which FIAT implements. The roles each of these plays will become clear shortly.
Important note¶
Throughout, FIAT defines the reference elements based on the interval \((-1,1)\) rather than the more common \((0,1)\). So, the one-dimensional reference element is \((-1,1)\), the three vertices of the reference triangle are \((-1,-1),(1,-1),(1,-1)\), and the four vertices of the reference tetrahedron are \((-1,-1,-1),(1,-1,-1),(-1,1,-1),(-1,-1,1)\).
Instantiating elements¶
FIAT uses a lightweight object-oriented infrastructure to define
finite elements. The Lagrange
module contains a class
Lagrange
modeling the Lagrange finite element family. This
class is a subclass of some FiniteElement
class contained in
another module (polynomial
to be precise). So, having imported
the Lagrange
module, we can create the Lagrange element of
degree 2
on triangles by
shape = shapes.TRIANGLE
degree = 2
U = Lagrange.Lagrange( shape , degree )
Here, shapes.TRIANGLE
is an integer code indicating the two
dimensional simplex. shapes
also defines
LINE
and TETRAHEDRON
. Most of the
upper-level interface to FIAT is dimensionally abstracted over element
shape.
The class FiniteElement
supports three methods, modeled on the
abstract definition of Ciarlet. These methods are
domain_shape()
, function_space()
, and dual_basis()
.
The first of these returns the code for the shape and the second
returns the nodes of the finite element (including information related
to topological association of nodes with mesh entities, needed for
creating degree of freedom orderings).
Quadrature rules¶
FIAT implements arbitrary-order collapsed quadrature, as discussed in
Karniadakis and Sherwin~cite{}, for the simplex of dimension one,
two, or three. The simplest way to get a quadrature rule is through
the function `make_quadrature(shape,m)`
, which takes a shape code
and an integer indicating the number of points per direction. For
building element matrices using quadratics, we will typically need a
second or third order integration rule, so we can get such a rule by
>>> Q = quadrature.make_quadrature( shape , 2 )
This uses two points in each direction on the reference square, then
maps them to the reference triangle. We may get a
Numeric.array
of the quadrature weights with the method
Q.get_weights()
and a list of tuples storing the quadrature
points with the method Q.get_points()
.
Tabulation¶
FIAT provides functions for tabulating the element basis functions and
their derivatives. To get the FunctionSpace
object, we do
>>> Ufs = U.function_space()
To get the values of each basis function at each of the quadrature
points, we use the tabulate()
method
>>> Ufs.tabulate( Q.get_points() )
array([[ 0.22176167, -0.12319761, -0.11479229, -0.06377178],
[-0.11479229, -0.06377178, 0.22176167, -0.12319761],
[-0.10696938, 0.18696938, -0.10696938, 0.18696938],
[ 0.11074286, 0.19356495, 0.41329796, 0.72239423],
[ 0.41329796, 0.72239423, 0.11074286, 0.19356495],
[ 0.47595918, 0.08404082, 0.47595918, 0.08404082]])
This returns a two-dimensional Numeric.array
with rows for each
basis function and columns for each input point.
Also, finite element codes require tabulation of the basis functions’
derivatives. Each FunctionSpace
object also provides a method
tabulate_jet(i,xs)
that returns a list of Python dictionaries.
The i``th entry of the list is a dictionary storing the values of
all ``i``th order derivatives. Each dictionary maps a multiindex
(a tuple of length ``i
) to the table of the associated partial
derivatives of the basis functions at those points. For example,
>>> Ufs_jet = Ufs.tabulate_jet( 1 , Q.get_points() )
tabulates the zeroth and first partial derivatives of the function space at the quadrature points. Then,
>>> Ufs_jet[0]
{(0, 0): array([[ 0.22176167, -0.12319761, -0.11479229, -0.06377178],
[-0.11479229, -0.06377178, 0.22176167, -0.12319761],
[-0.10696938, 0.18696938, -0.10696938, 0.18696938],
[ 0.11074286, 0.19356495, 0.41329796, 0.72239423],
[ 0.41329796, 0.72239423, 0.11074286, 0.19356495],
[ 0.47595918, 0.08404082, 0.47595918, 0.08404082]])}
gives us a dictionary mapping the only zeroth-order partial derivative to the values of the basis functions at the quadrature points. More interestingly, we may get the first derivatives in the x- and y- directions with
>>> Ufs_jet[1][(1,0)]
array([[-0.83278049, -0.06003983, 0.14288254, 0.34993778],
[-0.14288254, -0.34993778, 0.83278049, 0.06003983],
[ 0. , 0. , 0. , 0. ],
[ 0.31010205, 1.28989795, 0.31010205, 1.28989795],
[-0.31010205, -1.28989795, -0.31010205, -1.28989795],
[ 0.97566304, 0.40997761, -0.97566304, -0.40997761]])
>>> Ufs_jet[1][(0,1)]
array([[ -8.32780492e-01, -6.00398310e-02, 1.42882543e-01, 3.49937780e-01],
[ 7.39494156e-17, 4.29608279e-17, 4.38075188e-17, 7.47961065e-17],
[ -1.89897949e-01, 7.89897949e-01, -1.89897949e-01, 7.89897949e-01],
[ 3.57117457e-01, 1.50062220e-01, 1.33278049e+00, 5.60039831e-01],
[ 1.02267844e+00, -7.29858118e-01, 4.70154051e-02, -1.13983573e+00],
[ -3.57117457e-01, -1.50062220e-01, -1.33278049e+00, -5.60039831e-01]])
Lower-level API¶
Not only does FIAT provide a high-level library interface for users to evaluate existing finite element bases, but it also provides lower-level tools. Here, we survey these tools module-by-module.
shapes.py¶
FIAT currenly only supports simplicial reference elements, but does so in a fairly dimensionally-independent way (up to tetrahedra).
jacobi.py¶
This is a low-level module that tabulates the Jacobi polynomials and their derivatives, and also provides Gauss-Jacobi points. This module will seldom if ever be imported directly by users. For more information, consult the documentation strings and source code.
expansions.py¶
FIAT relies on orthonormal polynomial bases. These are constructed by
mapping appropriate Jacobi polynomials from the reference cube to the
reference simplex, as described in the reference of Karniadakis and
Sherwin~cite{}. The module expansions.py
implements these
orthonormal expansions. This is also a low-level module that will
infrequently be used directly, but it forms the backbone for the
module polynomial.py
.
quadrature.py¶
FIAT makes heavy use of numerical quadrature, both internally and in the user interface. Internally, many function spaces or degrees of freedom are defined in terms of integral quantities having certain behavior. Keeping with the theme of arbitrary order approximations, FIAT provides arbitrary order quadrature rules on the reference simplices. These are constructed by mapping Gauss-Jacobi rules from the reference cube. While these rules are suboptimal in terms of order of accuracy achieved for a given number of points, they may be generated mechanically in a simpler way than symmetric quadrature rules. In the future, we hope to have the best symmetric existing rules integrated into FIAT.
Unless one is modifying the quadrature rules available, all of the
functionality of quadrature.py
may be accessed through the
single function make_quadrature
.
This function takes the code for a shape and the number of points in
each coordinate direction and returns a quadrature rule. Internally,
there is a lightweight class hierarchy rooted at an abstract
QuadratureRule
class, where the quadrature rules for
different shapes are actually different classes. However, the dynamic
typing of Python relieves the user from these considerations. The
interface to an instance consists in the following methods.
get_points()
, which returns a list of the quadrature points, each stored as a tuple. For dimensional uniformity, one-dimensional quadrature rules are stored as lists of 1-tuples rather than as lists of numbers.get_weights()
, which returns aNumeric.array
of quadrature weights.integrate(f)
, which takes a callable objectf
and returns the (approximate) integral over the domain- Also, the
__call__
method is overloaded so that a quadrature rule may be applied to a callable object. This is syntactic sugar on top of theintegrate
method.
polynomial.py¶
The polynomial
module provides the bulk of the classes
needed to represent polynomial bases and finite element spaces.
The class PolynomialBase
provides a high-level access to
the orthonormal expansion bases; it is typically not instantiated
directly in an application, but all other kinds of polynomial bases
are constructed as linear combinations of the members of a
PolynomialBase
instance. The module provides classes for
scalar and vector-valued polynomial sets, as well as an interface to individual
polynomials and finite element spaces.
PolynomialBase¶
PolynomialSet¶
The PolynomialSet
function is a factory function interface into
the hierarchy
[BDM] | Brezzi, Franco; Douglas, Jim, Jr.; Marini, L. D. “Two families of mixed finite elements for second order elliptic problems”. Numerische Mathematik. vol 47. no. 2. June 1985. 217—235. doi:10.1007/BF01389710 |
[Sundance] | http://www.math.ttu.edu/~klong/Sundance/html/index.html |
[FFC] | https://bitbucket.org/fenics-project/ffc/src/master/ |
[PETSc] | https://www.mcs.anl.gov/petsc/ |