Introduction
Computational modeling and simulation is a critical aspect of modern science and industrial design, underpinning a diverse range of applications, from numerical simulations of physical phenomena, such as fluid dynamics to crowd simulation to computer graphics and animation. The basis of these simulations is a discrete representation of the underlying geometry, i.e. spatial domain. Modern mesh generation is concerned with the development of efficient and automatic algorithms to construct and maintain highquality meshes for complex spatial domains.
The EikMesh library focuses on unstructured meshes which are highly flexible and therefore, require less elements, compared to structured meshes. Since unstructured meshes can be highly irregular, they require more memory per element (vertex, face, edge) because the connectivity of elements can not be saved implicitly. Regardless of the higher memory consumption per element, the overall memory is often significantly lower. Addressing and changing the mesh requires more complicated algorithms but, due to the reduced number of mesh elements, the overall computation time is often lower.
EikMesh is a 2D mesh generator for unstructured meshes. Its design goal is to provide a fast, light and userfriendly meshing tool
with parametric input and advanced visualization capabilities. This library generates exact Delaunay triangulations (DT),
constrained Delaunay triangulations (CDT), conforming Delaunay triangulations (CCDT), Voronoi diagrams,
and highquality unstructured and conforming triangular meshes. The user can introduce new data types.
Mesh elements work like nodes of a data collection, i.e. mesh elements can carry data which can be addressed in O(1) time for a given mesh element.
Therefore, a mesh is also an abstract data type like a List
.
The goal of this page is to give the user a basic understanding of how to use the EikMesh library. We will also scratch the surface of algorithms and concepts introduced in the field of computational geometry which are not necessarily required for using the library. Feel free to skip those parts.
Mesh generation begins with the domain to be meshed. In 2D this domain is a subset of the Euclidean space. The EikMesh library uses two descriptions of the domain which we introduce in the following section.
Signed distance functions
One way is to define the domain is by analytical singed distance functions d such that, d gives the distance to the domain boundary and is negative inside and positive outside the domain. For example,
d(x) = x  1,
defines a disc. Especially for curved geometries where the gradient of d is unique, geometries defined by a distance function is a powerful technique. For more information we refer to the work of Persson and Strang.
Using distance functions can lead to some disadvantages. Especially for domains containing sharp boundaries, it can be difficult to enforce geometrical conformity at positions at which the gradient is not unique. One solution is to add additional fix points at those positions. Another problem might be the computational expensive evaluation of more complex distance functions. This can be solved by an approximated distance function d_a which is priorcomputed on background mesh. Using this technique d is computed at points of the coarse background mesh and is interpolated elsewhere.
Planar straight line graphs
Another wellknown domain representation are planar straight line graphs (PSLG)s. A PSLG consist of vertices and segments. The generated mesh has to contain all segments and points of the PSLG. We distinguish between segments of the PSLG and edges of the mesh. A segment can be represented by multiple edges of the generated mesh. The PSLG is planar, i.e. segments only intersect at a shared vertex.
*Remark: We can split intersecting segments to enforce the intersection criterion.
For mesh generation, a PSLG must be segmentbounded. In 2D this means that there exist a set of segments in the PSLG which form a simple polygon which contains all points and therefore all other segments. This special polygon separates the mesh domain from the exterior domain. Furthermore, the mesh domain can contain holes which have to be segmentbounded too. PSLGs are especially useful for describing noncurved domains but they can lead to a high number of small segments if a curve has to be approximated.
Remark: Note that we can transform any simple polygon into a distance function thus we can transform the PSLG description into a distance function description but we lose all segments not part of a polygon (the bounding polygon or any hole).
The mesh data structure
Streams, monads and lambdas
Java 8 introduced some concept known from functional programming to the Java programming language. For example the concept of lazy evaluation and monads such as the Optional
monad.
Additionally, the Stream
library is a powerful extension to write shorter and more readable code especially for iterating, filtering or reducing a collection of elements (in parallel).
To use a Stream
one also has to get used to the new lambda expressions, like p > p.x + p.y
, which is syntactical sugar for writing anonymous classes.
The EikMesh library excessively uses these language features to establish a userfriendly notation.
Implementation
The IMesh
data structure contains and manages all elements of the actual mesh, that is edges, faces, vertices and all data which is stored on those geometrical elements.
Therefore, to create, delete, change any element the user has to call methods offered by IMesh
.
For example, to create a new vertex we can call mesh.createVertex(x,y)
.
IMesh
is an implementation of the edge based halfedge data structure also called doubly connected edge list (DCEL).
A DCEL is able to manage planar straightline graphs (PSLG).
Note that a triangulation is a special case of a PSLG. Each edge of the mesh is represented by two
counter clockwise (CCW) oriented halfedges (one for each face of the edge). The figure above shows how mesh elements (vertices, edges / halfedges and faces)
are connected and how one can access neighbouring elements in constant time.
We decided to use the DCEL since it is very flexible and all local operations are rather fast (in theory O(1)). For example to compute the degree of a vertex the following code suffice:
E e = getEdge(v);
E n = edge;
int degree = 0;
do {
n = getTwin(getNext(n))
degree++;
} while(n != e)
Creating a mesh by hand requires a lot of code since all relations, e.g. a halfedge and its next, have to be set. However, this is not necessary to use our algorithms like the computation of the Delaunay triangulation or EikMesh or some other mesh generation method contained in the EikMesh library.
Pointer and arraybased implementation
We implemented two versions of the DCEL. For the pointerbased implementation each relation between two objects in IMesh
such as a IHalfEdge
and the IVertex
in which the IHalfEdge
ends, is realized via Java references, i.e. the object of type IHalfEdge
possesses a reference pointing to the object of type IVertex
representing its end point.
Classes in this category start with a P, e.g. PMesh
, PHalfEdge
, PFace
, PVertex
.
This version is easier to debug and follows strongly the objectoriented paradigm. For examples in this document, we use these classes, but they can be exchange effortlessly.
The arraybased version replaces pointers by indices.
Elements of the same type are contained in an array and can be identified by their index, i.e. their position in the array.
Classes in this category start with an A, e.g. AMesh
, AHalfEdge
, AFace
, AVertex
.
The advantage of this approach is that those arrays can be sorted with respect to some spatial criterion.
Consequently, elements geometrically close can be sorted in such a way that they are close in memory which
improves the performance of algorithms iterating over geometrically local elements.
Properties
Each mesh element (vertex, halfedge, face) can have different userdefined properties. Each property has to have a unique String name
. Properties are managed by the mesh, i.e. defining, inserting and accessing properties can be done via the mesh object by the following methods:

setData(V vertex, String name, CV value)
, whereCV
is a generic type 
setData(E edge, String name, CE value)
, whereCE
is a generic type 
setData(F edge, String name, CF value)
, whereCF
is a generic type 
CV
getData(V vertex, String name, Class), where
CV` is a generic type 
CE
getData(E edge, String name, Class), where
CE` is a generic type 
CF
getData(F face, String name, Class), where
CF` is a generic type
Remark: The user is responsible for the correctness of the types that is inserting the objects / primitives of the same type for one key / name and using the correct Class<>
. The following code, for example, would cause an Exception
:
mesh.setData(vertex, "velocity", 3.0);
mesh.getData(vertex, "velocity", Boolean.class);
In Java there is nothing like typedef
thus the programmer has always to fully specify the data type which might lead to long lines of code like:
IMesh<PVertex, PHalfEdge, PFace> mesh = ...
However, we wanted to make it possible that the user can integrate his own XMesh
implementation and at the same time, it should be easy to use the meshing library.
Therefore, we introduce some classes which predefine some types like PMesh
. We also encourage the user to use the localvariable type inference introduced in Java 10
by replacing long type definitions by the ``var` keyword if possible.
The code above and below are semantically the same.
PMesh mesh = ...
Implicit assumptions
Many operations make assumptions about the mesh it is operating on. For example the operation splitting a triangle into 3 triangles assumes that the split point lies inside the triangle and that the triangle is in fact a valid triangle. We insert assertions to test many of those assumptions but these assertions should only be active for testing because the performance overhead can be huge. The assumptions are listed in the source code documentation (JavaDoc). One very important one is that the mesh is always counterclockwise (CCW) oriented.
Working with meshes
Building a mesh by hand
The following code example creates a mesh consisting of one square defined by the point set containing (0,0), (1,0), (1,1) and (0,1).
Furthermore, we save Double
objects on its halfedges and faces.
Remark: The order of points matter, i.e. they have to form a simple counterclockwise (CCW) oriented polygon.
var mesh = new PMesh();
mesh.toFace(new VPoint(0,0),new VPoint(1,0),new VPoint(1,1),new VPoint(0,1));
To use the arraybased implementation instead, we only have to exchange PMesh
by AMesh
:
var mesh = new AMesh();
mesh.toFace(new VPoint(0,0),new VPoint(1,0),new VPoint(1,1),new VPoint(0,1));
Result
Read a PSLG file
The EikMesh libarary is able to transform a feasible segmentbounded PSLG file
into a IMesh
or internal geometry object like VPolygon
and VLine
.
The following code reads the A.poly
file from an InputStream
:
InputStream inputStream = ...
PSLG pslg = PolyGenerator.toVShapes(inputStream);
To transform a PSLG file into a IMesh
the following code suffice:
IMesh<...> mesh = ...
InputStream inputStream = ...
PolyGenerator.toPMesh(inputStream, mesh);
Here we assume that the mesh
is empty. All segments of the PSLG will be inserted.
Accessing mesh elements
Direct Access
Given any mesh element (vertex, halfedge or face), to access other adjacent mesh elements, the IMesh
object has to be used.
Read the following code from right to left.
mesh.getFace(mesh.getTwin(mesh.getEdge(mesh.getFace())));
First we get access to some (arbitrary) face f of the mesh. Via the second call we access some (arbitrary) edge e of the face f. Then we get the twin halfedge et of this edge. And finally we access the face ft of the twin halfedge et. Therefore, we access a neighboring face of f.
Iterators and streams
To simplify certain access pattern such as accessing all edges of a specific face the library offers different Iterator
s and Stream
s to go over
 all halfedges, vertices / points, neighbouring faces of a specific face
 all halfedges ending at a vertex / point
 all faces surrounding a specific vertex / point
 all faces, edges, vertices of a mesh
 ...
The following code iterates over all points of the mesh which are at the border, i.e. those points are connected to at least one halfedge which has only one neighboring face:
for(VPoint point : mesh.getPointIt(mesh.getBorder())) {
...
}
The following gives a parallel stream of the same elements:
mesh.streamPoints(mesh.getBorder()).parallel();
Remark: Note that those iterators and streams rely on the connectivity of the mesh itself. Therefore, changing the connectivity while iterating will cause unpredictable behaviors and will possibly destroy the validity of the mesh!
Containers
As already mentioned each mesh element (vertex, halfedge, face) can carry some data of the type defined by the use. In the following we save on each face its area.
for(PFace face : mesh.getFaces()) {
mesh.setData(face, "area", mesh.toTriangle(face).getArea());
}
Then we can, for example, compute the area which is triangulated:
double areaSum = dt.getMesh()
.streamFaces()
.mapToDouble(f > dt.getMesh().getData(f, "area", Double.class).get())
.sum();
double averageArea = areaSum / dt.getMesh().getNumberOfFaces();
double triangulatedArea = (100 * (areaSum / (width * height)));
System.out.println("Triangulated area = " + areaSum);
System.out.println("Average triangle area = " + averageArea);
System.out.println("Area triangulated = " + triangulatedArea + " %");
Result
Triangulated area = 97.10080319694295
Average triangle area = 0.04939003214493538
Area triangulated = 97.10080319694295 %
Visualization
During the development of the EikMesh library we created some useful utility classes to accelerate the development process. Visualizing the mesh helped us to understand and improve certain algorithms and to write documents like this.
Tikz generator
The first utility enables the user to convert any valid IMesh
into a Tikz file to generate high quality vector graphic figures.
The following code generates Tikz code where each non acute triangle of a Delaunay triangulation is painted in red:
TexGraphGenerator.toTikz(
dt.getMesh(),
f > dt.getMesh().toTriangle(f).isNonAcute() ? red : Color.WHITE,
1.0f)
Live visualization
Furthermore, we can display the same mesh on a Java Swing
canvas:
var panel = new PMeshPanel(
dt.getMesh(),
500,
500,
f > dt.getMesh().toTriangle(f).isNonAcute() ? red : Color.WHITE);
panel.display("Delaunay triangulation");
Preliminary algorithms
Point location
Given a triangular unstructured 2D mesh and a point p, a very important operation in computational geometry is to find the face f which contains it. Every following algorithm uses this basic operation multiple times. A fast implementation is key for fast unstructured mesh generation. To solve the point location problem different walking strategies, described in Devillers et al., 2001, are implemented:
 straight walk (default)
 orthogonal walk
 probabilistic walk
Furthermore we implemented different point localization algorithms and data structures:
 Jump and Walk (default) Devroye et al., 2004
 DelaunayTree Guibas et al., 1992
 DelaunayHierarchy Devillers, 2002
 plain walk (no additional strategy),
For a random insertion order of n points, using the DelaunayTree or the DelaunayHierarchy leads to a time complexity of
O(log(n)) for each point location and using the Jump and Walk strategy requires O(n^(1/4)) time.
However, surprisingly the Jump and Walk algorithm does not require any additional data structure thus offers the most flexibility.
Additionally, it performs very well in practice, often better than its two alternatives.
Basically a ITriangulation
is the combination of a triangular mesh IMesh
and a point location algorithm IPointLocator
.
Remark: For bookkeeping reasons, the EikMesh library does not support point removal if the DelaunayTree or the DelaunayHierarchy is used.
The first call of the following code example uses the fact that our mesh is a triangulation, i.e. dt
is of type ITriangulation
.
The call starts the so called Jump and Walk algorithm which is fast.
The second one checks each face in a brute force manner until it finds the face containing the point, which is slow:
dt.locateFace(5,5);
dt.getMesh().locate(5,5);
The Delaunay triangulation
The Delaunay triangulation (DT) is one of the most important structures in computational geometry and the bases for a whole set of algorithms that generate unstructured meshes. The so called Delaunay criterion states that a valid triangulation is Delaunay if and only if there is no point contained in any circumcircle of any triangle. We will refer to DT(V) as the Delaunay triangulation of the vertex set V.
DistMesh computes the Delaunay triangulation multiple times to construct the connectifity of the mesh while EikMesh avoids this repetitive computation. To be flexible, i.e. to insert and remove points after the initial triangulation has finished we implemented the incremental method of Lawson. The following code constructs a Delaunay triangulation of 100 random points uniformly distributed inside a square:
// (1) generate a point set
Random random = new Random(0);
int width = 10;
int height = 10;
int numberOfPoints = 100;
var supply = () > new VPoint(random.nextDouble()*width,random.nextDouble()*height);
Stream<VPoint> randomPoints = Stream.generate(supply);
List<VPoint> points = randomPoints
.limit(numberOfPoints)
.collect(Collectors.toList());
// (2) compute the Delaunay triangulation
var dT = new PDelaunayTriangulator(points);
var triangulation = dT.generate();
Result
Point removal
The next code snippet removes the first point p of the face located at q = (5.0, 5.0). We implemented the algorithm presented by Devillers in Devillers, 2002.
var mesh = triangulation.getMesh();
var face = triangulation.locateFace(new VPoint(5,5)).get();
var deletePoints = mesh.getPoints(face);
triangulation.remove(deletePoints.get(0));
The following images visualizes the triangulation before and after the removal.
The constrained Delaunay triangulation
Another wellknown triangulation is the so called constrained Delaunay triangulation which is often used instead of the Delauny trianuglaiton if certain segments, for example, defined by a planar straight line graph (PSLG), have to be part of the triangulation. We call CDT(V) a constrained Delaunay triangulation of V.
The algorithm we implemented was presented by Sloan: In the first step we compute the Delaunay triangulation for the PSLG G. In the second step we restore edges of G by edge flips. For more details we refer to Sloan, 1993.
The conforming Delaunay triangulation
The conforming Delaunay triangulation CCDT is constrained as well as a Delaunay triangulation. That is, each constrained segments is the union of edges in the CCDT and the triangulation is in fact a Delaunay triangulation as well.
The following code generates the CDT of the PSLG. If we replace conforming = false
by conforming = true
the CCDT will be computed.
final InputStream inputStream = ...
boolean conforming = false;
PSLG pslg = PolyGenerator.toPSLGtoVShapes(inputStream);
var cdt = new PContrainedDelaunayTriangulator(
pslg,
conforming);
var triangulation = cdt.generate();
The following illustrated the difference between the DT (left), the CDT (middle) and the CCDT (right) of a PSLG.
Mesh generation
Each meshing algorithm in this section expects at least an edge length function h and a segmentbounded planar straight line graph (PSLG) G as input.
The edge length function
The goal of unstructured mesh generation is to produce high quality meshes containing as less elements as possible. To control the size of the elements, the user has to define an edge length function. More precisely,
h : D > R^+
is the edge length function defined on the mesh domain D. In the following examples you will find different edge length functions.
Rebay's algorithms
The eikmesh library offers the implementation of two algorithms introduced by Rebay in Rebay, 1993. The first one, the so called VoronoiVertex point insertion method, is a Delaunayrefinement technique, i.e. additional so called Steiner vertices are inserted at the circumcenter of some Delaunay triangle. The second one, called VoronoiSegment point insertion method, is a FrontalDelaunay algorithm. FrontalDelaunay algorithms are a hybridization of advancingfront and Delaunayrefinement techniques, in which a Delaunay triangulation is used to define the topology of a mesh while new Steiner vertices are inserted in a manner consistent with advancingfront techniques. Both algorithms have a time complexity of O(n log(n)), where n is the number of vertices of the generated triangulation. For both the user can control the element size by an edge length function h. Both algorithm expect a segmentbounded PSLG G as input and terminates if h is satisfied everywhere.
In our implementation we start by constructing the constrained Delaunay triangulation of the PSLG G. In the second step we insert additional Steiner vertices based on the rules described by Rebay.
VoronoiVertex point insertion method
For the first algorithm, in each iteration k, the Steiner vertex v is the circumcenter of the triangle t_j in the current triangulation with the largest circumcenter radius. Note that if the triangulation is a Delaunay triangulation, v is in fact a vertex of the Voronoi diagram of the current vertex set. It might be the case that v lies outside of the segmentbound of G. If so, we ignore the vertex, i.e. it will not be inserted.
The following code generates a mesh using h(x) = 0.05:
PSLG pslg = ...
double h0 = 0.05;
var vviMethod = new PVoronoiVertexInsertion(
pslg,
p > h0);
var triangulation vviMethod.generate();
Result:
The gray code indicate the triangle quality. The result for h0 = 0.3, 0.1 and 0.05.
VoronoiSegment point insertion
For the second algorithm, in each iteration k, the Steiner vertex v lies on the segment of the Voronoi diagram of the current triangulation. Additional each inserted vertex is an attempt to generate a new triangle t_j such that h is satisfied, which is not always possible. Like before, c might be outside of the segmentbound of G. If so, we ignore the vertex, i.e. it will not be inserted
For more information we refer to Rebay, 1993.
The following code generates a mesh using h(x) = 0.05:
PSLG pslg = ...
double h0 = 0.05;
var vviMethod = new PVoronoiSegmentInsertion(
pslg,
p > h0);
var triangulation vviMethod.generate();
Result:
The gray code indicate the triangle quality. The result for h0 = 0.3, 0.1 and 0.05.
Ruppert's algorithm
Ruppert's algorithm Ruppert, 1993 for 2D quality mesh generation is probably the first theoretically guaranteed meshing algorithm to be truly satisfactory in practice. The algorithm allows the density of triangles to vary quickly over short distance. It is quite similar to Rebay's VoronoiVertex point insertion method, but avoids Steiner vertices outside the domain D. Furthermore, it is also a Delaunayrefinement method. Additionally to the edge length function h, the user can define a minimal allowed angle theta such that the final mesh will not contain any angle smaller than theta. For an excellent and extensive description we refer to shewchuk, 2002.
Remark: The termination of Ruppert's algorithm is only guaranteed for a minimal angle greater than theta = 20.7 degree and no angle smaller 60 degree in the given PSLG G. The extension in shewchuk, 2002 resolves the issue of small input angles in G but at this moment EikMesh only supports the algorithm proposed by Ruppert. A Cimplementation is accessible via Triangle.
The following code generates a mesh using h(x) = 0.05 and a minimal allowed theta angle of 20 degree:
PSLG pslg = ...
double h0 = 0.02;
double theta = 20.0;
var ruppert = new PRuppertsTriangulator(
pslg,
p > h0,
theta);
var triangulation = ruppert.generate();
Result:
The gray code indicate the triangle quality. The result for theta = 20 degree and h(x) = infinity, theta = 30 degree and h(x) = infinity, theta = 20 and h(x) = 0.02.
Ruppert's algorithm is able to mesh complex geometries defined by PSLGs containing sharp corners such as the greenland.poly
:
EikMesh
EikMesh uses our implementation of the DCEL as underlying data structure. EikMesh requires basically three ingredients:
 a bounding box
VPolygon
(containing all mesh points)  a signed distance function
d
which gives the distance to the topography border  an edge length function
h
which gives the desired edge length depending on the midpoint of the edge  an initial edge
h0
length which determines the minimal edge length i.e. h(x) <= h0 everywhere.
However if (2) is missing, EikMesh constructs d(x) by using only the bounding box (and maybe additional VShape
s). If (3) is missing EikMesh assumes h(x) = h0 everywhere leading to a mesh of balanced edge lengths.
The EikMesh algorithm is heavily based on DistMesh, a simple and effective mesh generator (in MATLAB) which was developed by PerOlof Persson and Gilbert Strang. EikMesh inherits from DistMesh the specification of the geometry via signed distance functions and the concept of iterative smoothing by converging towards a force equilibrium. However, EikMesh completely avoids the computation of the Delaunay triangulation, generates a different and cache friendly initial triangulation and treats boundary elements more carefully. Additionally, EikMesh supports geometries defined by a segment bounded planar straightline graphs (PSLG) and uses a more accessible and efficient data structure. For more details we refer to Zönnchen and Köster, 2018.
Examples
Example 1 (Uniform ring)
The following code snippet generates a mesh using the standard EikMesh algorithm for a simple curved geometry. The geometry is defined by the distance function d(x) = xc0.50.4, where c = (1,1) to be the center of the ring. The inner radius is 0.2 and the outer 1.0. The edge length function is h(x) = h0 = 0.1.
VRectangle bound = new VRectangle(0.1, 0.1, 2.2, 2.2);
IDistanceFunction d_r = IDistanceFunction.createRing(1, 1, 0.2, 1.0);
double h0 = 0.1;
PEikMesh meshImprover = new PEikMesh(d_r,h0,bound);
meshImprover.generate();
Result:
The result for h0 = 0.3, 0.1 and 0.05.
Example 2 (Subtract rectangle from disc)
Distance functions can be combined.
VRectangle bound = ...
VRectangle rect = new VRectangle(0.5, 0.5, 1, 1);
IDistanceFunction d_comb = IDistanceFunction.createDisc(0.5, 0.5, 0.5);
IDistanceFunction d_rec = IDistanceFunction.create(rect);
IDistanceFunction d = IDistanceFunction.substract(d_comb, d_rec);
double edgeLength = 0.03;
var meshImprover = new PEikMesh(
d,
p > edgeLength + 0.5 * Math.abs(d.apply(p)),
edgeLength,
bound,
Arrays.asList(rect),
(x, y) > new EikMeshPoint(x, y, false));
Line 1 defines a bounding box containing the mesh domain.
The second line defines a hole and is transformed into a distance function d_rec
in line 4.
The disc is defined by the distance function d_disc constructed in line 3.
By subtracting d_rect from d_disc we construct
d(x) = d_comb(x) = max(d_disc(x), d_rec(x))
Remark: Even if the geometry is defined by d_comb
, the VRectangle
rect
is an argument of EikMesh. This is optional but it will automatically inserts fix points at the corners of the rectangle which improves the mesh quality. Furthermore, in line 9 we use a nonconstant edge length function
h(x) = h_{min} + 0.3 * d_comb(x).
Result:
The result for h0 = 0.1, 0.03 and 0.01.
Example 3 (Combining distance functions)
We can also construct geometries by the union and intersection of distance functions. For example, the figure below illustrate the result of the following code.
// inner rectangle
VRectangle rect = new VRectangle(0.5, 0.5, 1, 1);
// outer rectangle
VRectangle boundary = new VRectangle(2,0.7,4,1.4);
// construction of the distance function, define the 2 discs
IDistanceFunction d1_c = IDistanceFunction.createDisc(0.5, 0, 0.5);
IDistanceFunction d2_c = IDistanceFunction.createDisc(0.5, 0, 0.5);
// define the two rectangles
IDistanceFunction d_r = IDistanceFunction.create(rect);
IDistanceFunction d_b = IDistanceFunction.create(boundary);
// combine distance functions
IDistanceFunction d_unionTmp = IDistanceFunction.union(d1_c, d_r)
IDistanceFunction d_union = IDistanceFunction.union(d_unionTmp, d2_c);
IDistanceFunction d = IDistanceFunction.substract(d_b,d_union);
// h_min
double edgeLength = 0.03;
var meshImprover = new PEikMesh(
d,
p > edgeLength + 0.5 * Math.abs(d.apply(p)),
edgeLength,
GeometryUtils.boundRelative(boundary.getPath()),
Arrays.asList(rect));
// generate the mesh
var triangulation = meshImprover.generate();
Here we define two circles and a rectangle and compute the union d_union
of these three distance functions in line 7. Finally we subtract d_union
from a larger rectangle.
Result:
Example 4 (Meshing a PSLG)
The following code sniped produces the mesh illustrated in below by constructing a distance functions based on the given PSLG:
PSLG pslg = ...
PEikMesh meshImprover = new PEikMesh(pslg.getSegmentBound(), 0.02, pslg.getHoles());
Result:
Example 5 (Smoothing a given mesh)
The following code uses a Delaunay triangulation of 1000 random points as the initial triangulation, also depicted below (left).
var dt = ... // the Delaunay triangulation
var eikMesh = new PEikMesh(p > 1.0 + Math.abs(bound.distance(p)),dt.getTriangulation());
var triangulation = eikMesh.generate();
The quality of the resulting triangulation is approximately 0.96 (using the same measure as Persson). This example shows that EikMesh can be used to improve even a very low quality initial triangulation.
Result:
Mesh before (left) and after (right) of the smoothing stage.
Example 6 (Urban environment)
In the last example we mesh an urban environment (600m x 250m) which is in fact a geometry used in Vadere.
// (1) read the PSLG from a file / an input stream
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/kaiserslautern.poly");
PSLG pslg = PolyGenerator.toPSLGtoVShapes(inputStream);
// (2) construct a distance function d out of the PSLG
Collection<VPolygon> holes = pslg.getHoles();
VPolygon segmentBound = pslg.getSegmentBound();
IDistanceFunction d = IDistanceFunction.create(segmentBound, holes);
// (3) use EikMesh to construct the mesh
double h0 = 5.0;
var meshImprover = new PEikMesh(
d,
p > h0 + 0.3 * Math.abs(d.apply(p)),
h0,
new VRectangle(segmentBound.getBounds2D()),
pslg.getHoles()
);
meshImprover.generate();
Remark: Since the distance function d is complicated the mesh generation is quite slow because d has to be evaluated for each iteration of the smoothing process for almost any mesh element. To reduce the complexity of the evaluation one can construct a background mesh using, for example, Ruppert's algorithm and approximate d via the interpolation on the background mesh.
Result:
Remark: The definition of the distance function d, the edge length function h, the minimum edge length h0 and the bounding box have to be carefully defined by the user. For example, if h0 is too large EikMesh might fail to construct a valid mesh.
Videos
To get a good understanding of how EikMesh works, we constructed some videos of the meshing process which can be found at vadere.org.