In this document, we discuss the solution of a 2D Poisson problem using oomphlib's
powerful mesh adaptation routines:
in the fishshaped domain , with homogeneous Dirichlet boundary conditions

The sharp corners in the domain create singularities in the solution (its derivatives are unbounded) and so accurate results can only be obtained if we use a sufficiently fine discretisation. Implementing this by uniform mesh refinement would create a huge number of elements in the interior of the domain where the fine discretisation is not required.
To avoid this problem, oomphlib
provides mesh adaptation routines that automatically adapt the mesh, based on a posteriori error estimates. Regions in which an error estimator indicates that the solution is not resolved to the required accuracy are refined; automatic unrefinement is performed in regions where the discretisation is unnecessarily fine.
We provide a detailed discussion of the driver code fish_poisson.cc which illustrates a variety of mesh refinement procedures. [The alternative driver code fish_poisson_no_adapt.cc solves the same problem without mesh adaptation. Its structure is very similar to that in the 2D Poisson problem considered earlier. It is provided mainly to illustrate how easy to it is incorporate adaptivity into a Problem
.]
In the current example we demonstrate how to use existing refineable meshes and elements. Two further examples will demonstrate how easy it is to create refineable meshes in domains with polygonal boundaries and in domains with curvilinear boundaries.
The namespace ConstSourceForPoisson
only contains the constant source function .
The main code is very short and calls two functions that illustrate two different adaptation strategies:
is performed automatically until the solution satisfies the required error bounds (or until the maximum permitted number of adaptation steps has been reached).
We start by creating the Problem object, using the refineable equivalent of the QPoissonElement
– the RefineableQPoissonElement
, which is templated by the dimension and the number of nodes along the element's edges; the RefineableQPoissonElement<2,3>
is a ninenode (biquadratic) quad element.
After creating the DocInfo
object, we document the (default) adaptivity targets:
These include
These default parameters can be changed by the user; see Comments and Exercises.
The fullyadaptive solution of the problem is very simple. We simply pass the maximum number of adaptations to the Newton solver and document the results. Done!
To allow the user more control over the mesh adaptation process, oomphlib
provides a number of functions that perform individual adaptation steps without recomputing the solution immediately. This allows the user to
The second driver function illustrates some of these functions. We start by setting up the problem, create the DocInfo
object and document the adaptivity targets, exactly as before:
Next, we solve the problem on the original, very coarse mesh and document the result:
We know that the result is unlikely to be very accurate, so we apply three levels of uniform refinement, increasing the number of elements from 4 to 256, and recompute:
The solution looks much smoother but we suspect that the corner regions are still underresolved. Therefore, we call the Problem::adapt()
function which computes an error estimate for all elements and automatically performs a single mesh adaptation (refinement/unrefinement) step. If this adaptation changes the mesh, we recompute the solution, using the "normal" Newton solver without automatic adaptation. We document the solution and continue the adaptation cycle until Problem::adapt()
ceases to change the mesh:
The progress of the adaptation is illustrated in the animated gif at the beginning of this document. The first frame displays the solution on the original fourelement mesh; the next frame shows the solution on the uniformly refined mesh; the final two frames show the progress of the subsequent, errorestimatedriven mesh adaptation.
The problem class is virtually identical to that used in the 2D Poisson problem without mesh refinement. In the present problem, we leave the function Problem::actions_before_newton_solve()
empty because the boundary conditions do not change. The function RefineableFishPoissonProblem::mesh_pt()
overloads the (virtual) function Problem::mesh_pt()
since it returns a pointer to a generic Mesh
object, rather than a pointer to the specific mesh used in this problem. This avoids explicit recasts in the rest of the code where member functions of the specific mesh need to be accessed.
We start by creating the mesh, using oomphlib's
RefineableFishMesh
object:
Next, we create an error estimator for the problem. The Z2ErrorEstimator
is based on Zhu and Zienkiewicz's flux recovery technique and can be used with all elements that are derived from the ElementWithZ2ErrorEstimator
base class (or with functions that implement the pure virtual functions that are defined in this class) – the RefineableQPoissonElement
is an element of this type.
Next we pin the nodal values on all boundaries, apply the homogeneous Dirichlet boundary conditions, pass the pointer to the source function to the elements, and set up the equation numbering scheme.
The postprocessing routine writes the computed result to an output file, labeled with the identifiers specified in the DocInfo
object.
The purpose of this example was to provide a highlevel overview of oomphlib's
mesh adaptation procedures. We demonstrated that the implementation of full adaptivity only required us to
FishMesh
and the QPoissonElement
objects by their refineable equivalents, RefineableFishMesh
and RefineableQPoissonElement
, respectively(Compare the Problem specification for the current problem to that of its nonrefineable equivalent, contained in the alternative driver code fish_poisson_no_adapt.cc.)
Since most of the "hard work" involved in the mesh adaptation is "hidden" from the user, we highlight some important aspects of the procedure:
The Problem::adapt()
function automatically determines the correct boundary conditions for newly created nodes on the Mesh boundary; it automatically updates the equation numbering scheme, and interpolates the solution from the original mesh onto the adapted mesh. This is important in nonlinear problems where the provision of a good initial guess for the Newton iteration is vital; and in timedependent problems where the solution at one timestep provides initial conditions for the next one. See the discussion of the adaptive solution of the unsteady heat equation for more details. Furthermore, the source function pointers are automatically passed to an element's four "son" elements when the element is subdivided. This allows the adaptation to proceed completely automatically, without any intervention by the "user". On return from Problem::adapt()
the problem can immediately be resolved.
In some special cases, certain actions may need to be performed before or after the mesh adaptation (e.g. if flux boundary conditions are applied by FaceElements
; this is explained in another example). To ensure that these steps are performed when the adaptation is controlled by the "blackbox" adaptive Newton solver, the Problem
class provides the two empty virtual functions
and
which are called automatically before and after the adaptation. The "user" can overload these in his/her specific Problem
class to implement such actions.
The mesh adaptation not only increases the number of elements but also produces a more accurate representation of the curvilinear domain boundary – new boundary nodes are placed exactly onto the analyticallydefined, curvilinear boundary, rather than on the boundaries of the "father" element, which only provides an approximate representation of the exact domain boundary. This is achieved by employing a MacroElementbased
representation of the Domain
– we will discuss this in more detail in another example.
Many adaptation routines in the Problem
class have equivalents in the RefineableMesh
class. It is important to appreciate the important differences between them: If adaptation is performed at the Problem
level, the adapted Problem
is fully functional, i.e. boundary conditions will have been assigned for newly created nodes on the mesh boundary, the equation numbering scheme will have been updated, etc. The adapted Problem
can therefore be resolved immediately. Conversely, if a mesh is refined directly, using the member functions of the RefineableMesh
class, many of these additional tasks need to be performed "by hand" before the adapted Problem
can be resolved.
To familiarise yourself with oomphlib's
mesh adaptation procedures we suggest the following exercises:
Problem::adapt()
does indeed interpolate the solution from the coarse mesh to the fine mesh – call Problem::doc_solution(...)
before and after its execution. Problem::refine_uniformly()
function has a counterpart Problem::unrefine_uniformly()
. Why does this function not simply unrefine every single element in the mesh? Explore the action of Problem::unrefine_uniformly()
by plotting the solution before and after a few executions of this function. Problem::refine_selected_elements(...)
.A pdf version of this document is available.