Creation of internal boundaries for Triangle based meshes

In this document we show how to add internal boundaries and create connections among them in 2D unstructured meshes with `oomph-lib`

. The mesh generation procedure is performed by using Jonathan Shewchuk's Triangle. This tutorial complements the discussion in another tutorial where we illustrated how to build `oomph-lib`

meshes using the `TriangleMeshPolygon`

and `TriangleMeshClosedCurve`

objects.

Additionally, we show via examples

- how to connect
`TriangleMeshCurveSection`

objects - how to add addition holes and extra regions, i.e. regions not defined via
`TriangleMeshPolygon`

or`TriangleMeshClosedCurve`

objects

An internal boundary is defined in `oomph-lib`

by the object `TriangleMeshOpenCurve`

. We can create and add as many internal boundaries as we need for defining our domain, see figure below.

The `TriangleMeshOpenCurve`

object allows the creation of complex internal boundaries as a combination of `TriangleMeshPolyLine`

and/or `TriangleMeshCurviLine`

objects. In order to create a `TriangleMeshOpenCurve`

it is necessary to create a `Vector`

of `TriangleMeshCurveSection`

objects that define the `TriangleMeshOpenCurve`

, see figure below.

Any internal boundary can be connected (via its ends) to any other boundary currently in the domain.

It is important to note that there are four different types of connections:

- Connection between a
`TriangleMeshPolyLine`

and a`TriangleMeshPolyLine`

- Connection between a
`TriangleMeshCurviLine`

and a`TriangleMeshPolyLine`

- Connection between a
`TriangleMeshPolyLine`

and a`TriangleMeshCurviLine`

- Connection between a
`TriangleMeshCurviLine`

and a`TriangleMeshCurviLine`

The next figure shows an sketch of the different types of connections.

Connections to `TriangleMeshPolyLine`

can only be made by using existing vertices of the destination `TriangleMeshPolyLine`

. Conversely, connections to `TriangleMeshCurviLine`

can be made by specifying the position along the `TriangleMeshCurviLine`

. The connection point is identified by the `TriangleMeshCurviLine's`

intrinsic coordinates. The figure below shows the case when a connection with a `TriangleMeshPolyLine`

is possible.

We will show how to add internal boundaries and how to create connections among them by using the same example as in another tutorial.

A sketch of the domain is shown on the next figure. The domain has seven boundaries, two of which define the outer boundary and the other five define internal boundaries.

The following code is extracted from the Problem constructor and shows how to create the domain/mesh.

The first task is to define the outer boundary which we represent by a `TriangleMeshClosedCurve`

comprising two `TriangleMeshPolyLines`

.

// Begin - Outer boundary

// Create storage for the vertices coordinates that define each boundary

const unsigned num_vertices_b0 = 4;

Vector<Vector <double> > vertices(num_vertices_b0);

vertices[0].resize(2);

vertices[1].resize(2);

vertices[2].resize(2);

vertices[3].resize(2);

// Building the outer boundary (vertices)

vertices[0][0] = 0;

vertices[0][1] = 0;

vertices[1][0] = 0;

vertices[1][1] = 1;

vertices[2][0] = 0;

vertices[2][1] = 3;

vertices[3][0] = 3;

vertices[3][1] = 3;

// The outer boundary is represented by two TriangleMeshPolyLine objects,

// this is the first one

boundary_id = 0;

TriangleMeshPolyLine *boundary0_pt =

new TriangleMeshPolyLine(vertices, boundary_id);

// Create storage for the vertices coordinates that define each boundary

const unsigned num_vertices_b1 = 5;

vertices.resize(num_vertices_b1);

vertices[0].resize(2);

vertices[1].resize(2);

vertices[2].resize(2);

vertices[3].resize(2);

vertices[4].resize(2);

// More vertices for the outer boundary

vertices[0][0] = 3;

vertices[0][1] = 3;

vertices[1][0] = 3;

vertices[1][1] = 0.5;

vertices[2][0] = 3;

vertices[2][1] = 0;

vertices[3][0] = 1;

vertices[3][1] = 0;

vertices[4][0] = 0;

vertices[4][1] = 0;

// The outer boundary is represented by two TriangleMeshPolyLine objects,

// this is the second one

boundary_id = 1;

TriangleMeshPolyLine *boundary1_pt =

new TriangleMeshPolyLine(vertices, boundary_id);

// A vector for storing the outer boundary representation

Vector<TriangleMeshCurveSection*> outer_boundary_polyLines_pt(2);

outer_boundary_polyLines_pt[0] = boundary0_pt;

outer_boundary_polyLines_pt[1] = boundary1_pt;

TriangleMeshClosedCurve *outer_boundary_pt =

new TriangleMeshPolygon(outer_boundary_polyLines_pt);

// End - Outer boundary

Note that we defined additional vertices along the two boundaries to allow the creation of connections to the outer boundary, see figure below.

We create a `Vector`

of five `TriangleMeshOpenCurves`

which will represent the internal boundaries.

// Internal open boundaries

// Total number of open curves in the domain

unsigned n_open_curves = 5;

// We want internal open curves

Vector<TriangleMeshOpenCurve *> inner_open_boundaries_pt(n_open_curves);

** Simple internal boundaries without connections **

Our domain has two unconnected internal boundaries: a straight line, represented by a `TriangleMeshPolyLine`

and a curve line represented by a `TriangleMeshCurviLine`

.

// We start by creating the internal boundaries

// The boundary 2 is defined by its two vertices

// Open curve 1

vertices.resize(2);

vertices[0].resize(2);

vertices[1].resize(2);

vertices[0][0] = 0.5;

vertices[0][1] = 2.0;

vertices[1][0] = 0.5;

vertices[1][1] = 2.5;

boundary_id = 2;

TriangleMeshPolyLine *boundary2_pt =

new TriangleMeshPolyLine(vertices, boundary_id);

// Each internal open curve is defined by a vector of TriangleMeshCurveSection,

// on this example we only need one curve section for each internal boundary

Vector<TriangleMeshCurveSection *> internal_curve_section1_pt(1);

internal_curve_section1_pt[0] = boundary2_pt;

// The open curve that define this boundary is composed of just one

// curve section

inner_open_boundaries_pt[0] =

new TriangleMeshOpenCurve(internal_curve_section1_pt);

// *********************************************************************

// We define the curved boundary as a TriangleMeshCurviline

// Open curve 2

double x_centre = 2.0;

double y_centre = 2.0;

double r_circle = 0.5;

Circle * boundary_circle1_pt = new Circle(x_centre, y_centre, r_circle);

// Number of segments used for representing the curve boundary

unsigned n_segments = 20;

// The intrinsic coordinates for the beginning and end of the curve

double s_start = 0.0;

double s_end = MathematicalConstants::Pi;

boundary_id = 3;

TriangleMeshCurviLine *boundary3_pt =

new TriangleMeshCurviLine(boundary_circle1_pt,

s_start,

s_end,

n_segments,

boundary_id);

Vector<TriangleMeshCurveSection *> internal_curve_section2_pt(1);

internal_curve_section2_pt[0] = boundary3_pt;

// The open curve that define this boundary is composed of just one

// curve section

inner_open_boundaries_pt[1] =

new TriangleMeshOpenCurve(internal_curve_section2_pt);

** Adding an internal boundary connected to an outer boundary **

We need to add three more internal boundaries which have connections with the outer boundary and with some internal boundaries; see the sketch of the domain.

In order to connect an internal boundary with any other boundary in the domain it is necessary to specify the following:

- Is the initial or the final end of the boundary the one that we want to connect?
- Do we want to connect the boundary with a
`TriangleMeshPolyLine`

or with a`TriangleMeshCurviLine`

? - If we want to connect with a
`TriangleMeshPolyLine`

then we need to identify the corresponding vertex number to which we want to connect. If we want to connect with a`TriangleMeshCurviLine`

we need to identify the corresponding intrinsic coordinate on the curve for connection.

Consider the creation of the internal boundary 4 shown on the sketch of the domain. This is an example of connecting a `TriangleMeshCurviLine`

with a `TriangleMeshPolyLine`

(the outer boundary). In this case both ends of the `TriangleMeshCurviLine`

are connected so we need to use the interface for connecting the initial and the final ends of the boundary to a `TriangleMeshPolyLine`

. There is only one thing left to identify, the vertex number of the destination boundary to which we want to connect. By looking at figure 1.5 we see that these values are vertex #3 for the initial end and vertex #1 for the final end. Once identified, we perform the connection as follow:

- Create the internal boundary as usual

// Open curve 3

x_centre = 0.0;

y_centre = 0.0;

r_circle = 1.0;

Circle * boundary_circle2_pt = new Circle(x_centre, y_centre, r_circle);

// Number of segments used for representing the curve boundary

n_segments = 20;

// The intrinsic coordinates for the beginning and end of the curve

s_start = 0.0;

s_end = MathematicalConstants::Pi / 2.0;

// State the vertex number for connection on the destination

// boundaries

unsigned vertex_to_connect_initial = 3;

unsigned vertex_to_connect_final = 1;

boundary_id = 4;

TriangleMeshCurviLine *boundary4_pt =

new TriangleMeshCurviLine(boundary_circle2_pt,

s_start,

s_end,

n_segments,

boundary_id);

- Specify connections using the destination boundary and the vertex number

// Do the connection with the destination boundary, in this case

// the connection is done with the outer boundary

boundary4_pt->connect_initial_vertex_to_polyline(

boundary1_pt,

vertex_to_connect_initial);

// Do the connection with the destination boundary, in this case

// the connection is done with the outer boundary

boundary4_pt->connect_final_vertex_to_polyline(

boundary0_pt,

vertex_to_connect_final);

- Finally, create the
`TriangleMeshOpenCurve`

object as usual

inner_open_boundaries_pt[2] =

new TriangleMeshOpenCurve(internal_curve_section3_pt);

** Adding an internal boundary connected with another internal boundary **

Boundary 5, a straight line connected to the outer boundary and to an internal boundary (the one that we created on the previous step).

The initial end of the boundary is connected to the outer boundary therefore we need to identify the vertex number to which it is connected; by looking at the sketch of the domain we see that the destination vertex number is 1. The final end is connected with a `TriangleMeshCurviLine`

therefore we need to specify the intrinsic coordinate in the object where we want to connect (the intrinsic coordinate is the arc-length along the circle).

// This boundary is connected to the outer boundary on the initial end

// and to an internal boundary on the final end

// Open curve 4

vertices.resize(3);

vertices[0].resize(2);

vertices[1].resize(2);

vertices[2].resize(2);

// We need to specify the vertices for the boundary on the first end

// The connection method performs a checking on the matching of the

// vertices used for connection

vertices[0][0] = 3.0;

vertices[0][1] = 0.5;

// These values are necessary when we subsequently perform the connection

// of the last boundary

vertices[1][0] = (3.0 + sqrt(3.0)) * 0.5;

vertices[1][1] = 0.5;

// We need to specify the vertices for the boundary on the last end

// The connection method performs a checking on the matching of the

// vertices used for connection

vertices[2][0] = sqrt(3.0) * 0.5;

vertices[2][1] = 0.5;

// State the vertex number for connection with the destination

// boundary

vertex_to_connect_initial = 1;

// State the s value for connection with the destination boundary

double s_connection_final = atan2(0.5, sqrt(3.0) * 0.5);

boundary_id = 5;

TriangleMeshPolyLine *boundary5_pt =

new TriangleMeshPolyLine(vertices, boundary_id);

// Do the connection with the destination boundary, in this case

// the connection is done with the outer boundary

boundary5_pt->connect_initial_vertex_to_polyline(

boundary1_pt,

vertex_to_connect_initial);

// Do the connection with the destination boundary, in this case

// the connection is done with an internal boundary

boundary5_pt->connect_final_vertex_to_curviline(

boundary4_pt,

s_connection_final);

Vector<TriangleMeshCurveSection *> internal_curve_section4_pt(1);

internal_curve_section4_pt[0] = boundary5_pt;

// The open curve that define this boundary is composed of just one

// curve section

inner_open_boundaries_pt[3] =

new TriangleMeshOpenCurve(internal_curve_section4_pt);

** Adding an internal boundary connected with another internal boundary (the last one from our example) **

The fifth internal boundary is connected to two internal boundaries. The initial end of the new internal boundary is connected to a `TriangleMeshPolyLine`

and the final end is connected to a `TriangleMeshCurviLine`

. The vertex number for connection with the `TriangleMeshPolyLine`

was established in the definition of the previous internal boundary, vertex #1. The circles defining boundaries 4 and 6 meet each other at which corresponds to the intrinsic coordinate on the `TriangleMeshCurviLine`

object representing boundary 4.

// Open curve 5

x_centre = 1.5;

y_centre = 0.0;

r_circle = 1.0;

Circle * boundary_circle3_pt = new Circle(x_centre, y_centre, r_circle);

// Number of segments used for representing the curve boundary

n_segments = 20;

// These numbers can be easily obtained by computing the

// intersection of the circle (circle3) and the other boundaries

// (in this case, boundaries 4 and 5)

s_start = atan2(0.5, (3.0 + sqrt(3.0)) * 0.5 - 1.5);

s_end = atan2(0.25*sqrt(7.0), -0.75);

// State the vertex number for connection on the destination

// boundaries

vertex_to_connect_initial = 1;

// State the s value for connection with the destination boundary

s_connection_final = atan2(0.25*sqrt(7.0), 0.75);

boundary_id = 6;

TriangleMeshCurviLine *boundary6_pt =

new TriangleMeshCurviLine(boundary_circle3_pt,

s_start,

s_end,

n_segments,

boundary_id);

// Do the connection with the destination boundary, in this case

// the connection is done with an internal boundary

boundary6_pt->connect_initial_vertex_to_polyline(

boundary5_pt,

vertex_to_connect_initial);

// Do the connection with the destination boundary, in this case

// the connection is done with an internal boundary

boundary6_pt->connect_final_vertex_to_curviline(

boundary4_pt,

s_connection_final);

Vector<TriangleMeshCurveSection *> internal_curve_section5_pt(1);

internal_curve_section5_pt[0] = boundary6_pt;

// The open curve that define this boundary is composed of just one

// curve section

inner_open_boundaries_pt[4] =

new TriangleMeshOpenCurve(internal_curve_section5_pt);

We build the mesh specific parameters using a `TriangleMeshParameters`

object.

//Create mesh

// Use the TriangleMeshParameters object for helping on the manage of the

// TriangleMesh parameters. The only parameter that needs to take is the

// outer boundary.

TriangleMeshParameters triangle_mesh_parameters(outer_boundary_pt);

// Specify the internal closed boundaries

triangle_mesh_parameters.internal_closed_curve_pt() = inner_boundaries_pt;

// Specify the internal open boundaries

triangle_mesh_parameters.internal_open_curves_pt() = inner_open_boundaries_pt;

// One can define an specific maximum element area

double element_area = 0.2;

triangle_mesh_parameters.element_area() = element_area;

// Pass the TriangleMeshParameters object to the TriangleMesh one

Problem::mesh_pt() = new TriangleMesh<ELEMENT>(triangle_mesh_parameters);

It is possible to define holes in the domain using the methods discussed in another tutorial; where we explained how to create holes by using closed curves. Holes could also be created by specifying additional parameters on the `TriangleMeshParameters`

object. The specification of regions or layers follows the same principle.

For defining holes it is only necessary to specify a point that lies inside a closed boundary (created by the connection of internal open boundaries). For example, the region bounded by boundaries 4, 5 and 6 can be turned into a hole by specifying the coordinate , see figure below.

// Adding a hole on the domain

Vector<Vector <double> > additional_hole(1);

// Define the coordinates on the domain

additional_hole[0].resize(2);

additional_hole[0][0] = 1.5;

additional_hole[0][1] = 0.75;

// Pass information about the additional holes coordinates

triangle_mesh_parameters.extra_holes_coordinates() = additional_hole;

// Pass the TriangleMeshParameters object to the TriangleMesh one

Problem::mesh_pt() = new TriangleMesh<ELEMENT>(triangle_mesh_parameters);

The definition of regions follows the same principle that defining holes, it means, one only need to specify the coordinates of a point that lies inside a region. If we would like to define the region showed on light grey on the figure below we use the method `add_region_coordinates`

of the `TriangleMeshParameters`

object. You can specify any region id for the defined region.

// Adding a region on the domain

Vector<double> region(2);

// Define the coordinates of the regions on the domain

region[0] = 0.5;

region[1] = 0.5;

// Pass information about the defined regions

triangle_mesh_parameters.add_region_coordinates(1, region);

// Specify different target area of region 1:

triangle_mesh_parameters.set_target_area_for_region(1,0.01);

// Pass the TriangleMeshParameters object to the TriangleMesh one

Problem::mesh_pt() = new TriangleMesh<ELEMENT>(triangle_mesh_parameters);

By default the complete domain belongs to region zero. Therefore, one could leave parts of the domain without a defined region or explicitly define all the regions on the domain. NOTE: Be sure to use the region ids when documenting the solution by its specified regions.

- Create two new internal boundaries without connections, one using a
`TriangleMeshPolyLine`

object and other using a`TriangleMeshCurviLine`

object.

- Define a new vertex on the new
`TriangleMeshPolyLine`

and create a connection to that vertex.

- Once one performs a connection to a
`TriangleMeshPolyLine`

or to a`TriangleMeshCurviLine`

there is an internal checking on the connection vertices coordinates for the`TriangleMeshPolyLine`

and on the connection intrinsic coordinate for the`TriangleMesheCurviline`

case, they must be close enough in order to allow the connection. The tolerance values are`tolerance_for_vertex_connection`

= for the`TriangleMeshPolyline`

case and`tolerance_for_zeta_connection`

= for the`TriangleMeshCurviline`

case. One can explicitly define the allowed tolerance by adding an extra argument when using the methods`connect_initial_vertex_to_polyline`

,`connect_final_vertex_to_polyline`

`connect_initial_vertex_to_curviline`

or`connect_final_vertex_to_curviline`

. As an exercise, use different values for the allowed connection tolerance.

- The source files for this tutorial are located in the directory:

demo_drivers/meshing/mesh_from_inline_triangle_internal_boundaries/

- The driver code is:

demo_drivers/meshing/mesh_from_inline_triangle/mesh_from_inline_triangle_internal_boundaries.cc

- The additional driver code

demo_drivers/meshing/mesh_from_inline_triangle_internal_boundaries/mesh_from_inline_triangle_internal_boundaries_extra.cc

shows how to define a hole by connecting internal boundaries.

A pdf version of this document is available.