Can one use MeshRefinementFunction to obtain a nonobtuse triangulation?

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.

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