Return to Question

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

xy[ρ_, t_] = ρ xy[t];
R = ParametricRegion[xy[ρ, t], {{t, 0, 2 Pi}, {ρ, 0, 1}}];

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

nr3 = ToNumericalRegion[R];
mesh3 = ToElementMesh[nr3];
AreaExact - Total[mesh3["MeshElementMeasure"], 2]

0.0033998

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[399]:= Needs["NDSolve`FEM`"]    
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.)

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

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

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 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 (“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.)

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

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

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

In[399]:= Needs["NDSolve`FEM`"]    
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