Mathematica Asked on November 30, 2021
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
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
Answered by Michael E2 on November 30, 2021
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]]
]
Answered by halmir on November 30, 2021
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"
]
Answered by Henrik Schumacher on November 30, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP