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} *)
timeStepis equivalent to (the faster for more nodes)SeedRandom[2]; p0 = RandomReal[{0, 10}, {100, 2}];– gpap Mar 28 at 12:08