This is a follow-up question to an earlier question.

Context: A curved boundary allows a better representation of a region, as demonstrated by comparing the area of exact and numerical regions. See Documentation (“Region Approximation Quality”).

Consider a weakly non-circular boundary:

In[400]:= r[t_] = (1 - .1 Cos[2 t]);
In[438]:= \[ScriptCapitalR] = ImplicitRegion[x^2 + y^2 <= r[ArcTan[x, y]]^2, {x, y}];
nr = ToNumericalRegion[\[ScriptCapitalR]];
mesh1 = ToElementMesh[nr];
In[442]:= AreaExact = Integrate[r[\[Theta]]^2/2, {\[Theta], 0, 2 Pi}]
Out[442]= 3.1573
In[443]:= AreaExact - Total[mesh1["MeshElementMeasure"], 2]
Out[443]= 2.23124*10^-6

Because the region is defined implicitly, i.e. using ImplicitRegion, information about the curvature of the boundary is automatically passed during the call to ToElementMesh, and the mesh boundary is curved, giving a better approximation of the region.

Question: what if the boundary of the region is defined by a general space curve xy[t_] = {x[t], y[t]}? How can the region be generated to yield a mesh with a curved boundary, if ImplicitRegion cannot be used? (A related question has appeared before, but is more specialized.)

A crude approach is to extract a sequence of points, and use them as boundary nodes. This however produces a piecewise linear boundary (polygonal) with a large “error” in the area.

In[444]:= xy[t_] = r[t] {Cos[t], Sin[t]};
In[445]:= Nt = 100; Coords = 
 Table[N[xy[t]], {t, 0, 2 Pi - 2 Pi/Nt, 2 Pi/Nt}]; LElms = 
 Table[{Mod[i + 2, Nt, 1], i, i + 1}, {i, 1, Nt - 1, 2}]; bm = 
 ToBoundaryMesh["Coordinates" -> Coords, 
  "BoundaryElements" -> {LineElement[LElms]}]; mesh2 = 
 ToElementMesh[bm];
In[447]:= AreaExact - Total[mesh2["MeshElementMeasure"], 2]
Out[447]= 0.00879532

Note that the error scales quadratically with 1/Nt, whereas the implicit region error scales as the fourth order of "MaxBoundaryCellMeasure” parameter.

Assuming a curved boundary, one would guess that simply interpolating the points during the construction of the mesh boundary would be a significant improvement to the "polygonal" approximation since the order of interpolation could match the order of the mesh itself. Can this be done? If not, is there a reason why? (I understand that one sometimes requires non-differentiable points, e.g. corners, on the boundary, but surely such points can be accommodated.)

Note that I have explored using ParametricRegion:

\[ScriptCapitalR] = 
  ParametricRegion[xy[\[Rho], t], {{t, 0, 2 \[Pi]}, {\[Rho], 0, 1}}];

However, this seems to require a full coordinate mapping xy[\[Rho],t] depending on a "radial" coordinate \[Rho], i.e. it is not addressing the stated problem. Furthermore, the error seems no better than the polygonal approximation:

In[499]:= nr3 = ToNumericalRegion[\[ScriptCapitalR]];
In[500]:= mesh3 = ToElementMesh[nr3];
In[501]:= AreaExact - Total[mesh3["MeshElementMeasure"], 2]
Out[501]= 0.0033998