|
|
1. [Introduction](#Introduction)
|
|
|
2. [The mesh data structure](#The-mesh-data-structure)
|
|
|
3. [Preliminary algorithms](#Preliminary-algorithms)
|
|
|
|
|
|
# Introduction <a name="Introduction"></a>
|
|
|
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 high-quality meshes for complex spatial domains.
|
|
|
|
|
|
The EikMesh library focuses on unstructured meshes which are highly flexible and require less elements, compared to structured meshes, to mesh the same domain.
|
|
|
Since unstructured meshes can be highly irregular, they require more memory per element (vertex, face, edge) because the can not be saved implicitly.
|
|
|
Furthermore, addressing and changing the mesh require more complicated algorithms.
|
|
|
|
|
|
EikMesh is a 2D mesh generator for unstructured meshes. Its design goal is to provide a fast, light and user-friendly 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 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.
|
|
|
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
|
|
|
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.
|
|
|
|
|
|
## Signed distance functions
|
|
|
|
|
|
One way is to define the domain 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 approximating d by a prior-computed a background mesh.
|
|
|
Using this technique d is computed at points of the coarse background mesh and is interpolated elsewhere.
|
|
|
|
|
|
## Planar straight line graphs
|
|
|
Another well-known 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.
|
|
|
Note that we can split intersecting segments to enforce this criterion.
|
|
|
|
|
|
![pslg](uploads/7c16613f0ed611d0d4d0e27b8f6ea36f/pslg.png)
|
|
|
|
|
|
For mesh generation, a PSLG must be segment-bounded.
|
|
|
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 segment-bounded too.
|
|
|
PSLGs are especially useful for describing non-curved 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 <a name="The-mesh-data-structure"></a>
|
|
|
|
|
|
<img align="center" src="uploads/6ab4c46afd78c0ee7761112992394f01/halfEdgeDS.png" alt="Square" width="400"/>
|
|
|
|
|
|
## 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).
|
|
|
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.
|
|
|
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).
|
|
|
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:
|
|
|
|
|
|
```java
|
|
|
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 half-edge 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 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``
|
|
|
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.
|
|
|
|
|
|
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.
|
|
|
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.
|
|
|
|
|
|
### Generic data types
|
|
|
We use a generic implementation such that the user can define types which can be associated with, and accessed via, a specific mesh element.
|
|
|
Therefore the ``PMesh<P,CE,CF>`` defines three generic types:
|
|
|
|
|
|
- ``P`` is the type of objects which can be associated with a ``PVertex``. ``P`` has to be a point (``P extends IPoint``)
|
|
|
- ``CE`` is the type of objects which can be associated with a ``PHalfEdge`` and
|
|
|
- ``CF`` is the type of objects which can be associated with a ``PFace`.
|
|
|
|
|
|
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:
|
|
|
```java
|
|
|
IMesh<VPoint,Double,Double,
|
|
|
PVertex<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.
|
|
|
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.
|
|
|
|
|
|
```java
|
|
|
PMesh<VPoint,Double,Double> 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 ****counter-clockwise (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 half-edges and faces.
|
|
|
|
|
|
****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)``:
|
|
|
|
|
|
```java
|
|
|
var mesh = new PMesh<>((x,y) -> new VPoint(x,y));
|
|
|
mesh.toFace(new VPoint(0,0),new VPoint(1,0),
|
|
|
new VPoint(1,1),new VPoint(0,1));
|
|
|
```
|
|
|
|
|
|
To use the array-based implementation instead, we only have to exchange ``PMesh`` by ``AMesh``:
|
|
|
|
|
|
```java
|
|
|
var mesh = new AMesh<>((x,y) -> new VPoint(x,y));
|
|
|
mesh.toFace(new VPoint(0,0),new VPoint(1,0),
|
|
|
new VPoint(1,1),new VPoint(0,1));
|
|
|
```
|
|
|
|
|
|
#### Result
|
|
|
|
|
|
![square](uploads/69e6f81c6c2bce0bad4962ef905518a5/square.png)
|
|
|
|
|
|
### Read a PSLG file
|
|
|
|
|
|
The EikMesh libarary is able to transform a feasible segment-bounded PSLG file
|
|
|
into a ``IMesh`` or internal geometry object like ``VPolygon`` and ``VLine``.
|
|
|
The following code reads the ``A.poly`` file from an ``InputStream``:
|
|
|
|
|
|
```java
|
|
|
InputStream inputStream = ...
|
|
|
PSLG pslg = PolyGenerator.toVShapes(inputStream);
|
|
|
```
|
|
|
|
|
|
To transform a PSLG file into a ``IMesh`` the following code suffice:
|
|
|
|
|
|
```java
|
|
|
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, half-edge or face), to access other adjacent mesh elements, the ``IMesh`` object has to be used.
|
|
|
Read the following code from right to left.
|
|
|
|
|
|
```java
|
|
|
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 half-edge et of this edge. And finally we access the face ft of the twin half-edge 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 half-edges, vertices / points, neighbouring faces of a specific face
|
|
|
- all half-edges 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 half-edge which has only one neighboring face:
|
|
|
|
|
|
```java
|
|
|
for(VPoint point : mesh.getPointIt(mesh.getBorder())) {
|
|
|
...
|
|
|
}
|
|
|
```
|
|
|
|
|
|
The following gives a parallel stream of the same elements:
|
|
|
|
|
|
```java
|
|
|
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, half-edge, face) can carry some data of the type defined by the use.
|
|
|
In the following we save on each face its area.
|
|
|
|
|
|
```java
|
|
|
for(PFace<VPoint, Double, Double> face : mesh.getFaces()) {
|
|
|
mesh.setData(face, mesh.toTriangle(face).getArea());
|
|
|
}
|
|
|
```
|
|
|
|
|
|
Then we can, for example, compute the area which is triangulated:
|
|
|
|
|
|
```java
|
|
|
double areaSum = dt.getMesh()
|
|
|
.streamFaces()
|
|
|
.mapToDouble(f -> dt.getMesh().getData(f).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
|
|
|
|
|
|
```java
|
|
|
Triangulated area = 97.10080319694295
|
|
|
Average triangle area = 0.04939003214493538
|
|
|
Area triangulated = 97.10080319694295 %
|
|
|
```
|
|
|
|
|
|
We can also exchange ``VPoint`` by any point container extending ``IPoint``.
|
|
|
|
|
|
### 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.
|
|
|
|
|
|
![visualization](uploads/75fefa30fa768154bb787030fd8ab46c/visualization.png)
|
|
|
|
|
|
#### 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:
|
|
|
|
|
|
```java
|
|
|
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:
|
|
|
|
|
|
```java
|
|
|
var panel = new PMeshPanel<>(
|
|
|
dt.getMesh(),
|
|
|
500,
|
|
|
500,
|
|
|
f -> dt.getMesh().toTriangle(f).isNonAcute() ? red : Color.WHITE);
|
|
|
panel.display("Delaunay triangulation");
|
|
|
```
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Preliminary algorithms <a name="Preliminary-algorithms"></a>
|
|
|
|
|
|
## 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-2001], are implemented:
|
|
|
|
|
|
- straight walk (default)
|
|
|
- orthogonal walk
|
|
|
- probabilistic walk
|
|
|
|
|
|
![straightWalk](uploads/2144580cc84481e075ed586c0203e1ab/straightWalk.png)
|
|
|
|
|
|
Furthermore we implemented different point localization algorithms and data structures:
|
|
|
|
|
|
- Jump and Walk (default) [devroye-2004]
|
|
|
- Delaunay-Tree [guibas-1992]
|
|
|
- Delaunay-Hierarchy [devillers-2002]
|
|
|
- plain walk (no additional strategy),
|
|
|
|
|
|
For a random insertion order of n points, using the Delaunay-Tree or the Delaunay-Hierarchy 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 Delaunay-Tree or the Delaunay-Hierarchy 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:
|
|
|
|
|
|
```java
|
|
|
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:
|
|
|
|
|
|
```java
|
|
|
// (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<VPoint, Double, Double>(points,
|
|
|
(x,y) -> new VPoint(x,y));
|
|
|
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. The following images visualizes the triangulation before and after the removal.
|
|
|
|
|
|
![tri_removal](uploads/da17611f4915a1dc88a28842b69d6531/tri_removal.png)
|
|
|
|
|
|
```java
|
|
|
var mesh = triangulation.getMesh();
|
|
|
var face = triangulation.locateFace(new VPoint(5,5)).get();
|
|
|
var deletePoints = mesh.getPoints(face);
|
|
|
triangulation.remove(deletePoints.get(0));
|
|
|
```
|
|
|
|
|
|
## 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.
|
|
|
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 confirming 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.
|
|
|
|
|
|
```java
|
|
|
final InputStream inputStream = ...
|
|
|
boolean conforming = false;
|
|
|
PSLG pslg = PolyGenerator.toPSLGtoVShapes(inputStream);
|
|
|
var cdt = new PContrainedDelaunayTriangulator<VPoint, Double, Double>(
|
|
|
pslg,
|
|
|
(x,y) -> new VPoint(x,y),
|
|
|
conforming);
|
|
|
var triangulation = cdt.generate();
|
|
|
```
|
|
|
|
|
|
The following illustrated the difference between the DT (left), the CDT (middle) and the CCDT (right) of a PSLG.
|
|
|
|
|
|
![CDT](uploads/28bdb462fa809bee789aaca6e5566c4e/CDT.png) |
|
|
\ No newline at end of file |