Creating a 2D meshing algorithm in Mathematica
Update
The code from this post is part of the FEMAddOns and can be installed by evaluating
ResourceFunction["FEMAddOnsInstall"][]
End Update
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 = Graphics`Region`DelaunayMesh[ pts];
Developer`ToPackedArray[ 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}},
Compile`GetElement[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 =
NDSolve`FEM`AssembleMatrix[ {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 ] } ] ]
} ];
Developer`ToPackedArray@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[Compile`GetElement[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[Compile`GetElement[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,
dgradDim, dgradN, depsIM, dgrad,
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,
dgradx, dgrady, dgrad2},
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 = Global`triangle;
divisor = 3;
barSelectors = {{1, 2}, {2, 3}, {3, 1}};
,
3,
elementType = Global`tetrahedra;
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.)];
dgradDim = Dimensions[p[[ix]]];
dgradN = ConstantArray[ 0., dgradDim ];
depsIM = deps*IdentityMatrix[dim];
Do[
dgradN[[All, i]] =
(cdfs[
p[[ix]] + ConstantArray[ depsIM[[i]], {Length[p[[ix]]]}]] -
d[[ix]])/deps
, {i, Last[dgradDim]}];
dgrad = Total[ dgradN^2, {2} ];
p[[ix]] = p[[ix]] - (d[[ix]]*dgradN)/dgrad;
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]]]];
ImageAdjust@dist
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} *)
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.
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;
Switch[Head[region],
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,
Graphics`Mesh`InPolygonQ[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;
ftot = SparseArray[MapThread[
{#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]]];
Graphics`Mesh`InPolygonQ[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]]