Can one use MeshRefinementFunction to obtain a nonobtuse triangulation?

Question

I have a planar convex polygonal region Rc. I would like to triangulate it such that every triangle is nonobtuse. I was hoping that the MeshRefinementFunction option would indeed "refine the simplex" when the function (fob in my case) returns T, to quote the documentation.

 R = DiscretizeRegion[Rc
      , MeshRefinementFunction -> fob
      , MeshQualityGoal -> "Maximal"
      , PerformanceGoal -> "Quality"]

But the above code does not in fact refine every triangle. E.g., 60 of the 478 triangles below are obtuse—not terribly obtuse, but obtuse (max angle approx. $113^\circ$).

---

        

---

I am aware that it is not straightforward to refine for nonobtuseness. Does the meshing function just give up after a few tries?

By request, some ugly code for fob:

---

TriLengths[{a_, b_, c_}] := Module[{A, B, C},
   A = Sqrt[(b - a).(b - a)] // N;
   B = Sqrt[(c - b).(c - b)] // N;
   C = Sqrt[(a - c).(a - c)] // N;
   Return[{A, B, C}]
   ];
TriAngles[{a_, b_, c_}] := 
  Module[{A, B, C, \[Alpha], \[Beta], \[Gamma]},
   {A, B, C} = TriLengths[{a, b, c}];
   \[Beta] = ArcCos[ (C^2 + A^2 - B^2)/(2 C A)] // N;
   \[Gamma] = ArcCos[ (A^2 + B^2 - C^2)/(2 A B)] // N;
   \[Alpha] = ArcCos[ (B^2 + C^2 - A^2)/(2 B C )] // N;
   Return[ { \[Beta], \[Gamma], \[Alpha]}];
   ];
fob[vlist_, area_] := Module[{\[Alpha], \[Beta], \[Gamma]},
   {\[Alpha], \[Beta], \[Gamma]} = TriAngles[vlist] // N;
   If[Max[{\[Alpha], \[Beta], \[Gamma]}] <= \[Pi]/2,
    Return[False],
    Return[True](*obtuse:refine*)
    ]];
(*Test*)
{a, b, c} = {{0, 0}, {0, 1}, {-1, 2}};
Print[TriAngles[{a, b, c}]];
fob[{a, b, c}, 999]

Answer 1

The following works reasonably well, up to a minimum area cutoff that is needed to limit the evaluation time:

TriLengths[{a_, b_, c_}] := 
  Module[{A, B, C}, A = Sqrt[(b - a).(b - a)];
   B = Sqrt[(c - b).(c - b)];
   C = Sqrt[(a - c).(a - c)];
   {A, B, C}];
TriAngles[{a_, b_, c_}] := 
  Module[{A, B, C, α, β, γ}, {A, B, C} = 
    TriLengths[{a, b, c}];
   β = ArcCos[(C^2 + A^2 - B^2)/(2 C A)];
   γ = ArcCos[(A^2 + B^2 - C^2)/(2 A B)];
   α = ArcCos[(B^2 + C^2 - A^2)/(2 B C)];
   {β, γ, α}];

r = Disk[{0, 0}, {3, 2}];

DiscretizeRegion[r, 
 MeshRefinementFunction -> (#2 > .005 && 
     Evaluate@Max[TriAngles[#1]] > Pi/2 &)]

I tried to clean up the TriAngles function a little, but that formal issue wasn't the reason why it didn't work.

The main thing I did was to eliminate the Module wrapper around the calls to TriAngles and add Evaluate in front of the condition.

This seems to be necessary because the parsing of MeshRefinementFunction is unable to deal with the depth of scoping constructs in the original formulation. One could say this is a bug in the parsing.

For comparison, here is what the output of the above command looks like if you omit the Evaluate:

Note the obviously obtuse triangle in the upper left boundar region.

Answer 2

Here is another approach. Mathematica uses Triangle as it's 2D mesh generator. Triangle is very efficient and returns good results for numerical routines like the Finite Element Method or interpolation functions. However, to the best of my knowledge, there is not way to tell Triangle to use not generate triangles that have a angles that are larger than a specified values. (One can request the opposite however, elements that have no angle smaller than a request value)

As another option, I'd suggest to use a different mesh generator. You can find a port of Distmesh to Mathematica (contained in the FEMAddOns) and use that. FEMAddOns is easily installed by evaluating:

ResourceFunction["FEMAddOnsInstall"][]

Distmesh is different in that it generates extremely smooth meshes. Here is an example after you install the Distmesh port:

Needs["DistMesh`"]

r = Disk[{0, 0}, {3, 2}];
em = DistMesh[r];

Distmesh returns an ElementMesh but that is trivial to convert to MeshRegion if you want.

MeshRegion[em]

There are still a few elements that do not meet the criterion:

pos = Position[
  MeshPrimitives[mr, 2] /. 
   Polygon[p_] :> (Max[TriAngles[p]] < \[Pi]/2), False]
{{2}, {10}, {18}, {175}, {196}, {260}}

HighlightMesh[mr, {2, Flatten[pos]}]

But they are not too far off:

(Extract[MeshPrimitives[mr, 2], pos] /. 
   Polygon[p_] :> Max[TriAngles[p]] - \[Pi]/2)*180/\[Pi]
{5.925868726419833`, 7.006195348307825`, 1.9880526587495373`, \
2.3708695038708707`, 5.427334986621736`, 1.4924245044698656`}

If you have suggestions for improving / extending the Distmesh port,I'd love to see a pull request ;-)

Answer 3

Something about your function TriAngles is interfering with MeshRefinementFunction. If I evaluate TriAngles inside the MeshRefinementFunction, then MeshRefinementFunction is ignored, even when the output of TriAngles isn't used.

DiscretizeRegion[Disk[], MeshRefinementFunction -> Function[{vlist, area}, N@TriAngles[vlist]; If[area < 0.001, Return[False], Return[True]]]]

But if I comment it out, MeshRefinementFunction behaves as desired.

DiscretizeRegion[Disk[], MeshRefinementFunction -> Function[{vlist, area}, N(*@TriAngles*)[vlist];If[area < 0.001, Return[False], Return[True]]]]