The data structure

Introduction

oomph-lib is big! This document gives a "bottom up" overview of the library's data structure and discusses how the various objects interact. In addition to the detailed discussion provided below, the following doxygen-generated lists/indices provide quick access to the documentation of oomph-lib's classes: clicking on the hyperlinks associated with a class takes you directly to a detailed description of its inheritance structure and its members.

This rest of this document provides a "bottom up" overview of the data structure in oomph-lib and discusses how the various objects interact. For brevity, we usually replace the list of arguments to functions by `(...)' and explain the main input and output parameters in words. The full specifications of the interfaces may be found in the individual class documentation, accessible via the links at the top of this page.

Technical terms

Throughout this document, certain commonly used terms have a specific technical meaning:

Overview of the basic data structure

The main components of oomph-lib are Data, Node, GeneralisedElement, Mesh and Problem.

Data

The most elementary data structure in oomph-lib is Data (ha!).

Consider the solution of a scalar PDE (e.g. a Poisson equation) with certain boundary conditions. The numerical solution of this problem requires the computation of the function values (double precision numbers) at a finite number of spatial positions (the Nodes). Typically, these values fall into two categories: those that are known a priori (i.e. are enforced by boundary conditions), and those that must be determined as part of the solution.

Data stores a value — a double precision number. Typically, the values of the unknowns are determined by the solution of a system of algebraic equations. The solution of this system usually requires a (linear) numbering of the unknowns and associated equations. Hence, Data also stores a (long) integer that represents the number of the unknown in the global numbering scheme. Convention: If the Data value is pinned, we set the equation number to the static member data Data::Is_pinned, a negative number.

The number of an unknown is related to the number of the equation that `determines its value', so we use the terms `equation number' and `number of the unknown' interchangeably. In fact, because the term `number of the unknown' is rather tedious, we only use the term `equation number'.

Two observations motivate a straightforward extension of this basic data structure:

Therefore, Data allows the storage of multiple values (all of which can either be pinned or free, and all of which have their own global equation number); Data can also store a certain number of auxiliary (history) values used for timestepping. Finally, Data stores a pointer to a TimeStepper whose member functions relate the history values to the values' time-derivatives.

Direct, pointer-based read/write access to the Data values is provided by the functions

Data::value_pt(i)

which returns a pointer to the i-th value at the present time, and by

Data::value_pt(t,i)

which returns a pointer to the t -th history value associated with value i. Read-only access is also provided by the functions

Data::value(i)

and its time-dependent counterpart

Data::value(t,i)

We recommend using these functions instead of the pointer-based access functions for read access because the two Data::value(...) functions are overloaded in the Node class (discussed below) so that they return suitably constrained nodal values if a Node is hanging. The Data::value(...) functions cannot be used to set Data values. For this purpose we provide the functions

Data::set_value(i,val)

which sets the i -th Data value to the double precision number val; and its time-dependent counterpart

Data::set_value(t,i,val)

which sets the t -th history value associated with the i -th Data value to the double precision number val.

The general convention for all time-dependent data is that the index t=0 refers to values at the present time, whereas the values associated with t>0 correspond to history values. In many cases (e.g. BDF schemes) these history values are simply the values at previous timesteps, but this is not guaranteed. See the section Time-stepping for further details.

Nodes

In FE computations, most (but not all; see below) Data are associated with nodal points. Conversely, all Nodes in a finite element mesh have Data associated with them. Nodes are therefore derived from Data, but also store a spatial position, specified by a certain number of spatial (Eulerian) coordinates.

The nodal positions are accessed by the member function

Node::x(i)

which returns the current value of i-th nodal coordinate or

Node::x(t,i)

which returns the value of i-th nodal coordinate at the present (t=0) timestep or a history value, if (t>0); again, note that the history values are not necessarily positions at previous timesteps.

Advanced features:

Nodes have the following additional features:

Elements

Overview

Most (finite-)elements in oomph-lib have a four-level inheritance structure which separates:

  1. the basic functionality that is shared by all (generalised) elements.
  2. the functionality that is shared by all finite elements.
  3. the implementation of the finite element geometry (i.e. the shape-function-based mapping between local and global coordinates).
  4. the representation of the mathematics that describes a specific problem.

The distinction between geometry and `maths' greatly facilitates code-reuse as a (geometric) quad-element, say, can form the basis for elements that solve a variety of equations (e.g. Poisson, Advection–Diffusion, Navier–Stokes, ...). We shall now discuss the four levels of the element hierarchy in more detail.

Level 0: GeneralisedElement

The class GeneralisedElement forms the base class for all elements in oomph-lib. It incorporates the basic functionality that all elements must have. The interfaces at this level are so general that a GeneralisedElement can represent discrete algebraic constraints (or even finite-difference stencils).

The main role of elements is to provide contributions to a global residual vector and to the global Jacobian matrix. The two virtual functions

virtual void GeneralisedElement::get_residuals(...)
virtual void GeneralisedElement::get_jacobian(...)

specify the appropriate interfaces.

In multi-physics problems, the elemental residuals and Jacobian will be a combination of the residuals vectors and Jacobian matrices of the constituent single-physics elements. An obvious implementation is to use multiple inheritance and function overloading

class MultiPhysicsElement : public virtual SinglePhysicsOneElement,
public virtual SinglePhysicsTwoElement
{
[...]
void get_residuals(...)
{
SinglePhysicsOneElement::get_residuals(...);
SinglePhysicsTwoElement::get_residuals(...);
}
[...]
};

where the MultiPhysicsElement inherits from SinglePhysicsOneElement and SinglePhysicsTwoElement.

A problem with this implementation arises when we consider where to initialise the residuals vector. If the second single-physics get_residuals(...) function initialises the residuals vector, then the contribution of the first single-physics element will be negated. When writing a single-physics element, however, we cannot know whether it will ever be used as part of a multi-physics element and, if so, in which order the get_residuals(...) functions will be called. The solution adopted in oomph-lib is to provide the two additional virtual functions

virtual void GeneralisedElement::fill_in_contribution_to_residuals(...)
virtual void GeneralisedElement::fill_in_contribution_to_jacobian(...)

which must not initialise the residuals or Jacobian, but merely add the contribution of the element to the vector or matrix.

We then use the default implementation

virtual void GeneralisedElement::get_residuals(Vector<double> &residuals)
{
//Zero the residuals vector
residuals.initialise(0.0);
//Add the elemental contribution to the residuals vector
fill_in_contribution_to_residuals(residuals);
}

which permits a simple multi-physics re-implementation

class MultiPhysicsElement : public virtual SinglePhysicsOneElement,
public virtual SinglePhysicsTwoElement
{
[...]
void get_residuals(Vector<double> &residuals)
{
//Zero the residuals vector
residuals.initialise(0.0);
//Add the first elemental contribution to the residuals vector
SinglePhysicsOneElement::fill_in_contribution_to_residuals(residuals);
//Add the second elemental contribution to the residuals vector
SinglePhysicsTwoElement::fill_in_contribution_to_residuals(residuals);
}
[...]
};

The default implementation of fill_in_contribution_to_jacobian(...) uses finite differences to calculate the Jacobian matrix. Hence, the simplest possible implementation of a new element requires only the specification of fill_in_contribution_to_residuals(...).

Level 1: FiniteElement

The class FiniteElement is derived from GeneralisedElement and incorporates the basic functionality that all finite elements must have.

Level 2: Geometric Elements

At this level, we specify the element geometry and the mapping between local and global coordinates. Wherever possible, templating has been (and, for any newly developed elements, should be) used to formulate the elements in a dimension-independent way. For instance, Lagrange-type 1D line, 2D quad and 3D brick elements are implemented in the doubly-templated QElement<DIM,NNODE_1D> class. The template parameters indicate the spatial dimension of the element and the number of nodes along the element's one-dimensional edges. Hence, QElement<1,3> is a three-node line element with a quadratic mapping between local and global coordinates; QElement<3,2> is an 8 node brick element with a trilinear mapping between local and global coordinates. The dimension and the number of Nodes must be set by calling the appropriate set_ functions in the constructor, see above.

The most important member functions implemented at this level include

Finally, we specify a pointer to a spatial integration scheme (usually a Gauss rule). The order of the integration scheme is based on the order of the interpolation in the isoparametric mapping. If this is inappropriate for an element that is derived from a given geometric element, the default assignment can be over-written at a higher level. This is discussed in more detail in a separate document.

Level 3: `The Maths'

At this level, we implement the equations that are represented by the specific element. We implement the interpolation(s) for the unknown function(s), employing either the geometric shape functions that already exist on level 2, or employing additional shape functions defined at this level. This allows us to write further member functions such as interpolated_u(...), say, which compute the i-th velocity component at the local coordinate s like this

// Create fluid element
SomeFluidElement fluid_element;
[...]
// Vector of local coordinates
Vector<double> s(3);
// Vector of velocity components
Vector<double> u(3);
// Compute the velocity at local coordinate s
fluid_element.interpolated_u(s,u);

We introduce wrapper functions to access function values, so that we can formulate "The Maths" in generic terms. Rather than referring to the pressure at node n, via

unsigned n;
double press = Node_pt[n]->value(3);

(which forces us to remember that in this particular 3D fluid element, the pressure is stored as the fourth value at all nodes...), say, we provide an access function p_fluid(...) which allows us to write

unsigned n;
double press = fluid_element.p_fluid(n);

When writing these wrapper functions, direct access to the nodal values should be avoided to ensure that the element remains functional in the presence of hanging nodes. Hence, the wrapper functions should make use of the Node::value(...) functions as in this example

double SomeFluidElement::p_fluid(const unsigned& n)
{ return Node_pt[n]->value(3);}

[See section Hanging Nodes for a full description of hanging nodes.]

The functions

FiniteElement::assign_additional_local_eqn_numbers()
FiniteElement::get_residuals(...)
FiniteElement::get_jacobian(...)

that were defined (as virtual functions) in FiniteElement can now be implemented for the specific system of equations that are represented by this element. The function assign_additional_local_eqn_numbers() is called by FiniteElement::assign_local_eqn_numbers() and may be used to assign local equation numbers that correspond to particular physical variables. For example, in QCrouzeixRaviart "fluid" elements, the pressure is stored as 'internal Data', so an internal array P_local_eqn could be defined by

P_local_eqn[i] = internal_local_eqn(0,i)

but in QTaylorHood "fluid" elements, the pressure is stored as nodal Data, so that

P_local_eqn[i] = nodal_local_eqn(Pconv[i],DIM+1)

In the above, Pconv[i] is an array that returns the local node number at which the i -th pressure freedom is stored and DIM is the dimension of the element. The use of such an array introduces a memory overhead, however, because each element must permanently store these additional integers. In general, we prefer to use member functions for this purpose. In QCrouzeixRaviart elements, for example,

int p_local_eqn(const unsigned &n)
{
return internal_local_eqn(0,n);
}

Thus, in our standard equations the function assign_additional_local_eqn_numbers() is not used.

Finally, the virtual function

FiniteElement::required_nvalue(...)

should be implemented to specify the number of values that are stored at each of the element's local Nodes. The default number of values stored at a Node is zero.

`Advanced' features for the `Maths' level:

It often makes sense to subdivide the `Maths' level further into

  1. A class that contains the abstract FE formulation of the mathematical problem.
  2. A class that combines the mathematical formulation with a specific geometrical element.

An example is given by the QPoissonElements which inherit their maths from PoissonEquations (templated by the spatial dimension) and their geometry from QElement (templated by the spatial dimension and the number of nodes). PoissonEquations specifies the weak form of the Poisson equation in terms of (virtual) shape and test functions. The QPoissonElements turn this abstract formulation into a specific (isoparametric) element, by specifying both test and shape functions as the geometric shape functions defined in QElement.

Important convention regarding element constructors

To facilitate mesh generation and adaptation it is important that element constructors should not have any arguments! If arguments must be passed to an element (e.g. function pointers to source functions, etc.), this should be done after the mesh generation, usually in the Problem constructor. In adaptive mesh refinement procedures any function/data pointers in newly created elements are set to be the same as those of the father element. If this is not the desired behaviour, arguments should be passed to the elements in the function Problem::actions_after_adapt().

Meshes

At its most basic level, a Mesh is simply a collection of elements and Nodes, accessed by the vectors of pointers, Mesh::Element_pt and Mesh::Node_pt, respectively. To facilitate the application of boundary conditions, we also store vectors of pointers to the (Boundary)Nodes that lie on the boundaries of the mesh. They are accessible via the member function

Mesh::boundary_node_pt(i,j)

which returns a pointer to the j -th node on the i -th boundary of the mesh.

Mesh generation

The oomph-lib data structure is designed to make the mesh generation process generic so that meshes developed for one particular problem can easily be re-used in others. For this reason, it is generally assumed that a Mesh contains only elements of a single type and that this element type is passed to the Mesh constructor as a template parameter. Problems that require multiple element types should use multiple meshes, which are discussed in more detail below. It is possible to mix types of elements within a single Mesh, if desired; this can be advantageous when the element types are very closely related, for example "free surface" elements in a mesh of "Fluid" elements.

Mesh generation (usually performed in the Mesh 's constructor) then works as follows:

  1. Create the first element and add the pointer to it to the Mesh::Element_pt vector. (We know which element to build because we have passed its type as a template parameter.) Since element constructors do not take any arguments, this step is completely generic and ensures that a Mesh that was originally created to solve a Poisson equation, say, can also be used to solve Navier-Stokes equations (provided the element topology is the same, i.e. provided the elements are derived from the same type of geometric element (e.g. quad or triangle). The element now exists but does not know anything about its Nodes etc.
  2. Loop over the element's Nodes and for each Node:
    1. Create the Node using the element's
      FiniteElement::construct_node(...)
      member function. As discussed above, this function creates Nodes of exactly the right type and fills in the element's own Node_pt vector.
    2. FiniteElement::construct_node(...) returns a pointer to the newly created Node; add this to the Mesh::Node_pt vector.
    3. Assign the nodal coordinates.
    The Node is now fully functional and (by default) all values that are stored with it are free (i.e. not pinned).
  3. If the Node is located on a mesh boundary then it must be a BoundaryNode. BoundaryNodes can be created using the function
    FiniteElement::construct_boundary_node(...)
    in place of FiniteElement::construct_node(...). Alternatively, if the Node has already been created, it can be upgraded to a BoundaryNode by using the function
    Mesh::convert_to_boundary_node(...)
    The function
    Mesh::add_boundary_node(i,&node)
    should then be used to add (a pointer to) the BoundaryNode to the Mesh's boundary-storage scheme. In addition, the function Mesh::add_boundary_node() passes boundary information to the Node itself.
  4. Create the next element and loop over its Nodes. Some Nodes will already exist (because they have been created by the first element). For such Nodes, we merely add the pointer to the existing Nodes to the element's Node_pt vector. If a Node does not exist yet, we create it, as discussed above.
  5. Keep going until all the elements and Nodes have been built.

Equation/DOF numbering at the mesh level

Now that the Mesh is assembled, we can set up the numbering scheme for the unknowns that are associated with the Nodes and with the elements' `internal Data '. (For problems that involve `external Data ', i.e. Data that is not associated with Nodes and elements, a further step is required; see section Problems below). As discussed above, whenever a Data object is created (either as part of a Node in the mesh or as `internal Data ' inside an element), its values are assumed to be free (i.e. not pinned). Before we can set up the equation numbering scheme, we must pin all those Data values that are prescribed by boundary conditions. This is generally done at the Problem level (see below) and for the subsequent discussion, we assume that this step has already taken place.

The equation numbering scheme must achieve two things:

These tasks are performed in two steps. The function

Mesh::assign_global_eqn_numbers(...)

whose argument is a vector of pointers to doubles, dof_pt, assigns the (global) equation numbers for the element's internal Data and for all Data associated with the Mesh's Nodes. On return from this function, dof_pt[i] points to the value of the i-th global degree of freedom. The setup scheme (which is fully implemented in oomph-lib) works as follows:

  1. Loop over all Nodes in the Mesh.
    • Loop over all values that are stored at that Node.
      • If the value is free (i.e. not pinned), add the pointer to the value to the dof_pt vector.
  2. Loop over all elements.
    • In every element, loop over the `internal Data '
      • For every instance of `internal Data ', loop over its values.
        • If the value is free (i.e. not pinned), add pointer to the value to the dof_pt vector.

Once this has been done (for all Meshes, if there are multiple ones), the function

Mesh::assign_local_eqn_numbers(...)

loops over all elements and executes their

FiniteElement::assign_local_eqn_numbers();

member function to set up the elements' lookup table that translates between local and global equation numbers.

Problems

Finally, we reach the highest level of the oomph-lib hierarchy, the Problem itself. The generic components of Problem, provided in the base class of that name, are:

Multiple Meshes

The mesh generation process previously described was for meshes that contain only elements of a single type. The process could easily be generalised to meshes that contain multiple element types by providing multiple template parameters. When solving a fluid-structure interaction (FSI) problem we could, therefore, create a single mesh that discretises both the fluid and the solid domains with appropriate elements. However, to facilitate code re-use it is desirable to keep meshes as simple as possible so that a mesh that was originally developed for a pure fluids problem can also be used in an FSI context. For this reason, the Problem class allows a problem to have multiple (sub-)meshes. Pointers to each sub-mesh must be stored in Problem's Sub_mesh_pt vector by using the function

Problem::add_sub_mesh(Mesh* const &mesh_pt)

However, many of the generic operations within oomph-lib (equation numbering, solving, output of solutions,...) involve looping over all elements and Nodes in the problem. Therefore, if a Problem contains multiple (sub-)meshes, the sub-meshes must be combined into a single global Mesh whose Element_pt and Node_pt vectors provide ordered access to the elements and Nodes in all submeshes. The function

Problem::build_global_mesh()

combines the sub-meshes into the global Mesh and must be called once all the sub-meshes have been constructed and "added" to the Problem.

Important: Many operations (such as the shifting of history values in time-stepping) must be performed exactly once for each Node (or Data). Therefore, the vector of (pointers to) nodes in the global Mesh must not contain any duplicate entries. When copying (pointer to) Nodes from the submeshes into the global Mesh, the function Problem::build_global_mesh() ignores any Nodes that have already been copied from a previous submesh. [The Mesh::self_test() function checks for duplicates in the Mesh::Node_pt vector.]

Convention: Recall that (Boundary)Nodes are `told' about the number of the boundary they live on when the (sub-)meshes are constructed. In the context of multiple meshes, this raises the question if this number should continue to refer to the boundary number within the submesh or be updated to a boundary number within the global Mesh. We adopt the convention that boundary numbers remain those that were originally assigned when the submeshes were constructed.

Multiple meshes and adaptivity: Mesh adaptation is performed separately for each submesh. Following the adaptation of one or more submeshes (in the process of which various Nodes/elements will have been created/deleted), we must also update the global Mesh. This is done by calling

Problem::rebuild_global_mesh();

and this function is executed automatically if the mesh adaptation/refinement is performed by

Problem::refine_uniformly();

or

Problem::adapt();

Problem Construction

Here's an overview of how Problems are set up and solved in oomph-lib. For simplicity, we illustrate the process for a problem with a single Mesh that contains elements of a single type.

  1. Define the element type and the TimeStepper (if the problem is time-dependent).
  2. Build the Mesh, passing the element type as a template parameter and the TimeStepper as an argument to the Mesh constructor. (Typically Mesh constructors take a pointer to a TimeStepper as an argument since the TimeStepper needs to be passed to the element's FiniteElement::construct_node(...) function. If the Problem has no time-dependence, we can pass a pointer to the static Mesh::Steady timestepper; many Mesh constructors use a pointer to this TimeStepper as the default argument).
  3. Create the Problem 's global Data (if it has any).
  4. Apply the essential boundary conditions by pinning the appropriate values; all other values remain free.
  5. We now have a fully assembled Mesh and all elements know their constituent Nodes. However, because element constructors are not allowed to have any arguments, all but the simplest of elements will now have to be provided with additional information. For instance, we might need to set pointers to source functions, etc.
  6. Assign the equation numbers to all unknowns in the problem. This is accomplished in a two-stage process by
    Problem::assign_eqn_numbers(...)
    which
    1. Assigns the equation numbers for all global Data (if any)
    2. Loops over the submeshes to perform the global equation numbering for all values associated with the Meshes (i.e. Nodes and elements), using
      Mesh::assign_global_eqn_numbers(...)
    3. Loops over the submeshes again to perform the local equation numbering for all elements associated with the meshes using
      Mesh::assign_local_eqn_numbers(...)
    [Note that we have to assign the global equation numbers for all meshes before we assign any local equation numbers. This is because the external Data of elements in one mesh might be nodal Data in another mesh – think of fluid-structure interaction problems]. At the end of this step, we will have filled in the Problem's Dof_pt vector which holds the pointers to the unknowns. [Note: The equation numbering scheme that is generated by the above procedure is unlikely to be optimal; it can (and in the interest of efficiency probably should) be changed afterwards by a problem-specific renumbering function.]
  7. Before we can solve the Problem we will usually have to perform a few additional (problem-dependent) initialisation steps. For instance, we might want to assign initial and/or boundary conditions and provide initial guesses for the unknowns. For time-dependent problems, the Problem class provides a member function
    Problem::assign_initial_values_impulsive()
    which creates a past history for all time-dependent unknowns, assuming an impulsive start. If you want a `smooth' start from a given previous time-history, you will have to implement this yourself; in this case, consult the section Time-stepping which outlines the time-stepping procedures in oomph-lib.
  8. Now we can solve the Problem by calling
    Problem::newton_solve()
    which employs the Newton-Raphson method to solve the nonlinear system of equations that is specified (formally) by the global Jacobian matrix and the global residual vectors. The function Problem::newton_solve() employs a linear solver which may be specified by Problem::linear_solver_pt(). The default linear solver is SuperLU . The Newton iteration proceeds until the maximum residual falls below Problem::Newton_solver_tolerance. [If the number of iterations exceeds Problem::Max_newton_iterations or if the the maximum residual exceeds Problem::Max_residuals, the Newton solver throws an error]. When a solution has been found, all unknowns in the problem (which are accessible to the Problem via its Dof_pt vector) are up-to-date. The Problem class also provides the function
    Problem::unsteady_newton_solve(...)
    for time-dependent problems. Given the current values of the unknowns (at time t = Problem::Time_pt->time() ) and their past histories, this function determines the solution at the advanced time t+dt. See section Time-stepping for more details on the conventions used in timestepping.

The Problem class also has member functions which assemble the global residual vector and the global Jacobian matrix.

Important:

Since the unknowns in the Problem are accessed directly via pointers, their values are automatically updated during the Newton iteration. If the Problem has any auxiliary parameters that depend on the unknowns, their values need to be updated whenever an unknown might have changed (i.e. after every step of the Newton iteration). For such cases, the Problem class provides the four (empty) virtual functions

Problem::actions_before_newton_step()
Problem::actions_after_newton_step()
Problem::actions_before_newton_solve()

and

Problem::actions_after_newton_solve()

which are executed before/after every step of the Newton iteration and before/after the nonlinear solve itself, respectively. In addition, the virtual function

Problem::actions_before_newton_convergence_check()

is executed before the residuals are calculated in the Newton solver.

When you formulate your own Problem, you will have to decide what (if anything) should live in these functions. Typical examples of actions that should be taken before a solve are the update of any boundary conditions. If the boundary conditions depend upon variables in the problem, they must be updated before every Newton step and should therefore be placed in Problem::actions_before_newton_step(). Actions that take place after a Newton step or solve would include things like updating the nodal positions, writing output or any other post-processing/solution monitoring. For example, if the solution after each Newton step were to be documented this could be accomplished by calling a suitable output function in the Problem::actions_after_newton_step() function.

In many cases, only the Problem::actions_before_newton_convergence_check() function is required. On entry to the Newton solver, the initial residuals are computed and checked, so the function is executed before any Newton steps are taken. After each Newton step, the residuals vector is recomputed and checked, so the function is also called after every Newton step (or before the next Newton step). Nonetheless, we provide all five functions for the greatest possible flexibility.

Time-stepping

Time

Time-derivatives are generally evaluated by finite difference expressions, e.g. by BDF schemes. Hence, within the code, functions are only ever evaluated at discrete time levels. The class Time contains (a pointer to) the `current' value of the continuous time and a vector of current and previous timesteps so that the value of the continuous time at any previous timestep can easily be reconstructed. This is useful/necessary if there are any explicitly time-dependent parameters in the problem. The general convention within all timestepping procedures is to associate the values at the `present' time (i.e. the time for which a solution is sought) with time level `0', those at the previous time level (where the solution is already known) with time levels `1', `2' etc. The function

Time::time(t)

therefore returns the `current' value of the continuous time if t=0, the continuous time before the previous timestep if t=1, etc.

Basic Time-stepping

The base class TimeStepper provides the basic functionality required to evaluate time derivatives, and to keep track of the time histories of the unknowns. Primarily, a TimeStepper stores the coefficients (weights) that allow the evaluation of time-derivatives (up to a certain order) in terms of the history values stored in Data. Synchronisation of multiple TimeStepper s is ensured by providing them with pointers to the Problem's (single) Time object.

Here's an illustration of the time-stepping procedure for an implicit scheme.

It is important to understand how the shifting of the timesteps (in preparation for the next timestep) is performed because certain default operations will not be appropriate for all Elements/ Problems. Recall that in time-dependent problems, Data::value_pt(i) points to the current (and, as yet, unknown) Data values, while Data::value_pt(t,i) for t>0 points to the history values that the TimeStepper uses to work out time-derivatives. When we move to the next timestep, the history values need to be adjusted in the manner that is appropriate for the timestepping scheme. For instance, for BDF schemes, all values need to be pushed back by one time level. This operation is performed by the function

TimeStepper::shift_time_values(...)

In the case of Nodes, the shifting of values associated with the nodal positions is performed by the function

TimeStepper::shift_time_positions(...)

To ensure that all Data in the Problem is shifted once (and only once!) Problem::shift_time_values() performs the following operations:

  1. Shift the values of the time history stored in the Time object.
  2. Shift the time values in the global Mesh. This involves the following steps:
    1. Loop over all elements and call TimeStepper::shift_time_values(...) for the TimeStepper corresponding to each internal Data value. This leads to the slightly ugly construction
      internal_pt(i)->time_stepper_pt()->shift_time_values(internal_pt(i))
      but there appears to be no way to avoid this.
    2. Loop over all Nodes in the mesh and
      1. call TimeStepper::shift_time_values(...) for each Node 's TimeStepper.
      2. call TimeStepper::shift_time_positions(...) for each Node's positional TimeStepper.
  3. Shift the time values for all global Data.

Adaptive time-stepping

In adaptive time-stepping, the size of the timestep is adjusted automatically, so that the global (temporal) error estimate, computed by the Problem member function

Problem::global_error_norm()

remains below a preset threshold. However, the function Problem::global_error_norm() must be implemented for each specific problem. The error norm is usually constructed from the (estimated) errors of individual Data values. Estimates for these quantities are given by the differences between the actual value and a predicted value, as determined by

TimeStepper::error_in_value(...)

and the errors in positional Data values, found by

TimeStepper::error_in_position(...)

In moving mesh problems, a suitable norm is the root-mean-square of the errors in all positional coordinates at every Node. In fluid problems, the error is usually based on the velocity components, etc.

Once a suitable norm has been chosen, a single adaptive timestep is taken by the function

Problem::adaptive_unsteady_newton_solve(...)

This function returns a double precision number that is the value of dt that should be taken at the next timestep. A typical calling loop is thus

// Initial suggestion for timestep
double dt=0.001;
for(unsigned t=0;t<100;t++)
{
// Try to take a timestep dt forward; if the computed
// solution is not sufficiently accurate, reduce the
// the timestep and repeat... Return the suggestion for
// next timestep
double dt_next = adaptive_unsteady_newton_solve(dt,...);
dt = dt_next;
}

Within the Problem::adaptive_unsteady_newton_solve(..) function, if the global error norm is too large, the step is rejected, the timestep reduced and the step repeated. If the timestep falls below a preset tolerance Problem::Minimum_dt (which has the default value 1.0e-12), the program will terminate. It is also possible to set a maximum timestep by over-writing the (large) default for Problem::Maximum_dt (initialised to 1.0e12) with a smaller value.

Restarts

Time-dependent simulations can consume a lot of computer time, and it is essential to be able to restart simulations rather than having to re-do them from scratch if a computer crashes during the code execution. For this purpose the Problem class provides the two member functions

Problem::dump(...)

and

Problem::read(...)

which write/read the generic Problem components (i.e. the Data values and history values, the history of previous timestep increments, etc) to/from disk. These generic functions are typically called from within a specific Problem's own dump/read functions which also deal with any additional, problem-specific data that has to be recorded/reset to allow a proper restart.

Problem/Mesh adaptation

Adaptation at the Problem level

The ability to adaptively refine/unrefine the Problem's mesh(es) in regions in which the solution undergoes rapid/slow variations is of crucial importance for the overall efficiency of the solution process. Mesh-adaptation involves the following steps:

  1. Compute an estimate of the error in the computed solution for all elements in the mesh.
  2. Label the elements whose error estimate is the above the maximum (or below the minimum) permissible error.
  3. Perform the actual mesh adaptation: elements whose error is too large are subdivided into `son' elements; elements whose error is too small are merged with their `siblings', provided the `siblings' are also scheduled for de-refinement. Usually, certain bounds are imposed on the maximum and minimum refinement levels to prevent excessive mesh refinement (e.g. near singularities) or de-refinement. (While the use of very large elements might be permissible in regions in which the solution varies little, such elements will provide a poor representation of the domain shape and thus lead to unsatisfactory post-processing.)
  4. Represent the previously computed solution on the newly created Nodes / elements.
  5. Once all the Problem's submeshes have been adapted, update the Problem itself by updating the global Mesh, re-generating the equation numbering scheme, etc.
  6. Solve the adapted Problem.

The Problem class provides several functions that perform these tasks completely automatically. For instance, the function

Problem::refine_uniformly()

performs one uniform mesh refinement step for all (refineable) submeshes in the Problem. Similarly, the function

Problem::adapt(...)

performs one mesh adaptation step. In both functions, mesh adaptation is followed by the update of the global Mesh and the re-assignment of the equation numbers so that, on return from these functions, the Problem can immediately be solved again.

The Problem class also provides overloaded versions of the steady and unsteady Newton solvers

Problem::newton_solve(...)

and

Problem::unsteady_newton_solve(...)

that automatically perform mesh adaptations until the computed solution satisfies the specified error bounds on all submeshes (or until a max. number of adaptations has been performed). The (empty) virtual member function

Problem::set_initial_condition()

establishes the interface to the function that sets all Data to their initial conditions. This function must be overloaded if nontrivial initial conditions are to be applied. (If mesh adaptations are performed while the first timestep is computed, the initial conditions on the adapted mesh can usually be represented more accurately by re-assigning them, rather than by interpolation from the coarse initial mesh).

Adaptation at the Mesh level

The ability to perform the adaptation at the Problem level relies on the availability of standardised interfaces to functions that handle the adaptation on the Mesh level. These interfaces are provided in the class RefineableMeshBase. RefineableMeshBase is derived from Mesh and stores a pointer to a spatial error estimator, as well as double precision numbers representing the target error levels for the adaptation. The member function

RefineableMeshBase::refine_uniformly(...)

performs one uniform mesh refinement step. Similarly,

RefineableMeshBase::adapt(...)

adapts the mesh according to the specified error bounds, using the mesh's spatial error estimator to compute the elemental errors.

The details of the mesh adaptation process depend on the mesh topology; currently the virtual functions in RefineableMeshBase are implemented in a general form for quad and brick elements. We shall discuss the mesh adaptation process in detail for meshes of a particular type: The RefineableQuadMesh class implements the mesh adaptation procedures for two-dimensional, block-structured meshes which consist of the refineable variant of 2D QElements. The description provides a template for the development of mesh refinement procedures for meshes with different element topologies (e.g. triangular elements, or 3D QElements).

Mesh adaptation for 2D quadrilateral meshes

Refineable QElements

The abstract base class RefineableQElement<2> `upgrades' existing elements of (geometric) type QElement<2,NNODE_1D> to versions that are suitable for use in adaptive computations. 'Upgrading' is achieved via inheritance so that, e.g., refineable Poisson elements are defined as:

template <unsigned DIM, unsigned NNODE_1D>
class RefineableQPoissonElement : public QPoissonElement<DIM,NNODE_1D>,
public virtual RefineablePoissonEquations<DIM>,
public virtual RefineableQElement<DIM>
{
[...]
}

The abstract base class RefineableQElement<2> defines virtual interfaces for those FiniteElement member functions that might have to be overloaded in the refineable version. In most cases, these member functions must be re-implemented to deal with the possible presence of hanging nodes, see below.

QuadTree procedures for mesh refinement

Many of the mesh adaptation procedures for meshes of type RefineableQuadMesh use quadtree representations of the mesh. The quadtree navigation and search algorithms are based on those described in H. Samet's "The design and analysis of spatial data structures" (Addison-Wesley, 1990). [Note: Unfortunately, in the usual tree terminology, quadtrees are made of "nodes" which are, of course, completely unrelated to the Nodes in the finite element mesh! The context and – within this document – the different typefaces should make it clear which is which...] It is important to understand that each RefineableQElement<2> has an associated QuadTree and each QuadTree has an associated RefineableQElement<2>. This two-way "has a" relationship permits a "clean" implementation of the (generic) QuadTree algorithms, although it does incur the cost of two additional pointers.

To illustrate the way in which RefineableQuadMeshes are represented by QuadTrees, the figure below shows a simple finite element mesh together with its quadtree-representation. There are two different types of quadtree classes: QuadTrees and QuadTreeRoots, which inherit from QuadTrees. The overall structure of the quadtree is defined by pointers between its "nodes". Each "node" (an object of type QuadTree, shown in pink) in the quadtree has a pointer to its "father" (if this pointer is NULL, the "node" is the "root node") and a vector of pointers to its four "sons" (if the vector is empty, the "node" is a "leaf \e node"). This data structure is sufficient to navigate the quadtree (e.g. identify the quadtree's "leaf nodes", determine a "node"'s neighbours, etc.) Each QuadTree also stores a pointer to an associated "object" of type RefineableElement (shown in light blue). The finite element mesh that is used in the computations only comprises those RefineableElements that are associated with "leaf nodes". We refer to these elements as "active elements". In the diagram below, the active elements are identified by thick blue boundaries and the blue element numbers correspond to those in the mesh.

QuadTree.gif
Illustration of a QuadTree and the associated RefineableQuadMesh

Mesh generation

Any Mesh that is designed to contain 2D QElement s forms a suitable basis for a RefineableQuadMesh mesh and the initial stages of the mesh generation process can be performed exactly as described in section Meshes above. Typically, the constructor for a RefineableQuadMesh will produce a relatively coarse background mesh (sometimes referred to as a "root mesh") which can subsequently be refined as required. As discussed before, the type of element will typically be passed as a template parameter and it is assumed that in any concrete instantiation of the RefineableQuadMesh class, the element is derived from the class RefineableQElement<2>.

Following the generic setup of the mesh (creating elements, Nodes, etc), the constructor of the RefineableQuadMesh must associate each RefineableQElement<2> in the mesh with a QuadTreeRoot and vice versa. The association between RefineableQElement<2>s and QuadTreeRoots is established by the QuadTreeRoot constructor which takes the pointer to the RefineableQElement<2> as its argument. The different QuadTrees must then be combined into a QuadTreeForest, whose constructor establishes each QuadTree 's N/S/W/E neighbours and their relative orientation. Here is an illustration:

QuarterCircleMesh.gif
Mesh consisting of three RefineableQElement<2>s/ QuadTrees. Element/QuadTree numbers are blue, and the QuadTree's `northern' direction is indicated by the red arrows. The QuadTrees are combined to a QuadTreeForest which establishes their adjacency (QuadTree 2 is QuadTree 0's northern neighbour etc.) and their relative orientation (QuadTree 1's North coincides with QuadTree 2's West, etc.) Note that the RefineableQElement<2>'s local coordinates s_0 and s_1 are aligned with the QuadTree 's W-->E and S-->N directions, respectively.

The virtual member function

RefineableQuadMesh::setup_quadtree_forest()

creates the QuadTreeForest automatically and should be called at the end of the RefineableQuadMesh's constructor. Finally, the mesh must be given a pointer to a spatial error estimator.

Hanging Nodes

Before explaining the details of the mesh adaptation process we discuss how hanging nodes are treated in oomph-lib. The figure below shows a RefineableQuadMesh that originally consisted of four 4-node RefineableQElements and nine nodes (nodes 0 to 8). The mesh was adapted twice: During the first adaptation, the top right element was subdivided into four son elements and five new Nodes (nodes 9 to 13) were created. Then one of the newly created elements was subdivided again and in the process Nodes 14 to 18 were created.
On this mesh, inter-element continuity is not ensured unless the hanging Nodes (=internal Nodes not shared by an adjacent element – here Nodes 9, 10, 14, 15, 17 and 18) are suitably constrained. For instance, the nodal values at Node 10 must be linear combinations of the values at nodes 4 and 7 if the solution is to remain continuous across the eastern edge of element 2. We refer to nodes 4 and 7 as the "master nodes" of Node


  1. The hanging node status of the Nodes can be determined by the following procedure:

If the difference in the refinement levels of two adjacent elements is greater than one, some of the master nodes determined by the above procedure will themselves be hanging nodes and these now need to be eliminated recursively. For instance, Nodes 9 and 11 are initially classified as master nodes for Node 17. However, Node 9 is a hanging node which depends on Nodes 4 and 5. Hence, the ultimate master nodes for Node 17 are Nodes 4, 5 and 11.

HangingNodes.gif
Illustration of hanging nodes
The presence of hanging nodes must not affect the node-based representations for the unknown function(s) and for the geometric mapping within the elements. Hence, within any element $ E $, a scalar function $ u $ should continue to be representable as

\[ u = \sum_{j} U_{J(j,E)} \psi_j(s) \mbox{\hspace{5cm} (1)} \]

while the mapping between local and global coordinates should retain the form

\[ {\bf x} = \sum_{j} {\bf X}_{J(j,E)} \psi_j(s), \mbox{\hspace{5cm} (2)} \]

where the sums are taken over the nodes $ j $ of the element. $ J(j,E) $ represents the global node number of local node $ j $ in element $ E $ and $ U_J $ and $ {\bf X}_J $ represent the function value at and the position vector to global node $ J $, respectively. To ensure inter-element continuity of $ u $ and $ {\bf x} $, we constrain the nodal values and positions of the hanging nodes so that for every hanging node $ J $ we have

\[ U_J = \sum_{K} U_{K} \omega_{JK} \]

and

\[ {\bf X}_J = \sum_{K} {\bf X}_{K} \omega_{JK} \]

where the sum is taken over the hanging node's master nodes $ K $ and the $ \omega_{JK} $ are suitable weights. It is precisely this representation of the nodal positions and values that is implemented in Node::value(...) and Node::position(). [Note that different nodal values can have different hanging node constraints; e.g. in Taylor-Hood elements where the pressure and velocities are interpolated by linear and quadratic basis functions.]

For simply hanging nodes (e.g. Nodes 9, 10, 14 and 15 in the above sketch) the weights $ \omega_{JK} $ are determined as follows:

For multiply hanging nodes (Nodes 17 and 18 in the above sketch), the weights of the ultimate master nodes are determined recursively, e.g. $ \omega_{17\ 4} = \omega_{17\ 9} \times \omega_{9 \ 4}. $

For RefineableQuadMesh meshes, the above procedures are fully implemented. Executing the function

RefineableQElement<2>::setup_hanging_nodes()

for each element in the mesh establishes which of the elements' Nodes are hanging and determines their primary master nodes and weights. Furthermore, it pins the values at the hanging nodes – because the values are constrained, they cannot be unknowns in the problem. When this function has been executed for all elements in the mesh, the recursive elimination of hanging master nodes is performed by calling

RefineableMesh::complete_hanging_nodes();

As mentioned above, the possible occurrence of hanging nodes needs to be reflected in the element's equation numbering scheme and in the functions that compute the elemental Jacobian matrix because the element residuals now potentially depend on Nodes outside the element. Therefore, RefineableQElement<2>s must re-implement various FiniteElement member functions, e.g., by re-implementing the virtual functions

RefineableQElement<2>::get_residuals()

and

RefineableQElement<2>::get_jacobian(...)

and various others, as specified in the RefineableQElement<2> class.

In practice, we again distinguish between the "geometry" and the "maths" by writing a general RefineablePoissonEquations<DIM> class that inherits from PoissonEquations<DIM> and re-implements the appropriate member functions.

Mesh adaptation

We can now discuss the details of the mesh adaptation process for RefineableQuadMesh meshes, although the general procedure is, in fact, generic: Once a solution has been computed, Problem::adapt() loops over all refineable submeshes, and uses their error estimator functions to compute the vector of elemental errors. The vector of errors is then passed to the submeshes'

RefineableMesh::adapt(...)

function which performs the actual mesh adaptation and entails the following steps:

Mesh adaptation in domains with curved boundaries

The mesh refinement procedures outlined above are perfectly adequate for meshes in polygonal domains. In such meshes the generation of the new nodal positions and the transfer of the solution from the old to the new mesh can be performed by simple interpolation, using the `father' element's geometric mapping and shape functions. However, in problems with curvilinear mesh boundaries we must ensure that the refined meshes provide a progressively more accurate representation of the curvilinear domain boundary.

To facilitate these steps, we now introduce the concept of GeomObjects, Domains and MacroElements, which allow a convenient and generic representation of domains that are bounded by time-dependent and/or curvilinear boundaries.

Geometric Objects

Here are two examples of curvilinear boundaries that are frequently encountered in computations on moving domains:

The common feature of these two examples is that, in both cases, the boundary is represented by a parametrised position vector. The GeomObject base class provides an abstract interface for such situations. In its most basic form, a `geometric object' simply provides a parametrisation of its shape by implementing the GeomObject's pure virtual member function

GeomObject::position(xi,r)

which computes the position vector r at the coordinates specified by the vector xi. GeomObject also provides a large number of additional interfaces in the form of (broken) virtual functions. The most important of these is the time-dependent version of the above function

GeomObject::position(xi,t,r)

which computes the position vector r at the coordinates specified by the vector xi at the previous discrete time level t. We follow the usual convention that

By default, the virtual member function GeomObject::position(xi,t,r) calls the steady version GeomObject::position(xi,r), so it only needs to be overloaded for genuinely time-dependent geometric objects (by default, the code execution terminates if the time-dependent version is called for $ {\tt t} \ne 0; $ the function needs to be overloaded if this is not appropriate).

Further virtual member functions provide interfaces for the determination of the first and second spatial and temporal derivatives of the position vector, allowing the computation of boundary velocities, accelerations and curvatures, etc. These interfaces are provided as broken virtual functions, so they only need to be overloaded if/when the functionality is actually required in a particular application.

Typically, the shape of a GeomObject depends on a certain number of parameters (in the above examples, the radius $ R(t) $ of the cylinder and the displacements of the shell element, respectively) which can be unknowns in the problem. We therefore store these parameters as (geometric) Data, whose values can be pinned or free.

Domains and Macro Elements

For the purposes of mesh generation, we represent curvilinear domains by objects that are derived from the base class Domain and use GeomObjects to represent the curvilinear boundaries. For instance, the QuarterCircleSectorDomain sketched in the figure below is bounded by the geometric object Ellipse that parametrises the domain's curved boundary, shown in red.

Consider the coarse discretisation of the domain shown in the Fig. (a) and assume that element 2 (a four-node quad element from the QElement family) is scheduled for refinement into four son elements. As discussed above, by default the Nodes of the son elements are established/created as follows:

MacroElementSketch.gif
Illustration of mesh refinement with and without macro elements

For the element 0, this procedure would be perfectly adequate, as the domain boundary after refinement would (still) be captured exactly. However, when refining element 2 by this procedure, the new Node on the boundary is positioned on the straight line between the two boundary Nodes in the father element and not on the curved boundary itself, as shown in Fig. (b). Hence repeated mesh refinement will not lead to a better representation of the domain geometry and preclude convergence to the exact solution.

To overcome this problem, objects of type Domain provide a decomposition of the domain into a number of so-called MacroElements, as sketched in Fig. (c). Boundaries of the MacroElements are either given by (parts of) the (exact) curvilinear boundaries (as specified by the Domain's GeomObjects) or by (arbitrary and usually straight) internal edges (or, in 3D, faces). In 2D, MacroElements provide a mapping $ {\bf r}(S_0,S_1) $ from a local coordinate system $ (S_0, S_1) $, with $ S_i \in [-1..1]^2 $, to the points inside the MacroElement via their member function

MacroElement::macro_map(S,r)

The mapping needs to be chosen such that for $ S_1\in [-1..1]$ the position vector $ {\bf r}(S_0=-1,S_1) $ sweeps along the `southern' edge of the MacroElement, etc.; see Fig. (e). The form of the macro-element mapping is obviously not unique and it depends on the MacroElement's topology. The class QMacroElement<2> provides a mapping that is suitable for 2D quadrilateral elements and can therefore be used with RefineableQElement<2> s.

The constructors for objects of type Domain typically require a vector of pointers to the GeomObjects that define its boundaries. The Domain constructor then usually employs function pointers to the GeomObject::position(...) function to define the MacroElement boundaries. Once built, the MacroElements are accessible via the member function

Domain::macro_element_pt(...).

We illustrate the macro-element based mesh generation and adaptation process for the case of RefineableQuadMesh meshes. Assume that the domain is represented by an object of type Domain. We build a coarse initial mesh which contains as many RefineableQElement<2>s as there are QMacroElements in the Domain. We associate each RefineableQElement<2> with one of the QMacroElements by storing a pointer to the QMacroElement in FiniteElement::Macro_elem_pt. Next, we use the QMacroElements' macro map to determine the RefineableQElement<2>'s nodal positions such that, e.g., the RefineableQElement<2>'s SW node is placed at the QMacroElement 's local coordinates $(S_0,S_1)=(-1,-1)$, etc. The fraction of the QMacroElement that is spanned by each RefineableQElement<2> is represented by the maximum and minimum values of the QMacroElement's local coordinates; the RefineableQElement<2> constructor initially sets these values to the defaults of +1.0 and -1.0, respectively, indicating that the RefineableQElement<2> spans the entire area of the corresponding QMacroElement.

With this additional information, we modify the RefineableQElement<2>::build() process as follows:

Macro-element-based node updates in moving mesh problems

Once a Mesh is associated with a Domain, the function Mesh::node_update() updates the nodal positions in response to a change in the shape of the GeomObjects that define the Domain boundaries. [Note: This function updates all nodal positions first and then updates the positions of the hanging nodes to subject them to the hanging node constraints.] Alternative procedures for the update of the nodal positions in response to changes in the domain boundary are implemented in the AlgebraicMesh, SpineMesh and SolidMesh classes.

PDF file

A pdf version of this document is available.