distmesh2d: A simple mesh generator for non-convex regions

In geometry: Mesh Generation and Surface Tessellation

A simple mesh generator for non-convex regions

Description

An unstructured simplex requires a choice of mesh points (vertex nodes) and a triangulation. This is a simple and short algorithm that improves the quality of a mesh by relocating the mesh points according to a relaxation scheme of forces in a truss structure. The topology of the truss is reset using Delaunay triangulation. A (sufficiently smooth) user supplied signed distance function (fd) indicates if a given node is inside or outside the region. Points outside the region are projected back to the boundary.

Usage

distmesh2d(
  fd,
  fh,
  h0,
  bbox,
  p = NULL,
  pfix = array(0, dim = c(0, 2)),
  ...,
  dptol = 0.001,
  ttol = 0.1,
  Fscale = 1.2,
  deltat = 0.2,
  geps = 0.001 * h0,
  deps = sqrt(.Machine$double.eps) * h0,
  maxiter = 1000,
  plot = TRUE
)

Arguments

fd	Vectorized signed distance function, for example mesh.dcircle or mesh.diff, accepting an n-by-2 matrix, where n is arbitrary, as the first argument.
fh	Vectorized function, for example mesh.hunif, that returns desired edge length as a function of position. Accepts an n-by-2 matrix, where n is arbitrary, as its first argument.
h0	Initial distance between mesh nodes. (Ignored of p is supplied)
bbox	Bounding box cbind(c(xmin,xmax), c(ymin,ymax))
p	An n-by-2 matrix. The rows of p represent locations of starting mesh nodes.
pfix	nfix-by-2 matrix with fixed node positions.
...	parameters to be passed to fd and/or fh
dptol	Algorithm stops when all node movements are smaller than dptol
ttol	Controls how far the points can move (relatively) before a retriangulation with delaunayn.
Fscale	“Internal pressure” in the edges.
deltat	Size of the time step in Euler's method.
geps	Tolerance in the geometry evaluations.
deps	Stepsize \Delta x in numerical derivative computation for distance function.
maxiter	Maximum iterations.
plot	logical. If TRUE (default), the mesh is plotted as it is generated.

Details

This is an implementation of original Matlab software of Per-Olof Persson.

Excerpt (modified) from the reference below:

‘The algorithm is based on a mechanical analogy between a triangular mesh and a 2D truss structure. In the physical model, the edges of the Delaunay triangles of a set of points correspond to bars of a truss. Each bar has a force-displacement relationship f(\ell, \ell_{0}) depending on its current length \ell and its unextended length \ell_{0}.’

‘External forces on the structure come at the boundaries, on which external forces have normal orientations. These external forces are just large enough to prevent nodes from moving outside the boundary. The position of the nodes are the unknowns, and are found by solving for a static force equilibrium. The hope is that (when fh = function(p) return(rep(1,nrow(p)))), the lengths of all the bars at equilibrium will be nearly equal, giving a well-shaped triangular mesh.’

See the references below for all details. Also, see the comments in the source file.

Value

n-by-2 matrix with node positions.

Wishlist

• Implement in C/Fortran

• Implement an nD version as provided in the Matlab package

• Translate other functions of the Matlab package

Author(s)

Raoul Grasman

References

http://persson.berkeley.edu/distmesh/

P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004

See Also

tri.mesh, delaunayn, mesh.dcircle, mesh.drectangle, mesh.diff, mesh.union, mesh.intersect

Examples

# examples distmesh2d
fd <- function(p, ...) sqrt((p^2)%*%c(1,1)) - 1
     # also predefined as `mesh.dcircle'
fh <- function(p,...)  rep(1,nrow(p))
bbox <- matrix(c(-1,1,-1,1),2,2)
p <- distmesh2d(fd,fh,0.2,bbox, maxiter=100)
    # this may take a while:
    # press Esc to get result of current iteration

# example with non-convex region
fd <- function(p, ...) mesh.diff(p , mesh.drectangle, mesh.dcircle, radius=.3)
     # fd defines difference of square and circle

p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3,maxiter=4)
p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3, maxiter=10)
     # continue on previous mesh

geometry documentation built on Sept. 11, 2024, 9 p.m.