Inline mesh generation and adaptation based on Triangle

In this document we demonstrate how to generate unstructured 2D meshes for oomph-lib, using Jonathan Shewchuk's Triangle library from within an oomph-lib driver code. This complements the discussion in another tutorial where we illustrated how to build oomph-lib meshes using the output generated when Triangle is used to create the mesh externally.

In addition, we demonstrate


Here we restrict ourselves to the solution of problems in fixed domains. Other tutorials discuss more advanced applications involving moving meshes, e.g.

Overview of mesh generation procedures

Meshes with polygonal boundaries

If the domain has polygonal outer (and possibly internal) boundaries, the mesh generation process follows a similar pattern to that required in the external use of Triangle , but using oomph-lib classes to represent the required data. We define the polygonal boundaries using TriangleMeshPolyLines, which represent distinct mesh boundaries and are defined by a series of vertex coordinates. Multiple TriangleMeshPolyLines can be combined to define (closed) TriangleMeshPolygons. These are then used to create an unstructured mesh with a pre-determined target area for all elements, via an interface to Triangle .

One particular feature of oomph-lib's interface to Triangle is that each closed TriangleMeshPolygon must contain at least two distinct TriangleMeshPolyLines, each with its own boundary ID. This is because every Node can only have a single-valued boundary coordinate but if the boundary is closed there would be a discontinuity in the coordinate value. (For example, a standard representation would have $ \zeta = [0,2\pi] $ , but the node at $\zeta=0$ should also have $ \zeta=2\pi $ and this is not possible.) The sketch below shows a representative domain as well as two legal and one illegal representations of the domain boundaries. Note that the boundaries can be enumerated in an arbitrary fashion.

Sketch of a polygonal domain (top) and two legal (bottom left and middle) and one illegal (bottom right) representations of the boundaries in terms of TriangleMeshPolyLines.

Meshes with curvilinear boundaries

It is also possible to discretise domains with curvilinear boundaries as shown in the sketch below. Assuming that each curvilinear boundary is represented by a GeomObject that specifies the position vector $ {\bf R}(\zeta) $ to a point on the curvilinear boundary as a function of some surface coordinate $ \zeta $, we split each closed boundary into (at least) two distinct TriangleMeshCurviLines – the curvilinear equivalents of TriangleMeshPolyLines. Each TriangleMeshCurviLine is constructed from a pointer to the GeomObject and the start and end values of the boundary coordinate $ \zeta $ along the relevant part of the curvilinear boundary. The TriangleMeshCurviLines are then combined to a TriangleMeshClosedCurve – the general closed curve.

Sketch of a domain bounded by a curvilinear boundary, containing two holes with curvilinear (hole 1) and polygonal (hole 2) boundaries, respectively.

The mesh is then created in a two-stage process: All TriangleMeshCurviLines are sampled at a certain number of points (specified by the user in the constructor) to create the vertices for a polyline representation of the boundary. This polygonal representation of the boundaries is used by Triangle to generate the mesh. Finally, nodes on curvilinear boundaries are "snapped" onto the actual curvilinear boundary.

Overview of mesh adaptation methodology

The methodology employed to adapt oomph-lib's unstructured meshes differs from that used for structured meshes. Specifically, rather than sub-dividing elements in which the error estimate exceeds a threshold and merging elements in which the solution is "too accurate", we completely re-generate the mesh and project the solution from the old to the new mesh. This is because we originally developed the underlying methodology to solve free-boundary problems in which the domain undergoes such large deformations that re-meshing is required. The ability to adjust the element sizes guided by spatial error estimates when re-meshing the domain is a simple fringe benefit.

A number of issues are important:

We recommend using the Problem::actions_after_adapt() function to re-assign boundary conditions and to complete the build of all elements after the adaptation.

Apart from these issues, the user interfaces to the mesh adaptation functions are exactly the same as for structured meshes. Specifically, it is possible to specify maximum and minimum element sizes and target values for the error such that the mesh is refined in regions where the error estimate is too large'' and unrefined where it is too small''.

Typically, the most computationally expensive stage of the mesh regeneration procedure is the multi-domain setup procedure which identifies corresponding points in the old and new meshes. In "cheap" problems, such as the Poisson problem discussed below, the cost of the mesh regeneration can exceed the cost of the subsequent solve, but in most "hard" problems (such as the ones listed at the beginning of this tutorial) the cost of the mesh regeneration is modest (and, in the case of large-displacement free-boundary problems, unavoidable anyway).

An example: The adaptive solution of Poisson's equations on an unstructured mesh

As an example we consider the adaptive solution of Poisson's equation

\[ \frac{\partial^2 u}{\partial x_i^2} = f(x_1,x_2) \]

in a circular domain that contains an elliptical and a polygonal hole. As in many previous examples, we apply Dirichlet boundary conditions on all domain boundaries and choose the boundary values and the source function $ f(x_1,x_2) $ such that the exact solution of the problem is given by

\[ u(x_1,x_2) = \tanh(1-\alpha(x_1 \tan\Phi - x_2)), \]

which approaches a step function, oriented at an angle $ \Phi $ within the $ (x_1,x_2) $ plane, as $ \alpha $ becomes large. The figure below shows contour plots of the solution for $ \alpha = 5 $ for various angles $ \Phi $. It illustrates how the mesh adaptation adjusts the mesh such the smallest elements are located in the region where the solution undergoes rapid change.

Contour plot of the solution.

Global parameters and functions

Following our usual practice, we use a namespace to define the source function and the exact solution.

//===== start_of_namespace=============================================
/// Namespace for exact solution for Poisson equation with "sharp step"
namespace TanhSolnForPoisson
/// Parameter for steepness of "step"
double Alpha=5.0;
/// Parameter for angle Phi of "step"
double TanPhi=0.0;
/// Exact solution as a Vector
void get_exact_u(const Vector<double>& x, Vector<double>& u)
/// Source function required to make the solution above an exact solution
void get_source(const Vector<double>& x, double& source)
source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
/// \short Zero function -- used to compute norm of the computed solution by
/// computing the norm of the error when compared against this.
void zero(const Vector<double>& x, Vector<double>& u)
} // end of namespace

The driver code

We start by processing command line arguments which allow us to run the code in self-test mode and build the problem with "projectable" six-noded triangular Poisson elements.

///Driver code for demo of inline triangle mesh generation
int main(int argc, char **argv)
// Store command line arguments
// Define possible command line arguments and parse the ones that
// were actually specified
// Validation?
// Parse command line
// Doc what has actually been specified on the command line
// Create problem

We then perform a parameter study, solving the problem for various orientations of the "step" and allowing a certain number of spatial adaptations per solve. (If the code is run in self-test mode, we perform fewer steps and allow for less adaptation to speed up the computation.)

// Loop over orientation of step
unsigned nstep=5;
if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
for (unsigned i=0;i<nstep;i++)
// Solve with spatial adaptation
unsigned max_adapt=3;
if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
// Doc the solution
std::stringstream comment_stream;
comment_stream << "Solution for tan(phi) = " << TanhSolnForPoisson::TanPhi;
// Rotate orientation of solution
} //End of main

The problem class

The problem class contains the usual member functions. As discussed above, the boundary conditions and the source function have to be re-specified after every mesh adaptation since the adapted mesh contains completely new elements. This is done in the function complete_problem_setup(), discussed below, which is called from the Problem constructor and from actions_after_adapt(). We re-assign the Dirichlet boundary conditions in actions_before_newton_solve(), using a second helper function apply_boundary_conditions():

/// Class definition
template<class ELEMENT>
class UnstructuredPoissonProblem : public virtual Problem
/// Constructor
/// Destructor
/// Actions before adapt. Empty
/// \short Actions after adapt:
/// Setup the problem again -- remember that the mesh has been
/// completely rebuilt and its element's don't have any
/// pointers to source fcts etc. yet
/// Update after solve (empty)
/// Update the problem specs before solve: Re-apply boundary conditons
/// Doc the solution
void doc_solution(const std::string& comment="");
/// Doc info object for labeling output
DocInfo Doc_info;
/// Helper function to apply boundary conditions
/// \short Helper function to (re-)set boundary condition
/// and complete the build of all elements
/// Pointers to specific mesh
RefineableTriangleMesh<ELEMENT>* My_mesh_pt;
/// Trace file to document norm of solution
ofstream Trace_file;
}; // end_of_problem_class

The Problem constructor

Most of the problem constructor is concerned with the specification of the mesh boundaries. We start by generating a GeomObject that describes the circular outer boundary of the domain.

/// Constructor
template<class ELEMENT>
// Intrinsic coordinate along GeomObject
Vector<double> zeta(1);
// Position vector on GeomObject
Vector<double> posn(2);
// Ellipse defining the outer boundary
double x_center = 0.0;
double y_center = 0.0;
double A = 1.0;
double B = 1.0;
Ellipse * outer_boundary_ellipse_pt = new Ellipse(A,B);

This GeomObject is now used to describe the outer boundary in terms of a TriangleMeshClosedCurve object, a base class which can represent polygonal and curvilinear boundaries. We start by providing a pointer to this (yet-to-be-built) object.

// Pointer to the closed curve that defines the outer boundary
TriangleMeshClosedCurve* closed_curve_pt=0;

As discussed above, the closed outer boundary must be broken up into (at least) two distinct sub-boundaries to allow oomph-lib to automatically refine and setup boundary coordinates. We therefore create two TriangleMeshCurviLines, specifying

// Provide storage for pointers to the two parts of the curvilinear boundary
Vector<TriangleMeshCurveSection*> outer_curvilinear_boundary_pt(2);

We choose five boundary segments for the first TriangleMeshCurviLine which represents the upper half of the boundary which we label as boundary 0,

// First bit
double zeta_start=0.0;
double zeta_end=MathematicalConstants::Pi;
unsigned nsegment=5;
unsigned boundary_id=0;
outer_curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(

and eight segments for the lower half which we label as boundary 1:

// Second bit
outer_curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(

We then combine the two TriangleMeshCurviLines to a TriangleMeshClosedCurve which describes the outer boundary.

// Combine to curvilinear boundary and define the
// outer boundary
new TriangleMeshClosedCurve(outer_curvilinear_boundary_pt);

Next we deal with the two inner (hole) boundaries

// Now build the holes
Vector<TriangleMeshClosedCurve*> hole_pt(2);

The first hole is a polygon whose 12 vertices we distribute along a circle of radius $ 0.1 $ , centred at $ (x_1,x_2)=(0,0.5). $ As above, we break the closed boundary into two distinct sub-boundaries – this time represented by TriangleMeshPolyLines:

// Build polygonal hole
// Build first hole: A circle
x_center = 0.0;
y_center = 0.5;
A = 0.1;
B = 0.1;
Ellipse* polygon_ellipse_pt=new Ellipse(A,B);
// Number of segments defining upper and lower half of the hole
unsigned n_seg = 6;
double unit_zeta = MathematicalConstants::Pi/double(n_seg);
// This hole is bounded by two distinct boundaries, each
// represented by its own polyline
Vector<TriangleMeshCurveSection*> hole_polyline_pt(2);

We create the vertex coordinates for the upper half of the polygonal hole,

// First boundary of polygonal hole
// Vertex coordinates
Vector<Vector<double> > bound_hole(n_seg+1);
for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
// Resize the vector
// Get the coordinates

and build the TriangleMeshPolyLine, specifying a boundary ID:

// Specify the hole boundary id
unsigned boundary_id=2;
// Build the 1st hole polyline
hole_polyline_pt[0] = new TriangleMeshPolyLine(bound_hole,boundary_id);

We repeat the exercise for the lower half which we turn into boundary 4:

// Second boundary of polygonal hole
for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
// Resize the vector
// Get the coordinates
// Specify the hole boundary id
// Build the 2nd hole polyline
hole_polyline_pt[1] = new TriangleMeshPolyLine(bound_hole,boundary_id);

Finally, we build the polygonal hole itself, specifying its constituent TriangleMeshPolyLines and the coordinate of a point inside the hole, which is required by Triangle:

// Build the polygonal hole
// Inner hole center coordinates
Vector<double> hole_center(2);
hole_pt[0] = new TriangleMeshPolygon(hole_polyline_pt, hole_center);

The construction of the second, curvilinear internal boundary (an ellipse centred at the origin) is virtually identical to the steps taken for the construction of the outer boundary, apart from the fact that, as an internal boundary, it again requires the specification of a point inside the hole.

// Build curvilinear hole
// Build second hole: Another ellipse
A = 0.2;
B = 0.1;
Ellipse* ellipse_pt=new Ellipse(A,B);
// Build the two parts of the curvilinear boundary
Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt(2);
// First part of curvilinear boundary
double zeta_start=0.0;
double zeta_end=MathematicalConstants::Pi;
unsigned nsegment=10;
curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
// Second part of curvilinear boundary
curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
// Combine to hole
Vector<double> hole_coords(2);
Vector<TriangleMeshClosedCurve*> curvilinear_hole_pt(1);
new TriangleMeshClosedCurve(curvilinear_boundary_pt,

Construct the mesh

To facilitate the construction of the mesh TriangleMesh object we use the object TriangleMeshParameters. The only necessary argument for creating this object is the outer boundary. The definition of holes, internal boundaries and regions is explained in another tutorial. The object can also be used to control whether additional refinement may be performed on the mesh boundaries. The default behaviour is that such refinement will occur so that the highest quality mesh is obtained. In some cases, e.g. periodic boundary conditions, you may wish to ensure that each input boundary segment corresponds to a single element edge in the final mesh. This can be achieved for the outer boundary by calling TriangleMeshParameters::disable_boundary_refinement(). However, Triangle will still add additional points to any internal boundaries unless the additional function TriangleMeshParameters::disable_internal_boundary_refinement() is also called. The functions TriangleMeshParameters::enable_boundary_refinement() and TriangleMeshParameters::enable_internal_boundary_refinement() can be used to return to the default behaviour. Note that it is not currently possible to suppress refinement on the internal boundaries, but refine the outer boundary.

We can specify a target for the element sizes and pass it to the TriangleMeshParameters object as showed next:

// Now build the mesh
// Use the TriangleMeshParameters object for helping on the manage of the
// TriangleMesh parameters
TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
// Specify the closed curve using the TriangleMeshParameters object
triangle_mesh_parameters.internal_closed_curve_pt() = hole_pt;
// Specify the maximum area element
double uniform_element_area=0.2;
triangle_mesh_parameters.element_area() = uniform_element_area;
// Create the mesh
// Store as the problem's one and only mesh

We specify a spatial error estimator and limit the maximum and minimum element sizes,

// Set error estimator for bulk mesh
Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
// Set element size limits

before completing the problem setup (see below) and assigning the equation numbers.

// Set boundary condition and complete the build of all elements
// Open trace file
char filename[100];
// Setup equation numbering scheme
oomph_info <<"Number of equations: "
<< this->assign_eqn_numbers() << std::endl;
} // end_of_constructor

Completing the problem setup

As discussed above, the helper function complete_problem_setup() starts by (re-)applying the boundary conditions by pinning the nodal values on all mesh boundaries,

/// Set boundary condition exactly, and complete the build of
/// all elements
template<class ELEMENT>
// Set the boundary conditions for problem: All nodes are
// free by default -- just pin the ones that have Dirichlet conditions
// here.
unsigned nbound=My_mesh_pt->nboundary();
for(unsigned ibound=0;ibound<nbound;ibound++)
unsigned num_nod=My_mesh_pt->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Get node
Node* nod_pt=My_mesh_pt->boundary_node_pt(ibound,inod);
// Pin one-and-only unknown value
} // end loop over boundaries

specifies the source function pointer for all elements,

// Complete the build of all elements so they are fully functional
unsigned n_element = My_mesh_pt->nelement();
for(unsigned e=0;e<n_element;e++)
// Upcast from GeneralisedElement to the present element
ELEMENT* el_pt = dynamic_cast<ELEMENT*>(My_mesh_pt->element_pt(e));
//Set the source function pointer
el_pt->source_fct_pt() = &TanhSolnForPoisson::get_source;

and then re-sets the boundary values:

// Re-apply Dirichlet boundary conditions (projection ignores
// boundary conditions!)

Assigning the boundary values

The function apply_boundary_conditions() does exactly what is says: It loops over all boundary nodes and assigns the value according the exact solution specified in the namespace TanhSolnForPoisson.

/// Helper function to apply boundary conditions
template<class ELEMENT>
// Loop over all boundary nodes
unsigned nbound=this->My_mesh_pt->nboundary();
for(unsigned ibound=0;ibound<nbound;ibound++)
unsigned num_nod=this->My_mesh_pt->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Get node
Node* nod_pt=this->My_mesh_pt->boundary_node_pt(ibound,inod);
// Extract nodal coordinates from node:
Vector<double> x(2);
// Compute the value of the exact solution at the nodal point
Vector<double> u(1);
// Assign the value to the one (and only) nodal value at this node
} // end set bc


We compare the computed solution against the exact solution:

/// Doc the solution
template<class ELEMENT>
std::string& comment)
ofstream some_file;
char filename[100];
// Number of plot points
unsigned npts;
some_file << "TEXT X = 22, Y = 92, CS=FRAME T = \""
<< comment << "\"\n";
// Output exact solution
// Output boundaries
// Doc error and return of the square of the L2 error
double error,norm,dummy_error,zero_norm;
// Doc L2 error and norm of solution
oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
oomph_info << "Norm of exact solution: " << sqrt(norm) << std::endl;
oomph_info << "Norm of computed solution: " << sqrt(dummy_error) << std::endl;
Trace_file << sqrt(norm) << " " << sqrt(dummy_error) << std::endl;
// Increment the doc_info number
} // end of doc

Comments and Exercises

Upgrading elements to become "projectable"

As discussed above, a key step in the "mesh-adaptation-by-mesh-regeneration" procedure is the projection of the solution from the old to the new mesh. The ability to perform this projection fully-automatically during the mesh adaptation requires the elements to be "wrapped" in the templated ProjectableElement<ELEMENT> class. This class is derived from the ProjectableElementBase base class which specifies a number of pure virtual functions that must be specified for each specific element type. For the Poisson problem considered here these functions are already provided in the ProjectablePoissonElement<ELEMENT> class. Similar wrappers exist for many other equations. If you want to "upgrade" your own elements to become projectable, inspect the prototypes for the relevant pure virtual functions which are defined in


while the specific implementation for the projectable Poisson elements is provided in


Of course, you can avoid the additional work by dispensing with adaptivity and simply generating a sufficiently fine (uniform) mesh. This is illustrated in the alternative driver code


How many vertices should I use to sample my curvilinear boundary?

As discussed above, meshes with curvilinear boundaries are created in a two-stage process. Initially, Triangle generates a mesh with polygonal boundaries using a user-specified number of vertices that are evenly distributed along the relevant GeomObject. The nodes on that boundary are then "snapped" onto the actual curvilinear boundary in a post-processing step. The decision of how many vertices to choose involves a compromise between two conflicting demands:


  1. Change the outer curvilinear boundary to a polygonal boundary (you can cheat – the relevant code is already contained in the driver code but it's "hidden" with ifdefs; this code also documents the re-distribution of straight-line segments between different polylines, a capability that is important in certain free-boundary problems).

  2. Vary the number of vertices used for the initial polygonal representation of the curvilinear hole to establish what number is required to avoid the inversion of elements during the "snap-nodes-to-the-curvilinear-boundary" phase.

  3. Create "projectable" advection diffusion elements to solve the the advection diffusion problem discussed in another tutorial, using spatial adaptation on an unstructured mesh.

Source files for this tutorial

Appendix A: Generalization for polylines and curvilines

The objects, TriangleMeshPolyLine and TriangleMeshCurviLine inherit the properties of a more general representation called TriangleMeshCurveSection. This allows one to define more general boundaries as a combination of TriangleMeshPolyLines and TriangleMeshCurviLines. Therefore, if we want to define a more general closed curve use the TriangleMeshClosedCurve object. You may notice the use of this object through the example code when defining the outer boundary.

For the interested reader, the class diagram showing the hierarchy of the mentioned objects is showed on the next figure.

The hierarchy of the TriangleMesh objects.

Note the TriangleMeshOpenCurve object, which allows to define internal boundaries on the domain, explained on another tutorial

PDF file

A pdf version of this document is available.