# duplicated vertices in mesh generation

I want to generate a MeshRegion from ParametricPlot3D, e.g. the hemisphere (the PlotPoints is used to get small numbers of vertices)

g = ParametricPlot3D[{Cos[ϕ] Sin[θ], Sin[ϕ] Sin[θ], Cos[θ]}, {θ, 0, π/2},
{ϕ, 0, 2 π}, PlotPoints -> {4, 4}, MaxRecursion -> 0, Mesh -> Full];
surf = DiscretizeGraphics[Normal@g];
{g, surf}


The number of vertices, as you may see is 10, but the generated MeshRegion has 42.

MeshCells[surf, 0] // Length


Many of vertices are duplicated. In fact they are distinct but their coordinates differ very little:

Counts@Round[MeshCoordinates[surf][[All]], .0001]
Counts@Round[MeshCoordinates[surf][[All]], 10.^-16]


Is there any way to avoid such duplications?

• @Anjan - how to paste greek letters in the post ? I have searched the forum a bit but did not find. Sep 23, 2017 at 14:30
• It's a plugin. See this mathematica.meta.stackexchange.com/a/1044/19742 Sep 23, 2017 at 17:33
• it's nice, thanks Sep 24, 2017 at 20:45

Summary

I find that discretised regions often run into difficulties here (along with the related problem that many of the higher dimensional cells have essentially negligible measure). I know of no way to avoid the issue, but the following function should fix it:

fixMesh[meshregion_] := Module[{meshpts, polys},
meshpts =
Mean /@ Gather[MeshPrimitives[meshregion, 0][[;; , 1]],
Norm[#1 - #2] < 10^-5 &];
polys =
Select[Polygon /@
Map[First@Nearest[meshpts, #] &,
MeshPrimitives[meshregion, 2][[;; , 1]], {2}],
RegionMeasure[#] > 10^-8 &];
DiscretizeGraphics[Graphics3D[polys]]
]


Then

fixMesh[surf]


will return the correct mesh region with no duplicated points.

An explanation

There are two steps in the fixMesh function. The first finds the "proper" vertices of the mesh. I have just taken the mean of all points deemed sufficiently close. There are certainly more accurate ways (there's no guarantee that this mean actually lies within the region being discretised, although it won't be out by very much), but it's nice and simple.

newmeshpts = Mean /@ Gather[MeshPrimitives[surf, 0][[;; , 1]], Norm[#1 - #2] < 10^-5 &];
Length@newmeshpts


10

The cut-off 10^-5 was chosen as the smallest threshold that yields the right number of points.

The second step is to replace all the vertices in the mesh polygons (MeshPrimitives[surf, 2]) by the Nearest corrected point in newmeshpts

newpolys = Polygon /@
Map[First@Nearest[newmeshpts, #] &, MeshPrimitives[surf, 2][[;; , 1]], {2}]]


We're not done yet -- some of these polygons now have zero measure because all their vertices are the same point:

RegionMeasure /@ newpolys


{0.122818, 0.177201, 0.122818, 0.378348, 0.37835, 0.378347, 0.177201, 0., 0.306921, 0.122818, 0.306922, 0.177201, 0., 0.436876, 0.306922, 0.43688, 1.86265*10^-9, 0., 0.436879, 0., 3.72529*10^-9}

So we just Select those with non-negligible measure:

newpolys = Select[newpolys, RegionMeasure[#] > 10^-8 &];
Length@newpolys


15

Now we have all our polygons, and just need to turn it back into a mesh. A simple way to do that is

newmesh = DiscretizeGraphics[Graphics3D[newpolys]]


Looks okay, but let's check

Length@MeshPrimitives[newmesh, 0]
Length@MeshPrimitives[newmesh, 2]


10

15

• Many thanks, for the function and its expalnation. It works. Sometimes it required the proper tuning for the size of the Norm[..] and RegionMeasure[..], dependently on the particular mesh. Sep 23, 2017 at 14:33
• @sebqas I'm glad it works. Don't forget that you can accept an answer if your question has been answered to your satisfaction. Oct 5, 2017 at 7:19