# Creating a 2D meshing algorithm in Mathematica

As what is proving to be a difficult, but entertaining task, I am attempting to adapt a 2D meshing algorithm created for MATLAB and port it to Mathematica. I understand meshing functions already exist in Mathematica so this is purely for fun/learning.

The paper here describes in great detail the algorithm created for MATLAB, in terms of the non code specific logic and each step of the code. I intend to modify this question as/if suggestions are made to show its progress.

To start off simple, I have a box that I'm creating a random selection of meshing points within. I modify these using Delaunay's Triangulation Algorithm and then create an initial triangular mesh. Next, imagining all these points are joined my 'magic springs', if the length, L, of a spring is less than a certain chosen scale length, L0, then it exerts a force on each joining point in the direction of increasing separation, proportional to L-L0. Thus you create a new set of points at a later artificial timestep due to the effect of these psuedoforces. Then if any points are outside of the original box, you drag 'em back in and repeat this process until some sort of equilibrium condition is met.

This is the code I've done so far:

SeedRandom[2];
bboxL = 10;
bboxH = 10;
Nnodes = 100;
p0 = Table[0, {i, Nnodes}, {i, 2}];
Do[
p0[[i, 1]] = RandomReal[{0, bboxL}];
p0[[i, 2]] = RandomReal[{0, bboxH}];
, {i, Nnodes}];

timeStep = 0.2;

Needs["ComputationalGeometry"];
dt = DelaunayTriangulation[p0];
toPairs[{m_, ns_List}] := Map[{m, #} &, ns];
edges = Flatten[Map[toPairs, dt], 1];
Graphics[GraphicsComplex[
p0, {Line[edges], Red, PointSize[Medium], Point[p0]}]]

L0 = 1;
displacement[{m_, ns_List}] :=
Total[Map[
If[Norm[p0[[#, All]] - p0[[m, All]]] <
bboxL*bboxH/Nnodes, -p0[[#, All]] + p0[[m, All]], {0, 0}] &,
ns]];

Do[
p0 = p0 + timeStep*Map[displacement, dt];
Do[
p0[[i, 1]] = Max[p0[[i, 1]], 0];
p0[[i, 2]] = Max[p0[[i, 2]], 0];
p0[[i, 1]] = Min[p0[[i, 1]], bboxL];
p0[[i, 2]] = Min[p0[[i, 2]], bboxH];
, {i, Nnodes}];
dt = DelaunayTriangulation[p0];
edges = Flatten[Map[toPairs, dt], 1];

, {i, 400}];

Graphics[GraphicsComplex[
p0, {Line[edges], Red, PointSize[Medium], Point[p0]}]]


You start off with this mesh:

Then after some iterations you get:

What I'd like to do next is instead of meshing to a boring old square I'd like to be able to mesh to any geometry whether its given in a functional form or an array of points. Can anybody suggest how to proceed? Maybe just a simple circle would be a nice next step for practice?

Everyone feel free to edit my post and join in the fun. Many thanks for any help :)

Also many thanks go to MarkMcClure for helping me on getting to this stage!

-
@gpap true, but I might want the bounding box to be rectangular instead of square so would want different ranges of random numbers for the the two columns. If there's a way to implement that in one line that'd be great. – Tom Mar 28 '13 at 16:24
assume a polygon boundary, include the bounding polygon vertices in the Delaunay traingulation and lock them during the iteration. You'll of course need to do an insideness test to keep your points in the polygon. – george2079 Mar 28 '13 at 18:46
@ThomasJebbSturges, I have ported distmesh to M- a while back, unfortunately, I do not have the time to clean it up and post it here. If you are interested though, I can send it to you via email. You can reach me at ruebenko AT wolfram com. – user21 Mar 28 '13 at 20:34
That is funny. I coded this Persson's algorithm in Mathematica myself some years ago. I remember it worked but it wasn't fast (at least as I coded it). I will try and find the code and share some example meshes. – BoLe Apr 11 '13 at 12:17
You say meshing functions already exist in Mathematica. Yeah? I don't know any. – BoLe Apr 11 '13 at 12:19

first part..i had lying around..

poly = Random[Real, {1, 2}] {Cos[#], Sin[#]} & /@  Sort[Table[Random[Real, {0, 2 Pi}], {5}]]
isLeft[P2_, {P0_, P1_}] := -Sign@Det@{P2 - P0, P1 - P0};
pinpoly[p_, poly_] := Module[{ed},(*winding rule*)
ed = Partition[Append[poly, poly[[1]]], {2}, 1];
Count[ed,pr_ /; (pr[[1, 2]] <= p[[2]] < pr[[2, 2]] &&isLeft[p, pr] == 1)]
- Count[ed, pr_ /; (pr[[2, 2]] <= p[[2]] < pr[[1, 2]] && isLeft[p, pr] == -1)] == 1];
seed = Select[ Table[RandomReal[{-2, 2}, 2], {2000}] ,  pinpoly[#, poly] &];
Show[{Graphics[Line[Append[poly, poly[[1]]]]],         Graphics[{Red, Point /@ seed}]}]
allp = Join[seed, poly];


Show[{Graphics[(Line /@ (allp[[#]] & /@ Thread[#])) & /@  DelaunayTriangulation[allp]],
Graphics[{Red, PointSize[.02], Point /@ seed}],
Graphics[{Blue, PointSize[.02], Point /@ poly}]}]


EDIT--step 2: relaxing seed node positions, note I'm not doing anythng to keep in bounds..so you see a few go out where the boundary is not convex.

lastmf = 0;
dmf = 1;
Monitor[While[dmf > 10^-3,
allp = Join[seed, poly];
tri = DelaunayTriangulation[allp][[;; Length[seed]]];
alledges = Union@(Sort /@ Flatten[ Thread /@ tri , 1]);
(*unit vector associated with each edge*)
alledgesunitv = ((allp[[#[[2]]]] - allp[[#[[1]]]]) // Normal) & /@  alledges;
(* edgemap is a list of positions on elledges at each seed node,
signed to indicate if the unit vector points toward or away from the node *)
edgemap =  Map[Function[t, (If[Length[(p = Position[alledges, #])] == 1,
p[[1, 1]], -Position[alledges, Reverse[#]][[1, 1]]]) & /@ Thread@ t], tri];
elen = Norm[(allp[[#[[1]]]] - allp[[#[[2]]]])] & /@ alledges;
targetlen = .8 Mean[elen];
(*note the mean does not include the bounding polygon  edge lengths*)
eforce =  1 - targetlen/elen; (*nonlinear spring.. need to play with this..*)
{mf, dmf, lastmf} = {#, Abs[# - lastmf], mf} &@ Max@Abs@eforce;
seed += (Total /@
(Map[Sign[#] .01  eforce[[Abs[#]]]  alledgesunitv[[Abs[#]]] & ,   edgemap , {2}]));
(* add code here to keep nodes in bounds..*)
g = Show[{Graphics[{Thick, Blue, Line[Append[poly, poly[[1]]]]}],
Graphics[Line /@ Map[allp[[#]] &, alledges, {2}]],
Graphics[{Red, PointSize[.02], Point /@ seed}],
Graphics[{Blue, PointSize[.02], Point /@ poly}]}]], {g, mf, dmf}]


I expect at the boundary you want to not only keep nodes inside, but look for nodes close inside and snap those to the edge.

Edit 3 ..

force the bounary polygon to have a minimum edge length:

   seedlen = .25
poly = Flatten[(kk = Ceiling[Norm[Subtract @@ # ]/(1.2 seedlen )] ;
Table[#[[1]] (  1 + ( 1 - ici )/kk)  + (ici - 1)/ kk #[[2]] , {ici, kk}]) &
/@  Partition[Append[poly, poly[[1]]], {2}, 1], 1]


and simply discard nodes that wander out of bounds , where i had "add code here" above..

   seed = Select[seed, pinpoly[#, poly] &];


works pretty good.. note when you drop nodes you need to render the plot after recomputing the triangulation.

-
Very nice! Thanks for posting this, will have a good look through this over the next few days. – Tom Mar 29 '13 at 19:23

Hm, I am late for the party but anyway here is my entry. This distmesh port works in 2D and 3D (though this has issues I should say) and does not need external code like qhull. It also has a quality control / max steps termination and boundary points to be included can be given. Note, however, this is a prototype code and it has issues. A word about distmesh before the code: distmesh is probably the simplest mesh generator (well if you have a Delaunay code) and produces high quality meshes. The downside, unfortunately, is that the algorithm is slow due to it's repeated re-delaunay-zation.

To understand the code you most likely will need some familiarity with the Person's Distmesh paper.

Preliminary Code:

Some computational geometry stuff:

ddiff[d1_, d2_] := Max[d1, -d2];
dunion[d1_, d2_] := Min[d1, d2];
dintersection[d1_, d2_] := Max[d1, d2];
dcircle[{x_, y_}, {cx_, cy_, r_}] := (x - cx)^2 + (y - cy)^2 - r^2
drectangle[{x_, y_}, {{x1_, y1_}, {x2_, y2_}}] := -Min[Min[Min[-y1 + y, y2 - y], -x1 + x], x2 - x]
dline[{x_,y_}, {{x1_, y1_}, {x2_,y2_}}] := ((x2 - x1)*(y2 - y) - (y2 - y1)*(x2 - x)) / Sqrt[(x2 - x1)^2 + (y2 - y1)^2 ]
dtriangle[{x_, y_}, {a : {x1_, y1_}, b : {x2_, y2_}, c : {x3_, y3_}}] :=
Max[ dline[{x, y}, {a, b}], dline[{x, y}, {b, c}],dline[{x, y}, {c, a}]]
dpoly[{x_, y_}, a : {{_, _} ..}] :=
Max @@ (dline[{x, y}, #] & /@ Partition[ a, 2, 1, {1, 1}])


Distmesh helper functions:

(* this is not quite correct, the coords in the mesh may not be the \
same as the ones origianly given as pts; e.g. when pts contain \
duplicate points *)
myDelaunay[ pts_, dim_] := Block[{tmp},
Switch[ dim,
2,
tmp = GraphicsRegionDelaunayMesh[ pts];
DeveloperToPackedArray[ tmp["MeshObject"]["MeshElements"]],
3,
TetGenDelaunay[ pts][[2]]
]
];

thresholdPosition[ p_, th_ ] :=
Flatten[ SparseArray[UnitStep[th], Automatic, 1] /.
HoldPattern[SparseArray[c___]] :> {c}[[4, 2, 2]] ]

getElementCoordinates =
Compile[{{coords, _Real, 2}, {part, _Integer, 1}},
CompileGetElement[coords, #] & /@ part,
RuntimeAttributes -> Listable];

mkCompileCode[vars_, code_, idx_] :=
With[{fVars = Flatten[vars]}, Compile[{{in, _Real, idx}},
Block[fVars,
fVars = Flatten[ in];
code]
, RuntimeAttributes -> Listable]
]

mkcForces[ chCode_] := Compile[{{p, _Real, 2}, {bars, _Integer, 2}, {fScale, _Real, 0}},
Block[{barVec, len, hBars, l0, force, forceVec},
barVec = p[[ bars[[ All, 1]] ]] - p[[ bars[[ All, 2 ]] ]];
len = Sqrt[ Plus @@@ Power[ barVec, 2]];
hBars = chCode /@ ((( Plus @@ p[[#]] ) & /@ bars )/2 );
l0 = hBars*fScale*
Sqrt[ (Plus @@ Power[ len, 2 ]) /( Plus @@ Power[ hBars, 2 ]) ];
force  = Max[ #, 0. ] & /@ (l0 - len);
forceVec = Transpose[ { #, -# } ] &[ (force / len ) * barVec ]
]
]

moveMeshPoints[ p_, bars_, cForces_, fScale_: 1.2 ] := Module[
{ totForce, forceVec},
forceVec = cForces[ p, bars, fScale];
totForce =
NDSolveFEMAssembleMatrix[ {bars,
ConstantArray[ Range[ Last[ Dimensions[p]]], {Length[ bars]}]},
Flatten[ forceVec, {{1}, {2, 3}} ], Dimensions[p]];
totForce
]

meshbarPlot[p_, bar_List] :=
Print[Show[{gr,
Graphics[Line /@ (Part[p, #] & /@ bar), AspectRatio -> Automatic,
PlotRange -> All]}]]

initialMeshPoints[ {{x1_, y1_}, {x2_, y2_} }, h0_ ] := Module[
{ xRow, yColumn, tmp },
xRow = Range[ x1, x2, h0 ];
yColumn = Range[ y1, y2, h0 *Sqrt[3]/2 ];
tmp = Transpose[ {
Flatten[
Transpose[
Table[ If[ OddQ[ i ], xRow, xRow + h0/2], { i,
Length[ yColumn ] } ] ] ],
Flatten[ Table[ yColumn, { Length[ xRow ] } ] ]
} ];
DeveloperToPackedArray@tmp
];

initialMeshPoints[ {{x1_, y1_, z1_}, {x2_, y2_, z2_} }, h0_ ] :=
Module[
{ xP, yP, zP },
xP = Range[ x1, x2, h0 ];
yP = Range[ y1, y2, h0 ];
zP = Range[ z1, z2, h0 ];
Flatten[ Outer[List, xP, yP, zP], 2]
];


Quality control functions [What is a good linear finite element]:

edgeLength[c : {{_, _} ..}] := EuclideanDistance @@@ Partition[c, 2, 1, 1];

TriangleAreaSymbolic[{c1_, c2_, c3_}] := 1/2*Det[Transpose[{c1, c2}] - c3];
triangleEdgeLength := edgeLength;
tri = Array[CompileGetElement[X, ##] &, {3, 2}];
source = 4*Sqrt[3]* TriangleAreaSymbolic[tri] / Total[triangleEdgeLength[tri]^2];

meshQuality[triangle] =
With[{code = source},
Compile[{{in, _Real, 2}}, Block[{X}, X = in; code],
RuntimeAttributes -> Listable]];

(*Shewchuk vertex ordering is a_,b_,c_,d_*)
(*this is to acomodate the mesh element ordering*)
TetrahedronVolumeSymbolic[{a_, c_, b_, d_}] :=
1/6*Det[Transpose[{a, b, c}] - d];
tetEdgeLength[{c1_, c2_, c3_, c4_}] :=
Norm[#, 2] & /@ {c1 - c2, c2 - c3, c3 - c1, c1 - c4, c2 - c4, c3 - c4}
tet = Array[CompileGetElement[X, ##] &, {4, 3}];
source = 6*Sqrt[2]*
TetrahedronVolumeSymbolic[tet]/
Sqrt[(1/6*Total[tetEdgeLength[tet]^2])]^3;

meshQuality[tetrahedra] =
With[{code = source},
Compile[{{in, _Real, 2}}, Block[{X}, X = in; code],
RuntimeAttributes -> Listable]];


Distmesh:

Options[distmesh] = {MaxIterations -> 150,
MeshElementQualityFactor -> 1/2,
MeshElementQualityFunction -> Min, ReturnBestMeshSofar -> True};

distmesh[ df_, rdf_, h0_, vars_, bbox_, pfix_, opts___ ] := Module[
{fopts, maxIterations, iterationNumber, meshQualityFactor,
meshQualityFunction,
dim,
p, ix, deps, cdfs, d,
pMid, t, bars, geps,
crdf, cForces,
mq, mqMax, elementType,
symElementCoords, qualityCode, cMeshQuality,
divisor, barSelectors,
saveBestMeshQ,
tBest, pBest,
dptol, ttol, fScale, deltat, r0, n, pOld, tmp, scaledR0, totForce,

dim = Length[ vars];

fopts = FilterRules[{opts}, Options[distmesh]];

{maxIterations, meshQualityFactor, meshQualityFunction,
saveBestMeshQ} =
{MaxIterations, MeshElementQualityFactor,
MeshElementQualityFunction, ReturnBestMeshSofar} /.
fopts /. Options[distmesh];

meshQualityFactor = Min[1, Max[0, meshQualityFactor]];
iterationNumber = 0;
mqMax = 0;

dptol = .001;
ttol = .1;
fScale = 1.2;
deltat = .2;
geps = .001*h0;

deps = Sqrt[10^-\$MachinePrecision]*h0;

(* compile the distance set fuction *)

cdfs = mkCompileCode[vars, df, 1];

(* compile raltive mesh coarsens function *)

crdf = mkCompileCode[vars, rdf, 1];
cForces = mkcForces[crdf];

(* select the mesh quality code *)
Switch[ dim,
2,
elementType = Globaltriangle;
divisor = 3;
barSelectors = {{1, 2}, {2, 3}, {3, 1}};
,
3,
elementType = Globaltetrahedra;
divisor = 4;
barSelectors = {{1, 2}, {2, 3}, {3, 1}, {1, 4}, {2, 4}, {3, 4}}
];
cMeshQuality = meshQuality[elementType];

p = initialMeshPoints[bbox, h0];
(* p = Select[ p, ( fd @@ #<geps)& ]; *)

p = p[[thresholdPosition[p, cdfs[p] - geps]]];

r0 = 1/crdf[p]^2;
scaledR0 = r0/Max[r0];
p = Join[pfix,
p[[thresholdPosition[p,
RandomReal[{0, 1}, {Length[p]}] - scaledR0 ]]]];
n = Length[p];

pOld = Infinity;
While[True,
iterationNumber++;

If[Max[Sqrt[Total[(p - pOld)^2, {2}]]/h0 ] > ttol,
pOld = p;
t = myDelaunay[p, dim];

pMid = Total[getElementCoordinates[ p, t ], {2}]/divisor;
(* t = t[[ Flatten[ Position[ fd @@@
pMid, _?(# < -geps &) ] ] ]]; *)

t = t[[thresholdPosition[ p, cdfs[ pMid] - geps] ]];
bars = Union[Sort /@ Flatten[t[[All, # ]] & /@ barSelectors, 1]];
,
Null;
];

totForce = moveMeshPoints[p, bars, cForces];
totForce[[Range[Length[pfix ]]]] =
ConstantArray[0., {Length[pfix], dim}];
p = p + deltat*totForce;
d = cdfs[p];
(*ix = Flatten[ Position[ d, _?(#>0.&) ] ];*)

ix = thresholdPosition[d, -(d + 0.)];

depsIM = deps*IdentityMatrix[dim];
Do[
(cdfs[
p[[ix]] + ConstantArray[ depsIM[[i]], {Length[p[[ix]]]}]] -
d[[ix]])/deps

If[ dim == 2 && meshBarPlotQ,  meshbarPlot[ p, bars ]; ];
mq = meshQualityFunction[
cMeshQuality[ getElementCoordinates[ p, t ] ]];
If[ (mq >= meshQualityFactor) || iterationNumber >= maxIterations,
If[ mq <= mqMax && saveBestMeshQ, mq = mqMax; p = pBest;
t = tBest; ];
Break[],
If[ mq > mqMax && saveBestMeshQ, mqMax = mq; pBest = p;
tBest = t; ];
];

];
Print["it num: ", iterationNumber, "  mq: ", mq];
Return[ { p, t } ];
]


Examples:

Example 1: Square with hole and different refinement on boundaries:

di = Function[ {x, y},
ddiff[ drectangle[ {x, y}, {{-1, -1}, {1, 1}} ],
dcircle[ {x, y}, {0, 0, 0.4} ] ] ];

h = Function[ {x, y}, Min[ 4*Sqrt[Plus @@ ({x, y}^2) ] - 1, 2 ] ];
h0 = 0.05;

bbox = {{-1, -1}, {1, 1}};
pfix = {{-1, -1}, {-1, 1}, {1, -1}, {1, 1}} // N;

AbsoluteTiming[ {coords, inci} =
distmesh[ di[x, y], h[x, y], h0, {x, y}, bbox, pfix ];]
(* it num: 87  mq: 0.66096 *)
(* {0.849295, Null} *)

Graphics[{EdgeForm[Black], FaceForm[], Polygon[(coords[[#]] & /@ inci)]}]


Example 2: Same example, finer base grid and animated:

di = Function[ {x, y},
ddiff[ drectangle[ {x, y}, {{-1, -1}, {1, 1}} ],
dcircle[ {x, y}, {0, 0, 0.4} ] ] ];
h = Function[ {x, y}, 1. ];
h0 = 0.1;
bbox = {{ -1, -1}, {1, 1} } // N;
pfix = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}} // N;
meshBarPlotQ = True;
gr = RegionPlot[
di[x, y] < 0, {x, bbox[[1, 1]], bbox[[2, 1]]}, {y, bbox[[1, 2]],
bbox[[2, 2]]}, AspectRatio -> Automatic, Frame -> False,
Axes -> False, Mesh -> False ];
distmesh[ di[x, y], h[x, y], h0, {x, y}, bbox, pfix ];
meshBarPlotQ = False;
ClearAll[gr]


Example 3: Same example, anisotropic refinement:

h = Function[ {x, y}, 1. + Abs[x*y]];
h0 = 0.03;
bbox = {{ -1, -1}, {1, 1} } // N;
pfix = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}} // N;
AbsoluteTiming[ {coords, inci} =
distmesh[ di[x, y], h[x, y], h0, {x, y}, bbox, pfix ]; ]


Example 4: A somewhat more complicated geometry:

di = Function[{x, y},
Evaluate[
Simplify[
dintersection[
ddiff[ ddiff[ dcircle[{x, y}, {0, 0, 2}],
dcircle[{x, y}, {3/2, 0, 1/5}]], dcircle[{x, y}, {0, 0, 1}]],
dtriangle[{x, y}, {{0, 0}, {2, 0}, {2, 1}}]]]]];

ci[{{x1_, y1_}, r_}, {m_, c_}] := Block[{tmpX, tmpH, tmpY},
tmpX = {(-(c*m) + x1 + m*y1 -
Sqrt[-c^2 + r^2 + m^2*r^2 - 2*c*m*x1 - m^2*x1^2 + 2*c*y1 +
2*m*x1*y1 - y1^2])/(1 + m^2), (-(c*m) + x1 + m*y1 +
Sqrt[-c^2 + r^2 + m^2*r^2 - 2*c*m*x1 - m^2*x1^2 + 2*c*y1 +
2*m*x1*y1 - y1^2])/(1 + m^2)};
tmpH = Function[{x}, m*x + c];
tmpY = tmpH /@ tmpX;
Transpose[{tmpX, tmpY}]
]

h = Function[ {x, y},
Min[ 5*Sqrt[Plus @@ ({x - 1.5, y}^2) ] - 0.2, 2 ] ];
h0 = 0.01;
bbox = { {0.8, 0}, {2, 1} } // N;
pfix = Join[ {{1., 0.}, {1.3, 0.}, {1.5, 0.2}, {1.7, 0.}, {2., 0.}},
ci[{{0, 0}, 2}, {1/2, 0}][[{2}]],
ci[{{0, 0}, 1}, {1/2, 0}][[{2}]] ] // N

AbsoluteTiming[ {coords, inci} =
distmesh[ di[x, y], h[x, y], h0, {x, y}, bbox, pfix,
MaxIterations -> 400 ];]
(*it num: 164  mq: 0.742647*)
(* {14.795667, Null} *)
Graphics[ {EdgeForm[Black], FaceForm[],
Polygon[ (coords[[#]] & /@ inci) ]} ]


Example 5: 3D:

di = Function[{x, y, z}, Sqrt[x^2 + y^2 + z^2] - 1];
h = Function[ {x, y, z}, 1. ];
h0 = 0.25;
bbox = {{ -1, -1, -1}, {1, 1, 1} } // N;
pfix = {} // N;
AbsoluteTiming[ {coords, inci} =
distmesh[ di[x, y, z], h[x, y, z], h0, {x, y, z}, bbox, pfix,
MaxIterations -> 20]; ]

TetrahedraWireframe[i_] :=
Line[ Flatten[
i[[All, #]] & /@ {{1, 2}, {2, 3}, {3, 1}, {1, 4}, {2, 4}, {3, 4}},
1]]
Graphics3D[GraphicsComplex[coords, TetrahedraWireframe[inci]]]


Example 6: Converting a black and white image into a mesh:

shape


;

Convert the shape into a distance function:

dist = DistanceTransform[Image[1 - ImageData[EdgeDetect[shape]]]];


Create an inteprolation function from it:

data = Transpose[Reverse[-ImageData[dist]*(2*ImageData[shape] - 1)]];
dataRange = Transpose[{{x, y}, {1, 1}, Dimensions[data]}];
if = ListInterpolation[ data,  dataRange[[All, {2, 3}]],
InterpolationOrder -> 1];


Call distmesh:

di = if;
h = Function[ {x, y}, 1];
h0 = 5;

bbox = Transpose[if["Domain"]];
pfix = {} // N;

AbsoluteTiming[ {coords, inci} =
distmesh[ di[x, y], h[x, y], h0, {x, y}, bbox, pfix]; ]
(* it num: 150  mq: 0.535491 *)
(* {3.215599, Null} *)


-
edgeLength[c : {{_, _} ..}] := EuclideanDistance @@@ Partition[c, 2, 1, 1] would be a bit more compact, no? – J. M. Apr 20 '13 at 19:15
@J.M., sure good suggestion. The code is a bit of a mess and there are probably many possible improvements that could be done... – user21 Apr 20 '13 at 19:19
This last example with meshing image shapes is very impressive. I successfully meshed them some years back but much less elegant. So you think distance function should define a region, not polygon line? – BoLe Apr 22 '13 at 10:47
@BoLe, thanks. I think both cases are important. The distance functions have the advantage that boolean operations like union, difference and intersection are very easy. However, the evaluation time for distance functions for complicated geometries can be too expensive for mesh generation. – user21 Apr 22 '13 at 11:51
being able to mesh images is really cool! i was thinking about how i could do this, so.. neat – Tom Apr 30 '13 at 14:14
Table[drawtriangulation[mesh @@ example, First@example,
AspectRatio -> Automatic],
{example, {circle, circle34, ellipseeye}}] // GraphicsRow


Calculating specifications for these examples:

(* distance function, bounding box, fixed points,
number of initial points, max iterations, min triangle quality *)

circle = {Sqrt[#1^2 + #2^2] - 1. &,
{{-1, -1}, {1, 1}}, {}, 100, 200, .6};

circle34 = {Max[
Sqrt[#1^2 + #2^2] - 1.,
Min[Min[Min[1. + #2, -#2], #1], 1. - #1]] &,
{{-1, -1}, {1, 1}}, {{0, -1}, {0, 0}, {1, 0}}, 100, 200, .6};

ellipseeye = {Max[
Sqrt[#1^2 + 1.6 #2^2] - 1.,
-(Sqrt[#1^2 + #2^2] - .6)] &,
{{-1, -.8}, {1, .8}}, {}, 300, 200, .6};


Following OP and the paper he's following (A simple mesh generator in MATLAB), I wrote the main code in mesh below which iteratively relaxes and moves the mesh points, and, as it is, stops at iterationmax or already when quality of the worst triangle goes over qmin; region defines geometry and can be specified as a distance function (see above; negative inside geometry only and zero at boundary) or as a list of polygon points (see update), box are two corner points of a rectangle enclosing the geometry, and npts is the number of initial points inside box.

mesh[region_List | region_Function, box_List, pfix_List,
npts_Integer, iterationmax_Integer, qmin_Real] :=
With[{ttol = .1, fscale = 1.2, deltat = .2, geps = .001},
Module[{p, dp, pold = {}, t, bars, barvec,
l, l0, fvec, ftot, x1, y1, x2, y2, outlined,
iteration = 0, nf, grad, test},
{{x1, y1}, {x2, y2}} = box;
List, nf = Nearest[region -> Automatic,
DistanceFunction -> EuclideanDistance];
outlined = outline[region, geps (x2 - x1)],
Function, grad[x_, y_] := Grad[region[x, y], {x, y}] // Evaluate];
test = If[Head[region] === List,
GraphicsMeshInPolygonQ[outlined, #] &, region @@ # < geps &];
p = pfix~Join~Select[Transpose[{
RandomReal[{x1, x2}, npts],
RandomReal[{y1, y2}, npts]}], test@# &];
While[iteration++ < iterationmax,
If[pold == {} || Max[Divide[
Norm /@ (p - pold), (x2 - x1)/Sqrt[npts]]] > ttol,
pold = p;
t = qDelaunayTriangulation[p];
t = Select[t, test@Mean[p[[#]]] &];
bars = DeleteDuplicates[Sort /@ Flatten[
Partition[#, 2, 1, 1] & /@ t, 1]]];
If[Min[quality @@@ (p[[#]] & /@ t)] > qmin, Break[]];
barvec = Subtract @@ p[[#]] & /@ bars;
l = Norm /@ barvec;
l0 = fscale Sqrt[Total[l^2]/Length@bars];
fvec = Map[Max[#, 0] &, l0 - l]/l barvec;
{#1 -> Hold[#2], Reverse[#1] -> Hold[-#2]} &,
{bars, fvec}]~Flatten~1];
If[pfix =!= {}, ftot[[Range@Length@pfix, All]] = 0];
dp = deltat Total /@ Map[ReleaseHold, ftot, {2}];
p += dp;
p = MapAt[If[Head[region] === List,
backtoedgelist[#, region, nf],
backtoedgefunction[#, region, grad]] &, p,
Position[p, {x_Real, y_Real} /; ! test[{x, y}]]];
]; p]]


You can see mesh calls qDelaunayTriangulation which is an external program, so the first thing is to mathlink Qhull which supplies this relatively speedy triangulator (@halirutan discusses it here, compiled and ready for Windows via @Oleksandr R.). It's about 3,000 times faster than DelaunayTriangulation in package ComputationalGeometry.

Providing Qhull is in the directory following the one where notebook is saved:

lnk = Install[NotebookDirectory[] <> "\\qhull64\\qh-math.exe"]


This gives access to qDelaunayTriangulation. There is also a call to InPolygonQ which is an undocumented function that determines if a point is inside a polygon or not.

Finally, the functions below: quality for testing triangle quality (1 for equilateral, 0 for a degenerate triangle with zero area), backtoedgefunction for bringing a wentover point back to boundary defined with distance function, backtoedgelist for polygon boundary, and drawtriangulation for drawing result.

quality[p1_, p2_, p3_] := With[{
a = Norm[-p1 + p2],
b = Norm[-p2 + p3],
c = Norm[-p1 + p3]},
(b + c - a) (c + a - b) (a + b - c)/(a b c)]

(* point to boundary region = 0, preevaluated gradient *)
backtoedgefunction[point_, region_Function, gradient_] :=
With[{r = point + t gradient @@ point},
r /. FindRoot[region @@ r,
{t, 0, .1}, Method -> "Secant"]]

(* point to polygon boundary, precalculated nf *)
backtoedgelist[point_, region_List, nf_NearestFunction] :=
With[{i = nf[point][[1]], d = Length@region},
Module[{a, b, c, test},
{a, b, c} = Extract[region, List /@ Which[
1 < i < d, {i - 1, i, i + 1},
i == 1, {-1, 1, 2},
i == d, {d - 1, d, 1}]];
test = And @@ Thread[{
VectorAngle[#3 - #1, #2 - #1],
VectorAngle[#3 - #2, #1 - #2]} < .5 Pi] &;
Which[
test[a, b, point], a + Projection[point - a, b - a],
test[b, c, point], b + Projection[point - b, c - b],
True, b]]]

drawtriangulation[pts_List,
region_List | region_Function, opt___] :=
With[{geps = .001},
Module[{
dt = qDelaunayTriangulation@pts,
outlined, test},
test = If[Head[region] === List,
outlined = outline[region,
geps {-1, 1}.Through[{Min, Max}[First /@ region]]];
GraphicsMeshInPolygonQ[outlined, #] &, region @@ # < geps &];
dt = Select[dt, test@Mean@pts[[#]] &];
Graphics[GraphicsComplex[pts,
{EdgeForm[Opacity[.2]], Map[{
ColorData["AlpineColors"][quality @@ pts[[#]]],
Polygon[#]} &, dt]}], opt]]]


# Update: Polygons

I extended the code to meshing polygons. (Testing whether a point is inside the region changes, also bringing the point back to the boundary changes.) Two examples:

heart = {
With[{n = 30, fx = 16 Sin[#]^3 &,
fy = 13 Cos[#] - 5 Cos[2 #] - 2 Cos[3 #] - Cos[4 #] &},
Table[{fx@t, fy@t}, {t, 0., 2 Pi (1 - 1/n), 2 Pi/n}]],
{{-16, -17}, {16, 12}}, {{0, 5}, {0, -17}}, 200, 200, .7};

slo = Module[{edge,
boundary = CountryData["Slovenia", "Polygon"][[1, 1]]},
edge = Take[boundary, {1, -1, 5}];
edge = # - Mean[edge] & /@ edge;
{edge, Transpose[{
Through[{Min, Max}[First /@ edge]],
Through[{Min, Max}[Last /@ edge]]}], {}, 200, 200, .7}];

Table[Graphics[Line@First@poly],
{poly, {heart, slo}}] // GraphicsRow


Table[drawtriangulation[mesh @@ example, First@example],
{example, {heart}}] // GraphicsRow


In both cases of region, whether it's a polygon or distance function, a geometric tolerance geps is used in testing whether a point lies inside the region. For this I "inflated" a polygon a tiny bit with these two functions:

moveout[{r1_, r2_, r3_}, d_] := Module[{
t = VectorAngle[r1 - r2, r3 - r2],
s = Normalize[r1 - r2], sign},
sign = If[Negative@Last@Cross[
(r2 - r1)~Append~0,
(r3 - r2)~Append~0], -1, 1];
r2 + Plus[d  (s /. {x_, y_} :> {y, -x}), sign d Tan[.5 (Pi - t)] s]]

outline[poly_, d_] :=
Map[Function[tri, moveout[tri, d]], Partition[poly, 3, 1, 1]]

-
+1 to implement the mesh quality! – user21 Apr 20 '13 at 6:54
+1 for tight code! – R Hall May 8 '13 at 21:14
@RHall It might be better to decide on geometry form, let's say polygon, not distance function (in my case here). That would shorten and clean the code. – BoLe May 9 '13 at 9:32
+1 for Slovenia :) – shrx Oct 27 '15 at 19:06

There is now a DistMesh port to the Wolfram Language on the Wolfram Research GitHub page.

-