MESH2D Automatic 2D Mesh Generation

MESH2D
is a MATLAB program which generates
unstructured meshes in 2D,
by Darren Engwirda.

The code is relatively simple, flexible and powerful. The user is able
to define a variety of geometric shapes, and desired mesh densities.

MESH2D is most useful because it allows a user to specify
a shape or region, which the program will then fill with a triangular mesh.
The density of the triangular mesh can be uniform, or the user can request
that smaller triangles be used near certain features of the region. The
program relies heavily on the features of the Delaunay triangulation,
which chooses, among all possible triangulations of a set of points, that
triangulation which best avoids small angles.

Interested users should refer to the copy of MESH2D
that is made available through the MATLAB Central File Exchange.
This copy is essentially my personal working copy, to which I may
have added comments, small coding changes, and extra tests and examples.

The MESH2D routines include a call to a function called wbar().
The call to this function fails on my system, and since it
seems to have no importance whatsoever, I commented it out.
(I believe it is intended to generate a "wait bar", similar
to the hourglass or beachball or wristwatch icons.

The MESH2D function "mytsearch()" was originally written to call
MATLAB's "tsearch()" function. The tsearch() function has since been
removed from MATLAB. Therefore, the call to tsearch() also
causes MESH2D to fail!

One alternative is a file called
tsearch_mex.c, which searches a triangulation to determine which triangle
contains each point. It does not require that the triangulation be
Delaunay. To use this function with MATLAB, you need to apply MATLAB's
MEX compiler...if you have never used the MEX compiler before, you may
have some difficulty, since you need to determine that you have the MEX
compiler, that you have a C or C++ compiler on your system, that MEX
knows where these compilers are, and that you know how to invoke MEX
to compile the function. That should be something like

mex tsearch_mex.c

(The second time you do something like this is, of course, a hundred
times easier and only half as mysterious!)

A second alternative is to replace the call to tsearch() by a
call to MATLAB's replacement function DelaunayTri; however, a simple
substitution of one call for the other does not quite work. There is,
apparently, some feature of tsearch() that is not available in
DelaunayTri(). In particular, it may be that tsearch() did not require
the triangulation to be Delaunay...

A third alternative is to replace the call to tsearch(x,y,t,px,py) by a
call to tsearchn([x y], t, [px py] ), which seems to work.

Usage:

[ p, t ] = mesh2d ( vertices, edge, hdata, options );

where:

vertices, required input, a V by 2 list of (X,Y) coordinates of
vertices of the boundary. If vertices is the only input argument,
then it must be the case that the vertices are listed consecutively.
Otherwise, assuming edge is supplied, the vertices can be given
in any order.

edge, optional input, a V by 2 list of pairs of indices in the
vertices array that constitute the edges of the polygonal boundary.
If vertices is actually already in order, then edge, if
specified, would contain the values [1,2; 2,3; 3,4; ... ; V,1].

options, optional input that allows the user to modify the default
behavior of the solver (see below).

p, (output), an N by 2 list of node coordinates. The number of
nodes generated, N, is determined in part by the size of the edges
along the boundary, and by other user input such as the maximum element
size, and the user size function, if supplied.

t, (output), an M by 3 list of node indices, forming counterclockwise
triangles. The number of triangles created depends on the number of nodes
created.

hdata, the element size information. This structure, if supplied,
can include the following information:

hdata.hmax, the maximum allowable global element size.

hdata.edgeh, an array of element sizes on specified geometry
edges, where e1 is an index into the edge array. The
edgeh component would have the form [e1,h1; e2,h2; ...], where
the user has specified a certain number of sizes.

hdata.fun, the name of a function preceded by an AT sign,
which is the user-defined size function. fun must have the form

h = fun ( x, y, args{} )

where x and y are vectors of point coordinates, and args
is an optional addition set of input set in hdata.args.
The function returns the user-desired elementsize at the given points.

options allows the user to modify the default behavior of the solver.
This structure, if supplied, can include the following information:

options.mlim is the convergence tolerance. The maximum relative
change in edge length per iteration must be less than this value, which
defaults to 0.02.

options.maxit, the maximum allowable number of iterations,
which defaults to 20.

options.dhmax, the maximum allowable (relative) gradient in the size
function, which defaults to 0.3.

options.output, a "logical" variable which displays the mesh
and the mesh statistics upon completion, and defaults to "TRUE", that is, 1.

Licensing:

Copyright (c) 2009, Darren Engwirda
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

tinterp.m
carries out linear interpolation at points within a triangle, given
function values at the vertices.

triarea.m
computes the area of one or more triangles, assuming their
vertices are given in counterclockwise order.

tsearch_mex.c,
which carries out a search of the triangulation; Matlab's "tsearch()" function
already does this, but is a) officially obsolete, resulting in lots of warnings;
b) requires that the triangulation be Delaunay.
In "mytsearch()", replace the call to "tsearch()" by a call to "tsearch_mex()".
To use this file, it must be compiled within MATLAB (once) using a command like

mex tsearch_mex.c

Examples and Tests:

facedemo.m
demonstrates two example polygonal geometries for input to MESHFACES.

REFINE_DEMO demonstrates how the refinement feature of MESH2D
can be controlled with a mask, which specifies that only certain
triangles should be refined, rather than all of them. The masked
triangles are split into 4 subtriangles. Triangles immediately neighboring
such triangles will be split in half. Remaining triangles will be untouched.

refine_demo2.png
the mesh after refining all the triangles whose centroid has a Y
coordinate greater than 1/2.

refine_demo3.png
the mesh after a second refinement, on all the triangles whose
centroid has an X coordinate less than 1/2.

SQUARE_BORDER demonstrates how a constrained mesh can be created
using MESHFACES. In this case, we have a square in a square. We want
to mesh the two resulting regions in such a way that they share the
nodes along the common edges.