This is a follow-up question to an [earlier question][1].

Context: A curved boundary allows a better representation of a region, as demonstrated by comparing the area of exact and numerical regions.  See [Documentation][2] (“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][3] 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

  [1]: https://mathematica.stackexchange.com/questions/149178/discrepancy-with-volume-of-two-differently-generated-finite-element-meshes
  [2]: http://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html
  [3]: https://mathematica.stackexchange.com/questions/125598/how-to-define-space-inside-a-closed-curve-as-a-region