## DistMesh2D.jl

Simple implementation of DistMesh algorithm allowing to generate meshes on 2D plane.
Author jstarczewski
Popularity
7 Stars
Updated Last
7 Months Ago
Started In
December 2020

# DistMesh2D.jl  The package provides Julia implementation of DistMesh algorithm developed by Per-Olof Persson and Gilbert Strang allowing to generate meshes on 2D plane, using DelDir.jl which is a Julia wrapper for Delaunay triangulations and Voronoi/Dirichlet tessellations. Before using this package I highly recommend reading this document covering basic MATLAB's use-cases, because this packages tries to provide similar runtime interface to original one.

## Installation

To use DistMesh2D.jl clone the repository and add it to the local package registry.

## Usage

To generate a mesh, define a signed distance function, desired edge length function, bounding box, distance between points in initial distribution and run the meshing algorithm.

```julia> fdistance(p) = sqrt(sum(p .^ 2)) - 1
fdistance (generic function with 1 method)

julia> fedgelength(p) = 1.0
fedgelength (generic function with 1 method)

julia> boundingbox = [-1.0 -1.0; 1.0 1.0]
2×2 Array{Float64,2}:
-1.0  -1.0
1.0   1.0

julia> initdistance = 0.2
0.2

julia> x, y = distmesh2d(fdistance, fedgelength, boundingbox, initdistance)
([-0.8408297336717067, -0.6993829085766673, NaN, -0.5315455685570857, -0.6993829085766673, NaN, -0.5822722080639808, -0.6993829085766673, NaN, -0.5822722080639808  …  NaN, -0.60618605865758, -0.36756281453508466, NaN, -0.60618605865758, -0.7846917737263828, NaN, -0.60618605865758, -0.5633752938295506, NaN], [-0.5412996931425073, -0.7147471921795109, NaN, -0.8470296986801119, -0.7147471921795109, NaN, -0.5521369576265936, -0.7147471921795109, NaN, -0.5521369576265936  …  NaN, 0.795322867083831, 0.92999869760509, NaN, 0.795322867083831, 0.6198861346071187, NaN, 0.795322867083831, 0.5800262182849156, NaN])```

The output of a given function are points that are ready to be plotted with one of the available plotting libraries.

```julia> using Gadfly

julia> Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))```

## Generating more complex meshes

Library provides simple generic signed distance functions to define popular shapes like `dcircle` and `drectangle`.

```function drectangle(
p,
x1::T,
x2::T,
y1::T,
y2::T,
)::T where {T<:Float64}
return -min(min(min(-y1 + p, y2 - p, -x1 + p, x2 - p)))
end

function dcircle(
p,
xc::T,
yc::T,
r::T
)::T where {T<:Float64}
return sqrt((p - xc) .^ 2 + (p - yc) .^ 2) - r
end```

Library provides simple operations on shapes like union, difference and intersections with via functions like `dunion` `ddiff` and `dintersect`. In each definition `p` is considered point. Other parameters are proper signed distance functions.

```function dunion(
p,
dfs::Function...,
)::Float64
return min([fun(p) for fun in dfs]...)
end

function ddiff(
p,
d1::Function,
dfs::Function...
)::Float64
return max(d1(p), [-fun(p) for fun in dfs]...)
end

function dintersect(
p,
dfs::Function...
)::Float64
return max([fun(p) for fun in dfs]...)
end```

Other useful functions are `huniform` which provides uniform height distribution for point `p` and `protate` which rotates given point by passed value.

```function huniform(p)::Float64
return ones(size(p, 1), 1)
end

function protate(p,phi)
return p*[cos(phi) -sin(phi); sin(phi) cos(phi)]
end```

### Non-uniform height function

Meshes with non uniform height functions also can be easily plotted. First we define signed distance functions which in this example is a rectangle with circle hole inside.

```julia> outrect(p) = drectangle(p, -100.0, 100.0, -100.0, 100.0)
outrect (generic function with 1 method)

julia> incircle(p) = dcircle(p, 0.0, 0.0, 40.0)
incircle (generic function with 1 method)

julia> fd(p) = ddiff(p, outrect, incircle)
fd (generic function with 1 method)```

Then we define height functions which is not going to return the same value for every input parameters, but will give finer resolution closer to the circle.

```julia> fh(p) = min(4 * sqrt(sum(p.^2)) - 100)
fh (generic function with 1 method)```

Next we define a bounding box `bbox` that can fit our rectangle and a value representing distance between points in initial distribution.

```julia> bbox = [-100.0 -100.0; 100.0 100.0]
2×2 Array{Float64,2}:
-100.0  -100.0
100.0   100.0

julia> h0 = 10.0
10.0```

To let the rectangle withstand the triangulation we need to defined its corner points as fixed.

```julia> pfix = [-100.0 -100.0; -100.0 100.0; 100.0 -100.0; 100.0 100.0]
4×2 Array{Float64,2}:
-100.0  -100.0
-100.0   100.0
100.0  -100.0
100.0   100.0```

Finally we define the `dptol` parameter which represents the limit value. Algorithms will loop as long as all movements, relative to bar lengths in given iteration are not below it.

```julia> dptol = 0.003
0.003

julia> x, y = distmesh2d(fd, fh, bbox, h0, pfix=pfix, dptol=dptol)
([-100.0, -76.71718215838415, NaN, -62.29844095886764, -76.71718215838415, NaN, -100.0, -76.71718215838415, NaN, -100.0  …  NaN, 100.0, 100.0, NaN, 100.0, 78.01048490854643, NaN, 100.0, 80.72224653575613, NaN], [-80.56252474196684, -77.31440723125387, NaN, -62.263859671319544, -77.31440723125387, NaN, -65.15571042094268, -77.31440723125387, NaN, -65.15571042094268  …  NaN, 79.78462487451092, 63.972558404391755, NaN, 79.78462487451092, 76.7946778947107, NaN, 79.78462487451092, 100.0, NaN])

julia> Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))```

### Composing shapes

Let's say we want to mesh a simplified version of R2D2 head. By "simplified" I mean a semicircle with a hole in the center. We are going to work on `bbox = [0.0 0.0; 20.0 10.0]` and `h0 = 0.25`. First we define the hole as a small circle.

```julia> sc(p) = dcircle(p, 10.0, 5.0, 3.0)
sc (generic function with 1 method)```

Then we define a bigger circle.

```julia> bc(p) = dcircle(p, 10.0, 0.0, 10.0)
bc (generic function with 1 method)```

Next we want the smaller circle to stay empty during triangulation so we calculate the difference between functions.

```julia> fds(p) = ddiff(p, bc, sc)
fds (generic function with 1 method)```

Finally we define a rectangle and calculate its intersection with the area represented by the `fds` function.

```julia> dr(p) = drectangle(p, 0.0, 20.0, 0.0, 10.0)
dr (generic function with 1 method)

julia> fd(p) = dintersect(p, dr, fds)
fd (generic function with 1 method)```

We also need two fixed points.

```julia> pfix = [0.0 0.0; 20.0 0.0]
2×2 Array{Float64,2}:
0.0  0.0
20.0  0.0```

Finally we mesh our distance, passing additional `logmi` parameter with value equal to `true`. This parameter slows down the algorithm, but turns on logging of the minimal move index value that occured in all triangulations so far. Monitoring minimal move index is helpful to determine whether the algorithm is going forward and closer to `dptol` which by default is equal to `0.001`.

```julia> x, y = distmesh2d(fd, huniform, bbox, h0, pfix=pfix, logmi=true)
[ Info: Minimal move index is 1.7848378957838162
[ Info: Minimal move index is 0.6341210762677251
[ Info: Minimal move index is 0.3189998493116862
.
.
.
.
[ Info: Minimal move index is 0.0010124268021825809
[ Info: Minimal move index is 0.0010102509818988507
[ Info: Minimal move index is 0.0010080823842259872
[ Info: Minimal move index is 0.0010059209713983509
[ Info: Minimal move index is 0.0010037667055315528
[ Info: Minimal move index is 0.0010016195522329304
([0.39168705623677563, 0.24949716765305802, NaN, 1.4066677906807563, 1.1484883086226838, NaN, 2.089037409389501, 1.9403068274377, NaN, 1.6763442046107901  …  NaN, 5.225956577772305, 5.45051039916194, NaN, 5.225956577772305, 4.9344362730766385, NaN, 5.225956577772305, 5.523495989156209, NaN], [2.063690666551136, 2.2198411283347967, NaN, 1.619386786583647, 1.6285058990860535, NaN, 0.0, 0.23721296745446183, NaN, 0.23008747967549248  …  NaN, 8.786837276630676, 8.730320688506946, NaN, 8.786837276630676, 8.622068433075366, NaN, 8.786837276630676, 8.942086553129595, NaN])```

## Error handling

Considering the fact that this package relies on DelDir.jl to triangulate points, some errors occurring in triangulation process can be hard to solve, because they can occur in original triangulation package itself. Those errors are caught, wrapped and exposed to the user with information about iteration, triangulated dataset and inner error itself. Supposing we changed `bbox` and `h0` from previous example to:

```julia> bbox = [0.0 0.0; 10.0 10.0]
2×2 Array{Float64,2}:
0.0   0.0
10.0  10.0

julia> h0 = 2.5
2.5```

Which means that the current area defined by `bbox` is too small to triangulate points. We are going to face the following error.

```julia> x, y = distmesh2d(fd, huniform, bbox, h0, pfix=pfix)
ERROR: Triangulation failed with given data:
Points = [0.125 0.21650635094610965; 0.375 0.21650635094610965; 0.25 0.4330127018922193; 0.375 0.649519052838329; 0.625 0.21650635094610965; 0.5 0.4330127018922193; 0.625 0.649519052838329; 0.875 0.21650635094610965; 0.75 0.8660254037844386; 1.125 0.21650635094610965; 1.0 0.8660254037844386]
Iteration = 0
Inner error = DomainError([0.0, 1.0, 0.0, 1.0], "Boundary window is too small")
Stacktrace:
 distmesh2d(::typeof(fd), ::typeof(huniform), ::Array{Float64,2}, ::Float64; pfix::Array{Float64,2}, dptol::Float64, ttol::Float64, geps::Float64, Fscale::Float64, deltat::Float64, deps::Float64, logmi::Bool) at /Users/jstarczewski/GitProjects/DistMesh2D.jl/src/distmesh.jl:33
 top-level scope at none:1
julia> ```

Note that when filling an issue it is highly recommended to check whether the error comes from meshing algorithm or directly from triangulation.