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

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.

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]]
]

• this is brilliant, it works on my non-axysimmetric surface, and it is also what I originally wanted to do, but couldn't quite wrap my head around. My question is 100% solved, thanks again Jul 22 '20 at 21:23

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

• this works just as well as the other answers! Do you know where to find documentation on such Method options? Jul 22 '20 at 21:55
• Many are undocumented. There are some Q&A on site that collect them, like this one, mathematica.stackexchange.com/a/173431/4999 Jul 22 '20 at 22:11

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;

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}}]
];