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

**Question: Given a general space curve `xy[t_] = {x[t], y[t]}`, how can a region be generated to yield a mesh with a curved boundary?  (A [related question][3] has appeared before, but is more specialized.)**

Background: An implicitly defined boundary (i.e. an exact curve) allows a better representation of a region, as compared to one generated the usual way (from a set of points connected by line segments).  This fact is demonstrated by comparing the area of exact and numerical regions.  See [Documentation][2] (“Region Approximation Quality”).

Consider, for simplicity, a circular boundary:

    Needs["NDSolve`FEM`"]    
    R = ImplicitRegion[x^2 + y^2 <= 1, {x, y}];
    nr = ToNumericalRegion[R];
    mesh1 = ToElementMesh[nr];
    AreaExact = Pi;
    AreaExact - Total[mesh1["MeshElementMeasure"], 2]

>`2.00118*10^-6`

Now compare this to building a custom mesh.  A crude approach is to extract a sequence of points, and use them as boundary nodes (let us take this number equal to the number of nodes of `mesh1`).  

    Nt = Length[DeleteDuplicates[Flatten[mesh1["BoundaryElements"][[1, 1]]]]];
    xy[t_] = {Cos[t], Sin[t]};
    coords = Table[N[xy[(2 Pi)/Nt i]], {i, 0, Nt - 1}];
    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];
    AreaExact - Total[mesh2["MeshElementMeasure"], 2]

>`0.00896404`

Because the region `R` is defined implicitly, i.e. using `ImplicitRegion`, information about the curvature of the boundary is automatically passed during the call to `ToElementMesh`, and `mesh1` is endowed with a curved boundary, making it a better approximation of the region.  For the custom mesh (`mesh2`), a piecewise linear boundary (polygonal) results, causing a large “error” in the area.

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

Assuming a `xy[t]` is smooth, 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`:

    R = 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[R];
    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