# Can one use MeshRefinementFunction to obtain a nonobtuse triangulation?

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]

• Giving a minimal working example would help. In particular, can you write down what fob is? – Jess Riedel Jun 5 '16 at 18:23
• @JessRiedel: I was afraid you'd ask. I'll have to untangle from my code... – Joseph O'Rourke Jun 5 '16 at 18:25
• You might find this tiny gem from Kahan useful: tang[v1_?VectorQ, v2_?VectorQ, v3_?VectorQ] := Module[{n1 = Normalize[v1 - v2], n2 = Normalize[v3 - v2]}, 2 ArcTan[Norm[n1 + n2], Norm[n1 - n2]]]. This determines the angle at vertex v2; you can then use tang @@@ Partition[pts, 3, 1, 2] to get all the angles. – J. M.'s ennui Jun 5 '16 at 20:29

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.

• Just to be clear: I didn't try it without the area cutoff because the computation ran too long for the chosen example in my answer. So the more fundamental question of whether it "can be done perfectly" in all cases may have to be left to the mathematicians. – Jens Jun 5 '16 at 20:12
• It should be mentioned that there still are obtuse triangles left in the refined mesh, so I'm not even sure if an obtuse-free meshing is possible: Position[MeshPrimitives[%, 2] /. Polygon[p_] :> (Max[TriAngles[p]] < π/2), False]. Of course, one can split the obtuse triangles right at the altitude perpendicular to their corresponding hypotenuses, but that might take more time to do. – J. M.'s ennui Jun 5 '16 at 20:24
• @J.M. Splitting obtuse triangles by the altitude requires then splitting the adjacent triangle sharing the same base, in order to achieve a proper triangulation mesh. And there can be cascades of splits caused by one split. – Joseph O'Rourke Jun 5 '16 at 20:30
• @JosephO'Rourke, the underlying triangulation algorithm is triangle. If it can be done based on that then it can be done in Mathematica. I don't think I'll implement anything outside of that. What is the application? I think there is bigger fish to fry. – user21 Jun 5 '16 at 21:54
• @JosephO'Rourke, I can not take claim for implementing Triangle - I did, however, implement the link to Triangle which you can find here. Also added another possible approach. Hope any of that helps. – user21 Jun 6 '16 at 1:14

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 ;-)

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]]]]

• I don't think the function Triangles is broken. It's just not written very cleanly, but that can't be why the function is essentially ignored. – Jens Jun 5 '16 at 19:59
• I'll rephrase "broken"-->"interfering". – Jess Riedel Jun 5 '16 at 20:00
• I have a different interpretation of what's going on, and with some changes to the MeshRefinementFunction call the original function definition can actually be made to work. – Jens Jun 5 '16 at 20:06
• @Jens: Thanks to you both for your attention. I don't see why TriAngles[] should interfere... Note that I ignore the area in the fob` function, which presumably is my prerogative. – Joseph O'Rourke Jun 5 '16 at 20:07