... | ... | @@ -22,20 +22,20 @@ EikMesh is a 2D mesh generator for unstructured meshes. Its design goal is to pr |
|
|
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 high-quality unstructured and conforming triangular meshes. The user can introduce new data types.
|
|
|
Mesh elements work like nodes of a data collection, i.e. mesh can carry data which can be addressed in O(1) for a given mesh element.
|
|
|
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 document is to give the user a basic understanding of how to use the EikMesh library.
|
|
|
We will also describe implementation details, algorithms and concepts introduced in the field
|
|
|
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.
|
|
|
Two descriptions of the domain which we introduce in the following sections.
|
|
|
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 by analytical singed distance functions d such that,
|
|
|
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,
|
|
|
|
... | ... | @@ -51,7 +51,7 @@ Especially for domains containing sharp boundaries, it can be difficult to enfor |
|
|
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 approximating d by a prior-computed a background mesh.
|
|
|
This can be solved by an approximated distance function d_a which is prior-computed on background mesh.
|
|
|
Using this technique d is computed at points of the coarse background mesh and is interpolated elsewhere.
|
|
|
|
|
|
## Planar straight line graphs
|
... | ... | @@ -61,7 +61,8 @@ 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.
|
|
|
Note that we can split intersecting segments to enforce this criterion.
|
|
|
|
|
|
****Remark***: We can split intersecting segments to enforce the intersection criterion.
|
|
|
|
|
|
![pslg](uploads/7c16613f0ed611d0d4d0e27b8f6ea36f/pslg.png)
|
|
|
|
... | ... | @@ -85,22 +86,22 @@ segments not part of a polygon (the bounding polygon or any hole). |
|
|
## 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, filter or reducing a collection of elements (in parallel).
|
|
|
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 user-friendly notation.
|
|
|
|
|
|
## Implementation
|
|
|
|
|
|
The ``IMesh`` data structure contains all manages of the actual mesh, that is edges, faces, vertices and all data which is stored on those geometrical elements.
|
|
|
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 half-edge data structure also called doubly connected edge list (DCEL).
|
|
|
DCEL is able to manage planar straight-line graphs (PSLG).
|
|
|
A DCEL is able to manage planar straight-line 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 half-edges (one for each face of the edge). The figure above shows how mesh elements (vertices, edges / half-edges 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. For example to compute the degree of a vertex the following code suffice:
|
|
|
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:
|
|
|
|
|
|
```java
|
|
|
E e = getEdge(v);
|
... | ... | @@ -119,13 +120,13 @@ of the Delaunay triangulation or EikMesh or some other mesh generation method co |
|
|
|
|
|
### Pointer- and array-based implementation
|
|
|
We implemented two versions of the DCEL. For the pointer-based 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 \ie the object of type ``IHalfEdge``
|
|
|
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 object-oriented paradigm. Examples in this document we use these classes, but they can be exchange effortlessly.
|
|
|
This version is easier to debug and follows strongly the object-oriented paradigm. For examples in this document, we use these classes, but they can be exchange effortlessly.
|
|
|
|
|
|
The array-based version replaces pointers by indices.
|
|
|
Elements of the same type are contained in an array and can be identified by their index \ie their position in the array.
|
|
|
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
|
... | ... | @@ -146,7 +147,7 @@ IMesh<VPoint,Double,Double, |
|
|
PHalfEdge<VPoint,Double,Double>,
|
|
|
PFace<VPoint,Double,Double>> mesh = ...
|
|
|
```
|
|
|
However, we wanted to make it possible that the user can integrate his own ``XMesh`` implementation and at the same time, its should be easy to use meshing library easy.
|
|
|
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 local-variable 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.
|
... | ... | @@ -169,7 +170,7 @@ The assumptions are listed in the source code documentation (JavaDoc). One very |
|
|
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 half-edges and faces.
|
|
|
|
|
|
****Remark****: the order of points matter, i.e. they have to form a ****simple counter-clockwise (CCW) oriented polygon****.
|
|
|
****Remark****: The order of points matter, i.e. they have to form a ****simple counter-clockwise (CCW) oriented polygon****.
|
|
|
|
|
|
Some operations offered by ``IMesh`` require the construction of new points ``P``, therefore it requires a ``IPointConstructor``, here defined by the lambda expression ``(x,y) -> new VPoint(x,y)``:
|
|
|
|
... | ... | @@ -390,9 +391,7 @@ var triangulation = dT.generate(); |
|
|
### 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](https://hal.inria.fr/inria-00167201/document). The following images visualizes the triangulation before and after the removal.
|
|
|
|
|
|
![tri_removal](uploads/da17611f4915a1dc88a28842b69d6531/tri_removal.png)
|
|
|
We implemented the algorithm presented by Devillers in [Devillers, 2002](https://hal.inria.fr/inria-00167201/document).
|
|
|
|
|
|
```java
|
|
|
var mesh = triangulation.getMesh();
|
... | ... | @@ -401,16 +400,20 @@ var deletePoints = mesh.getPoints(face); |
|
|
triangulation.remove(deletePoints.get(0));
|
|
|
```
|
|
|
|
|
|
The following images visualizes the triangulation before and after the removal.
|
|
|
|
|
|
![tri_removal](uploads/da17611f4915a1dc88a28842b69d6531/tri_removal.png)
|
|
|
|
|
|
## The constrained Delaunay triangulation
|
|
|
|
|
|
Another well-known 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.
|
|
|
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](https://doi.org/10.1016/0045-7949(93)90239-A).
|
|
|
|
|
|
## The confirming Delaunay triangulation
|
|
|
## 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.
|
... | ... | @@ -436,7 +439,7 @@ The following illustrated the difference between the DT (left), the CDT (middle) |
|
|
|
|
|
# Mesh generation
|
|
|
|
|
|
The meshing algorithms in this section expect as input at least an edge length function h and a segment-bounded planar straight line graph (PSLG) G.
|
|
|
Each meshing algorithm in this section expects at least an edge length function h and a segment-bounded planar straight line graph (PSLG) G as input.
|
|
|
|
|
|
## The edge length function
|
|
|
|
... | ... | @@ -450,7 +453,7 @@ is the edge length function defined on the mesh domain D. In the following examp |
|
|
|
|
|
## Rebay's algorithms
|
|
|
|
|
|
The eikmesh library offers the implementation of two algorithm introduced by Rebay in [Rebay, 1993](https://doi.org/10.1006/jcph.1993.1097). The first one, the so called *Voronoi-Vertex point insertion method*,
|
|
|
The eikmesh library offers the implementation of two algorithms introduced by Rebay in [Rebay, 1993](https://doi.org/10.1006/jcph.1993.1097). The first one, the so called *Voronoi-Vertex point insertion method*,
|
|
|
is a *Delaunay-refinement* technique, i.e. additional so called *Steiner vertices* are inserted at the circumcenter of some Delaunay triangle.
|
|
|
The second one, called *Voronoi-Segment point insertion method*, is a *Frontal-Delaunay algorithm*.
|
|
|
Frontal-Delaunay algorithms are a hybridization of advancing-front and Delaunay-refinement 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 *advancing-front techniques*.
|
... | ... | @@ -484,7 +487,7 @@ The gray code indicate the triangle quality. The result for h0 = 0.3, 0.1 and 0. |
|
|
|
|
|
### Voronoi-Segment point insertion
|
|
|
|
|
|
For the first 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 segment-bound of G. If so, we ignore the vertex, i.e. it will not be inserted
|
|
|
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 segment-bound of G. If so, we ignore the vertex, i.e. it will not be inserted
|
|
|
|
|
|
For more information we refer to [Rebay, 1993](https://doi.org/10.1006/jcph.1993.1097).
|
|
|
|
... | ... | @@ -546,14 +549,15 @@ Ruppert's algorithm is able to mesh complex geometries defined by PSLGs containi |
|
|
|
|
|
# EikMesh <a name="EikMEsh"></a>
|
|
|
|
|
|
EikMesh uses our implementation of the [DCEL](https://gitlab.lrz.de/vadere/vadere/wikis/eikmesh/The-Mesh-Data-Structure) as underlying data structure. EikMesh requires three basically three ingredients:
|
|
|
1. a bounding box `VPolygon` (containing all mesh points)
|
|
|
2. a signed distance function d(x) which gives the distance to the topography border
|
|
|
3. an edge length function h(x) which gives the desired edge length depending on the midpoint of the edge
|
|
|
4. an initial edge `h0` length which determines the minimal edge length i.e. h(x) <= h0 everywhere.
|
|
|
EikMesh uses our implementation of the [DCEL](https://gitlab.lrz.de/vadere/vadere/wikis/eikmesh/The-Mesh-Data-Structure) as underlying data structure. EikMesh requires basically three ingredients:
|
|
|
1. a bounding box ``VPolygon`` (containing all mesh points)
|
|
|
2. a signed distance function ``d`` which gives the distance to the topography border
|
|
|
3. an edge length function ``h`` which gives the desired edge length depending on the midpoint of the edge
|
|
|
4. 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.
|
|
|
|
|
|
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.
|
|
|
For more details we refer to [Zönnchen and Köster, 2018](https://www.sciencedirect.com/science/article/pii/S1877750318303193).
|
|
|
The EikMesh algorithm is heavily based on [DistMesh](http://persson.berkeley.edu/distmesh/), a simple and effective mesh generator (in MATLAB) which was developed by Per-Olof 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 straight-line graphs (PSLG) and uses a more accessible and efficient data structure. For more details we refer to [Zönnchen and Köster, 2018](https://www.sciencedirect.com/science/article/pii/S1877750318303193).
|
|
|
|
|
|
![headMesh](uploads/76ced2ddd4ef0e71d3f76043bcbcad83/headMesh.png)
|
|
|
|
... | ... | @@ -602,13 +606,14 @@ The second line defines a hole and is transformed into a distance function ``d_r |
|
|
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, -d_rec)
|
|
|
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 non-constant edge length function
|
|
|
****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 non-constant edge length function
|
|
|
|
|
|
h(x) = h_{min} + 0.3 * |d_comb(x)|.
|
|
|
|
|
|
#### Result:
|
|
|
|
|
|
The result for h0 = 0.1, 0.03 and 0.01.
|
|
|
|
|
|
![simple_combine](uploads/9156ad93e1147040f96032d76433d146/simple_combine.png)
|
... | ... | @@ -654,8 +659,6 @@ 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.
|
|
|
|
|
|
****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.
|
|
|
|
|
|
#### Result:
|
|
|
|
|
|
![eikmesh_d_combined_png](uploads/4f86de4de6bec3e4aa77870b63216685/eikmesh_d_combined_png.png)
|
... | ... | @@ -724,6 +727,8 @@ meshImprover.generate(); |
|
|
|
|
|
![kaiserslautern](uploads/af46213c0e9e02f254c740e71e0709d8/kaiserslautern.png)
|
|
|
|
|
|
****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](http://www.vadere.org/the-eikmesh-library/). |
|
|
\ No newline at end of file |