5
$\begingroup$

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

$\endgroup$
4
$\begingroup$
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]]

enter image description here

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

enter image description here

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]]
    ]
$\endgroup$
1
  • 1
    $\begingroup$ 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 $\endgroup$
    – at At
    Jul 22 '20 at 21:23
5
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ this works just as well as the other answers! Do you know where to find documentation on such Method options? $\endgroup$
    – at At
    Jul 22 '20 at 21:55
  • $\begingroup$ Many are undocumented. There are some Q&A on site that collect them, like this one, mathematica.stackexchange.com/a/173431/4999 $\endgroup$
    – Michael E2
    Jul 22 '20 at 22:11
3
$\begingroup$

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]

enter image description here

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

enter image description here

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"
  ]
$\endgroup$
3
  • $\begingroup$ This code is Brilliant, yet, to my understanding it would only work for axisymmetric surfaces, is it correct? What can be done for non axysimmetric surfaces? The minimal example I provided was indeed axysimmetric, but I am interested into the parametrization of a similar surface for which I have a polar expression, but that is not a surface of revolution $\endgroup$
    – at At
    Jul 22 '20 at 21:12
  • $\begingroup$ Erm. You're question was about surfaces of revolution, right? Anyways, the compiled function can generate the quads for reactangle, cylinder, and torus topologies. So just get the vertex positions somehow and set the Boolean input variables as needed. $\endgroup$ Jul 22 '20 at 21:18
  • $\begingroup$ MY BAD, somehow while posting I was rolled back to the original (incorrect) title I had written, containing the word "surfaces of revolution" instead of polar. I will try to edit now, yet your compiled function is 200% useful. $\endgroup$
    – at At
    Jul 22 '20 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.