ALBERTA - An adaptive hierarchical finite element toolbox
Concepts |
ALBERTA is an Adaptive multiLevel finite element toolbox using Bisectioning refinement and Error control by Residual Techniques for scientific Applications. Its design is based on appropriate data structures holding geometrical, finite element, and algebraic information. Using such data structures, abstract adaptive methods for stationary and time dependent problems, assembly tools for discrete systems, and dimension dependent tasks like mesh modifications can be provided in a library. This allows dimension-independent development and programming of a general class of applications. In ALBERTA, hierarchical 2d and 3d meshes are stored in binary trees. Several sets of finite elements can be used on the same mesh, either using predefined ones, or by adding new sets for special applications. Depending on the currently used finite element spaces, all degrees of freedom are automatically managed during mesh modifications.
Contents
1. Introduction next Contents
The core part of every finite element program is the problem-dependent assembly and solution of the discretized problem. This holds for programs which solve the discrete problem on a fixed mesh as well as for adaptive methods which automatically adjust the underlying mesh to the actual problem and solution. In an adaptive iteration, the solution of a discrete system is necessary after each mesh change. A general finite element toolbox must provide flexibility in problems and finite element spaces while on the other hand this core part can be performed efficiently.
Data structures are needed which allow an easy and efficient implementation of the problem-dependent parts and also allow to use adaptive methods, mesh modification algorithms, and solvers for linear and nonlinear discrete problems by calling library routines.
Starting point for our considerations is the abstract concept of a finite element space defined (similar to the definition of a single finite element by Ciarlet [2]) as a triple consisting of
This leads to the definition of three main groups of data structures:
Both geometrical and finite element information strongly depend on the space dimension. Thus, mesh modification algorithms and basis functions have to be implemented for one (1d), two (2d), and three (3d) dimensions separately. Everything besides that can be formulated in such a way that the dimension only enters as a parameter (like size of local coordinate vectors, e.g.). For usual finite element applications this results in a dimension-independent programming, where all dimension-dependent parts are hidden in a library. Hence, program development can be done in 2d, where execution is usually much faster and debugging is much easier (because of simple 2d visualization, e.g., which is much more involved in 3d). With no (or maybe few) additional changes, the program will then also work in 3d. This approach leads to an enormous reduction of program development time for 3d problems.
The most time consuming part of an (adaptive) finite element code is the solution of the discrete problem. This includes the assembly of matrices and right hand sides for the (non-) linear systems which describe the discrete problem. Additionally, full information about geometry, finite elements, and algebraic structures is involved here.
On one hand, large flexibility is needed in order to choose various kinds of finite element spaces, with higher order elements or combinations of different spaces for mixed methods or systems. On the other hand, the solution of the resulting discrete systems may profit tremendously from a simple vector-oriented storage of coefficient vectors and matrices. This also allows the use of optimized solver and BLAS libraries. Additionally, multilevel preconditioners and solvers may profit from hierarchy information, leading to highly efficient solvers for the linear (sub-) problems.
The finite element toolbox ALBERTA provides the whole abstract framework like finite element spaces and adaptive strategies, together with hierarchical meshes, routines for mesh adaptation, and the complete administration of finite element spaces and the corresponding degrees of freedom (DOFs) during mesh modifications. For the solution of a specific problem, just some problem-dependent aspects like routines for assembly of system matrices and (local) error estimators have to be provided.
It is based on data structures which allow a flexible handling of information about the underlying mesh, finite element basis functions, and degrees of freedom. Such information is also used in the problem-dependent parts.
Applications of ALBERTA to several problems like CFD and phase transition problems can be found in [6], [7]. The detailed description of algorithms and data structures implemented in ALBERTA is given in [5], [8].
2. The hierarchical mesh | next | previous | Contents |
---|
The underlying mesh is a conforming triangulation of the computational domain into simplices, i.e. intervals (1d), triangles (2d), or tetrahedra (3d). The simplicial mesh is generated by refinement of a given initial triangulation. Refined parts of the mesh can be de-refined, but elements of the initial triangulation may not be coarsened. The refinement and coarsening routines construct a sequence of nested meshes with a hierarchical structure. In ALBERTA, the recursive refinement by bisection is implemented, using the notation of Kossaczký [4].
During refinement, new degrees of freedom are created. A single degree of freedom is shared by all elements which belong to the support of the corresponding finite element basis function (compare Section 3). The mesh refinement routines must create a new DOF only once and give access to this DOF from all elements sharing it. Similarly, DOFs are handled during coarsening. The actual number and locations of DOFs which have to be created or deleted depends on the finite element spaces defined on the mesh and is obtained from the DOF administration tool.
The bisectioning refinement of elements leads naturally to nested meshes with the hierarchical structure of binary trees, one for every element of the initial triangulation. Every interior node of that tree has two pointers to the two children; the leaf elements are part of the actual triangulation, which is used to define the finite element space. The whole triangulation is a list of given macro elements together with the associated binary trees. The hierarchical structure allows the generation of most information by the hierarchy, which reduces the amount of data to be stored. Some information is stored on the (leaf) elements explicitly, other information is located at the macro elements and is transfered to the leaf elements while traversing through the binary tree. Element information about vertex coordinates, domain boundaries, and element adjacency can be computed easily and very fast from the hierarchy, when needed. Data stored explicitly at tree elements can be reduced to pointers to the two possible children and information about local DOFs (for leaf elements). Furthermore, the hierarchical mesh structure directly leads to multilevel information which can be used by multilevel preconditioners and solvers.
Access to mesh elements is available solely via routines which traverse the hierarchical trees; no direct access is possible. The traversal routines can give access to all tree elements, only to leaf elements, or to all elements which belong to a single hierarchy level (for a multilevel application, e.g.). In order to perform operations on visited elements, the traversal routines call a subroutine which is given to them as a parameter. Only such element information which is needed by the current operation is generated during the tree traversal.
3. Finite elements next previous Contents
The values of a finite element function or the values of its derivatives are uniquely defined by the values of its DOFs and the values of the basis functions or the derivatives of the basis functions connected with these DOFs. We follow the concept of finite elements which are given on a single element T in local coordinates: Finite element functions on an element T are defined by a finite dimensional function space P on a reference element T and the (one to one) mapping l^{T} from the reference element T to the element T. In this situation the non vanishing basis functions on an arbitrary element are given by the set of basis functions of P in local coordinates l^{T}. Also, derivatives are given by the derivatives of basis functions on P and derivatives of l^{T}.
Each local basis function on T is uniquely connected to a global degree of freedom, which can be accessed from T. The actual memory location where a single DOF is stored on an element depends on the current set of basis functions, as well as all other finite element spaces defined on the mesh, and is accessible via the DOF administration tool.
ALBERTA supports basis functions connected with DOFs, which are located at vertices of elements, at edges, at faces (in 3d), or in the interior of elements. DOFs at a vertex are shared by all elements which meet at this vertex, DOFs at an edge or face are shared by all elements which contain this edge or face, and DOFs inside an element are not shared with any other element. The support of the basis function connected with a DOF is the patch of all elements sharing this DOF.
For a very general approach, we only need a vector of the basis functions (and its derivatives) on T and a function for the communication with the DOF administration tool in order to access the degrees of freedom connected to local basis functions. By such information every finite element function (and its derivatives) is uniquely described on every element of the mesh.
During mesh modifications, finite element functions must be transformed to the new finite element space. For example, a discrete solution on the old mesh yields a good initial guess for an iterative solver and a smaller number of iterations for a solution of the discrete problem on the new mesh. Usually, these transformations can be realized by a sequence of local operations. Local interpolations and restrictions during refinement and coarsening of elements depend on the function space P and the refinement of T only. Thus, the subroutine for interpolation during an atomic mesh refinement is the efficient implementation of the representation of coarse grid functions by fine grid functions on T and its refinement. A restriction during coarsening is implemented using similar information.
Lagrange finite element spaces up to order four are currently implemented in one, two, and three dimensions. This includes the communication with the DOF administration as well as the interpolation and restriction routines.
4. Degrees of freedom next previous Contents
Degrees of freedom (DOFs) connect finite element data with geometric information of a triangulation. For general applications, it is necessary to handle several different sets of degrees of freedom on the same triangulation. For example, in mixed finite element methods for the Navier-Stokes problem, different polynomial degrees are used for discrete velocity and pressure functions.
During adaptive refinement and coarsening of a triangulation, not only elements of the mesh are created and deleted, but also degrees of freedom together with them. The geometry is handled dynamically in a hierarchical binary tree structure, using pointers from parent elements to their children. For data corresponding to DOFs, which are usually involved with matrix-vector operations, simpler storage and access methods are more efficient. For that reason every DOF is realized just as an integer index, which can easily be used to access data from a vector or to build matrices that operate on vectors of DOF data. This results in a very efficient access during matrix/vector operations and in the possibility to use libraries for the solution of linear systems with a sparse system matrix ([3], e.g.).
Using this realization of DOFs two major problems arise:
These problems are solved by a DOF administration tool. During refinement, it enlarges the ranges of indices, if no unused indices produced by a previous coarsening are available. During coarsening, a book-keeping about used and unused indices is done. In order to reestablish a contiguous range of used indices, a compression of DOFs can be performed; all DOFs are renumbered such that all unused indices are shifted to the end of the index range, thus removing holes of unused indices. Additionally, all vectors and matrices connected to these DOFs are adjusted correspondingly. After this process, vectors do not contain holes anymore and standard operations (BLAS, e.g.) can be applied.
In many cases, information stored in DOF vectors has to be adjusted to the new distribution of DOFs during mesh refinement and coarsening. Each DOF vector can provide pointers to subroutines that implements these operations on data (which usually strongly depend on the corresponding finite element basis, compare Section 3). Providing such a pointer, a DOF vector will be transformed during mesh modifications.
All tasks of the DOF administration are performed automatically during refinement and coarsening for every kind and combination of finite elements defined on the mesh.
5. Adaptive methods next previous Contents
The aim of adaptive methods is the generation of a mesh which is adapted to the problem such that a given criterion, like a tolerance for the estimated error between exact and discrete solution, is fulfilled by the finite element solution on this mesh. An optimal mesh should be as coarse as possible while meeting the criterion, in order to save computing time and memory requirements. For time dependent problems, such an adaptive method may include mesh changes in each time step and control of time step sizes. The philosophy implemented in ALBERTA is to change meshes successively by local refinement or coarsening, based on error estimators or error indicators, which are computed a posteriori from the discrete solution and given data on the current mesh.
Several adaptive strategies are proposed in the literature, that give criteria which mesh elements should be marked for refinement. All strategies are based on the idea of an equidistribution of the local error to all mesh elements. Babuska and Rheinboldt [1] motivate that for stationary problems a mesh is almost optimal when the local errors are approximately equal for all elements. So, elements where the error estimate is large will be marked for refinement, while elements with a small estimated error are left unchanged or are marked for coarsening. In time dependent problems, the mesh is adapted to the solution in every time step using a~posteriori information like in the stationary case. As a first mesh for the new time step we use the adaptive mesh from the previous time step. Usually, only few iterations of the adaptive procedure are then needed for adaptation of the mesh for the new time step. This may be accompanied by an adaptive control of time step sizes.
Given pointers to the problem-dependent routines for assembling and solution of the discrete problems, as well as an error estimator/indicator, the adaptive method for finding a solution on a quasi-optimal mesh can be performed as a black-box algorithm. The problem-dependent routines are used for calculation of discrete solutions on the current mesh and (local) error estimates. Here, the problem-dependent routines heavily make use of library tools for assembling system matrices and right hand side for an arbitrary finite element space, as well as tools for the solution of linear or nonlinear discrete problems. On the other hand, any specialized algorithm may be added if needed. The marking of mesh elements is based on general refinement and coarsening strategies relying on the local error estimates. During the following mesh modification step, DOF vectors are transformed automatically to the new finite element spaces as described in Sections 3 and 4.
Using such a black-box algorithm, an abstract definition of basis functions, quadrature formulas and the DOF administration tool, only few parts of a finite element code depend on the dimension. Usually, these dimension-dependent parts are hidden in a library. This allows a dimension-independent programming of applications to the greatest possible extent.
6. References previous Contents