How to discretize a parametric polar surface ensuring the correct periodicity of the mesh

Question

I am trying to discretize a polar surface using a combination of ParametricPlot3D and DiscretizeGraphics. Unfortunately, the mesh facets refuse to connect along the u=0, u=2Pi line. This is clearly shown by FindMeshDefects.

test = DiscretizeGraphics[
  ParametricPlot3D[{Cos[u], Sin[u], v}, {u, 0, 2 \[Pi]}, {v, 0, 1}, 
   PlotPoints -> {155, 20}, MaxRecursion -> 0, Mesh -> None, 
   MeshStyle -> None]]

test // FindMeshDefects

I know that specialized tools exist for surfaces of revolution, but this is meant to be just a minimal example. I also know of approaches using Implicit Regions, that I would like to avoid due to the bad quality of the meshing it produces.

I realize that the problem is due to the fact that DiscretizeGraphics computes distinct edges for u=0 and u=2Pi, yet I believe that it might be possible to identify said edges using a small distance criterion, yet I fail to do so algorithmically.

thanks in advance for the advice

Answer 1

mesh = DiscretizeGraphics[
   ParametricPlot3D[{Cos[u], Sin[u], v}, {u, 0, 2 \[Pi]}, {v, 0, 1}, 
    PlotPoints -> {155, 20}, MaxRecursion -> 0, Mesh -> None, 
    MeshStyle -> None]];

Get a mesh connectivity graph and find candidate edges (or points):

g = MeshConnectivityGraph[mesh, {1, 1}, 2];

bcells = Pick[VertexList[g], VertexDegree[g], 2];

bpoints = 
  DeleteDuplicates[
   Flatten[MeshPrimitives[mesh, bcells][[All, 1]], 1]]; 

Graphics3D[Point[bpoints]]

Then find pair of points that are close each other:

nfunc = Nearest[bpoints];

prules = Rule @@@ 
   DeleteDuplicates[
    With[{p = nfunc[#, 2][[2]]}, 
       If[Norm[# - p] < 10^-5, Sort[{#, p}], Nothing]] & /@ bpoints];

And construct new mesh:

nmesh = MeshRegion[MeshCoordinates[mesh] /. prules, 
  MeshCells[mesh, 2]];

FindMeshDefects[nmesh, "HoleEdges"]

You could make function to do this all together:

stitchMesh[mesh_, delta_:10^-5] :=
    Block[{g, bcells, bpoints, nfunc, prules},
        g = MeshConnectivityGraph[mesh,{1,1},2];
        bcells = Pick[VertexList[g],VertexDegree[g],2];
        bpoints = DeleteDuplicates[Flatten[MeshPrimitives[mesh,bcells][[All,1]],1]];
        nfunc = Nearest[bpoints];
        prules = Rule@@@DeleteDuplicates[With[{p=nfunc[#,2][[2]]},If[Norm[#-p]< delta,Sort[{#,p}], Nothing]]& /@ bpoints];
        MeshRegion[MeshCoordinates[mesh]/.prules, MeshCells[mesh,2]]
    ]

Answer 2

The gap in the edge can be eliminated in the plot with the method option ”BoundaryOffset”:

test = DiscretizeGraphics[
  ParametricPlot3D[{Cos[u], Sin[u], v},
   {u, 0, 2 \[Pi]}, {v, 0, 1}, 
   PlotPoints -> {155, 20}, MaxRecursion -> 0,
   Mesh -> None, MeshStyle -> None,
   Method -> “BoundaryOffset” -> False]]

test // FindMeshDefects

Answer 3

Instead of

f = 1. + .5 Sin[4 Pi #] &;
ParametricPlot3D[{f[v] Cos[u], f[v] Sin[u], v}, {u, 0, 2 \[Pi]}, {v, 
  0, 1}, PlotPoints -> {155, 20}, MaxRecursion -> 0, Mesh -> None, 
 MeshStyle -> None]

you can just do

f = 1. + .5 Sin[4 Pi #] &;
n = 155;
{x, y} = Transpose@Cases[
     Plot[f[v], {v, 0, 1}, PlotPoints -> 20],
     _Line,
     \[Infinity]
     ][[1, 1]];
m = Length[x];
\[Theta] = Most@Subdivide[0., 2. Pi, n];
pts = Join @@ Transpose[{
     Transpose[ConstantArray[x, Length[\[Theta]]]],
     KroneckerProduct[y, Cos[\[Theta]]],
     KroneckerProduct[y, Sin[\[Theta]]]
     },
    {3, 1, 2}
    ];
{q1, q2, q3, q4} = Transpose[getGridQuads[n + 1, m, True, False]];

R = MeshRegion[pts, Triangle[Join[Transpose[{q1, q2, q3}], Transpose[{q3, q4, q1}]]]]

where

getGridQuads = Compile[{
   {m, _Integer}, {n, _Integer},
   {xclosed, True | False}, {yclosed, True | False}
   },
  Block[{a1, a2, a3, a4, b1, b2, quads, qq, mm, nn},
   b1 = Boole[xclosed];
   b2 = Boole[yclosed];
   mm = m - b1;
   nn = n - b2;
   
   quads = Flatten[Table[
      qq = Table[
        a1 = mm (j - 1) + i;
        a2 = mm (j - 1) + i + 1;
        a3 = mm j + i;
        a4 = mm j + i + 1;
        {a1, a2, a4, a3},
        {i, 1, mm - 1}];
      
      If[xclosed,
       Join[qq,
        a1 = mm (j - 1) + mm;
        a2 = mm (j - 1) + 1;
        a3 = mm (j) + mm;
        a4 = mm (j) + 1;
        {{a1, a2, a4, a3}}
        ],
       qq
       ]
      ,
      {j, 1, nn - 1}], 1];
   
   If[yclosed,
    qq = Table[
      a1 = mm (nn - 1) + i;
      a2 = mm (nn - 1) + i + 1;
      a3 = i;
      a4 = i + 1;
      {a1, a2, a4, a3},
      {i, 1, mm - 1}];
    If[xclosed,
     a1 = mm nn;
     a2 = mm (nn - 1) + 1;
     a3 = mm;
     a4 = 1;
     qq = Join[qq, {{a1, a2, a4, a3}}]
     ];
    Join[quads, qq],
    quads
    ]
   ],
  RuntimeOptions -> "Speed"
  ]