(*
meshgenerator package
*)
BeginPackage["MeshGenerator`", {"NumericalMath`NLimit`", 
    "DiscreteMath`ComputationalGeometry`"}]
generateMesh::usage = "generateMesh[d_,h_,h0_,{x1_,y1_,x2_,y2_}] generates a mesh for a region defined by the distance function d(x,y) and size function h(x,y). The mesh sise is h0 and an initial rectangle region is defined by (x1,y1) and (x2,y2)."


dcircle::usage = "dcircle[{x,y},{cx,cy,r}] returns a distance function for a disck bounded by a circle with radius r and center point (cx, cy)."


drectangle::usage = "drectangle[{x,y},{x1,y1,x2,y2}] returns a distance function for a rectangle with lower left corner (x1,y1) and upper right corner (x2,y2)."


ddiff::usage = "ddiff[d1,d2] returns a distance function that takes shape d2 out of d1. d1 and d2 are distance functions defining the shapes."


Options[generateMesh] = {Fscale -> 1.2, DeltaT -> 0.2, Eps -> 0.2, 
    MaxSteps -> 50}
Begin["`Private`"] 
\!\(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]
\!\(seedMeshpoints[{x1_, y1_, x2_, y2_}, h0_] := Table[{x + \(\(1 - \((\(-1\))\)\^Quotient[y - y2, h0 \@ 3\/2]\)\/4\) h0, y}, {y, y1, y2, \(\@3\/2\) h0}, {x, x1, x2, h0}] // Flatten[#, 1] &\)
selectMeshpoints[d_, {x1_, y1_, x2_, y2_}, h0_] := 
  Select[seedMeshpoints[{x1, y1, x2, y2}, h0], d[#] < 0 &]
ddiff[d1_, d2_] := Max[d1, -d2]
NGrad[d_, {x0_, y0_}] := {ND[d[{x, y0}], x, x0], ND[d[{x0, y}], y, y0]}
backtoBoundary[d_, {x_, y_}] := 
  Module[{expr}, expr = {x, y} - s NGrad[d, {x, y}];
    expr /. (FindRoot[d[expr], {s, 0}])]
\!\(scale[x_List] := x\/Max[x]\)
\!\(scale[a_List, \ {p_, \ bars_}] := \ \@\(Plus @@ Map[barlength[p, #]\^2 &, bars]\/Plus @@ \(a\^2\)\)\[Times]a\)
\!\(middlepoint[p_, {a_, b_}] := \(p[\([a]\)] + p[\([b]\)]\)\/2\)
meshbar[t_] := 
  Union[Sort /@ (Thread[List[Sequence @@ #]] & /@ t // Flatten[#, 1] &)]
barlength[p_, {a_, b_}] := Norm[p[[a]] - p[[b]]]
ForceRule[{a_, b_}, {Fx_, Fy_}] := {{a, b} -> Hold[{Fx, Fy}],
    (* Newton's Third Law *)
     {b, a} -> Hold[{-Fx, -Fy}]}
\!\(forcemove[d_, h_, pp_, pfix_, \ Fscale_, \ deltat_] := Module[{p, t, bars, l0, F, Fvec, forcematrix, dp, newp, temp}, \[IndentingNewLine]p = Join[pfix, pp]; \[IndentingNewLine]t = DelaunayTriangulation[p]; \[IndentingNewLine] (*\ select\ only\ interia\ meshbars\ *) \[IndentingNewLine]\ bars = Select[meshbar[t], d[middlepoint[p, #]] < 0 &]; \ \ l0 = Fscale*scale[\(h[middlepoint[p, #]] &\) /@ bars, \ {p, bars}]; \ F = \(Max[#, 0] &\) /@ \((l0 - \((\(barlength[p, #] &\) /@ bars)\))\); Fvec = F \((\(Subtract @@ p[\([#]\)]\/barlength[p, #] &\) /@ bars)\); \[IndentingNewLine]forcematrix = \(SparseArray[#, {Length[p], Length[p]}] &\)@\((MapThread[ForceRule, {bars, Fvec}] // Flatten[#, 1] &)\); \[IndentingNewLine]If[pfix != {}, forcematrix[\([Range[Length[pfix]], All]\)] = 0]; \[IndentingNewLine]dp = deltat*Table[Plus @@ \((ReleaseHold /@ forcematrix[\([i]\)])\), {i, Length[p]}]; \[IndentingNewLine]newp = p + dp; \[IndentingNewLine]temp = Position[newp, _?\((d[#] > 0 &)\)]; \[IndentingNewLine]If[temp != {}, \[IndentingNewLine]newp = ReplacePart[newp, \((\(backtoBoundary[d, #] &\) /@ Extract[newp, temp])\), temp, Table[{i}, {i, Length[temp]}]]]; \[IndentingNewLine]Return[newp]\[IndentingNewLine]]\)
goodEnough[p1_, p2_, eps_] := 
  Max[Norm /@ Extract[(p2 - p1), Position[p2, _?(d[#] < 0 &)]]] < eps
\!\(generateMesh[d_, h_, h0_, {x1_, y1_, x2_, y2_}, pfix_: {}, opts___Rule] := Module[{pts, values, p, t, fscale, deltat, eps, maxsteps}, \[IndentingNewLine]fscale = \(Fscale /. {opts}\) /. Options[generateMesh]; \[IndentingNewLine]deltat = \(DeltaT /. {opts}\) /. Options[generateMesh]; \[IndentingNewLine]eps = \(Eps /. {opts}\) /. Options[generateMesh]; \[IndentingNewLine]maxsteps = \(MaxSteps /. {opts}\) /. Options[generateMesh]; \[IndentingNewLine]pts = selectMeshpoints[d, \ {x1, y1, x2, y2}, h0]; \[IndentingNewLine]values = scale[\((\(1\/h[#]\^2 &\) /@ pts)\)]; \[IndentingNewLine]p = Extract[pts, \ Position[values, \ _?\((# > Random[] &)\)]]; \[IndentingNewLine]Return[NestWhile[forcemove[d, h, #, pfix, fscale, deltat] &, p, \(! \((goodEnough[#1, #2, eps])\)\) &, 2, maxsteps]]\[IndentingNewLine]]\)
End[]
EndPackage[]