# Create a spherical Hex mesh from given surface points

For my calculations I got two spherical meshes. One is a surface mesh of quadrilateral elements and the other is the corresponding "inner" hex mesh. The problem with the meshes is, that the nodes on the surface do not fit exactly the coordinates of the other mesh. As you can see below.

So my question is, is there a simple way to just create an "inner" hex mesh from the given surface nodes?

Max

Edit:

Let me specify the problem a bit more detailed. I got the coordinates and the mesh connectivity of a boundary mesh.

bnodes = {{-100, 0, 0}, {-92.38796, -38.26835, 0}, {-92.38796,
0, -38.26835}, {-92.38796, 0, 38.26834}, {-92.38796, 38.26834,
0}, {-85.47635, -36.70011, -36.70011}, {-85.47635, -36.70011,
36.7001}, {-85.47635, 36.7001, -36.70011}, {-85.47635,
36.7001, 36.7001}, {-70.71068, -70.71068, 0}, {-70.71068,
0, -70.71068}, {-70.71068, 0, 70.71067}, {-70.71068, 70.71067,
0}, {-67.38874, -67.38874, -30.29055}, {-67.38874, -67.38874,
30.29054}, {-67.38874, -30.29055, -67.38874}, {-67.38874, \
-30.29055, 67.38873}, {-67.38874,
30.29054, -67.38874}, {-67.38874, 30.29054,
67.38873}, {-67.38874, 67.38873, -30.29055}, {-67.38874,
67.38873,
30.29054}, {-57.73503, -57.73503, -57.73503}, {-57.73503, \
-57.73503, 57.73502}, {-57.73503,
57.73502, -57.73503}, {-57.73503, 57.73502,
57.73502}, {-38.26835, -92.38796, 0}, {-38.26835,
0, -92.38796}, {-38.26835, 0, 92.38795}, {-38.26835, 92.38795,
0}, {-36.70011, -85.47635, -36.70011}, {-36.70011, -85.47635,
36.7001}, {-36.70011, -36.70011, -85.47635}, {-36.70011, \
-36.70011, 85.47634}, {-36.70011,
36.7001, -85.47635}, {-36.70011, 36.7001,
85.47634}, {-36.70011, 85.47634, -36.70011}, {-36.70011,
85.47634,
36.7001}, {-30.29055, -67.38874, -67.38874}, {-30.29055, \
-67.38874, 67.38873}, {-30.29055,
67.38873, -67.38874}, {-30.29055, 67.38873,
67.38873}, {0, -100,
0}, {0, -92.38796, -38.26835}, {0, -92.38796,
38.26834}, {0, -70.71068, -70.71068}, {0, -70.71068,
70.71067}, {0, -38.26835, -92.38796}, {0, -38.26835,
92.38795}, {0, 0, -100}, {0, 0, 100}, {0,
38.26834, -92.38796}, {0, 38.26834, 92.38795}, {0,
70.71067, -70.71068}, {0, 70.71067, 70.71067}, {0,
92.38795, -38.26835}, {0, 92.38795, 38.26834}, {0, 100,
0}, {30.29054, -67.38874, -67.38874}, {30.29054, -67.38874,
67.38873}, {30.29054, 67.38873, -67.38874}, {30.29054,
67.38873,
67.38873}, {36.7001, -85.47635, -36.70011}, {36.7001, \
-85.47635,
36.7001}, {36.7001, -36.70011, -85.47635}, {36.7001, \
-36.70011, 85.47634}, {36.7001, 36.7001, -85.47635}, {36.7001,
36.7001, 85.47634}, {36.7001, 85.47634, -36.70011}, {36.7001,
85.47634, 36.7001}, {38.26834, -92.38796, 0}, {38.26834,
0, -92.38796}, {38.26834, 0, 92.38795}, {38.26834, 92.38795,
0}, {57.73502, -57.73503, -57.73503}, {57.73502, -57.73503,
57.73502}, {57.73502, 57.73502, -57.73503}, {57.73502,
57.73502,
57.73502}, {67.38873, -67.38874, -30.29055}, {67.38873, \
-67.38874,
30.29054}, {67.38873, -30.29055, -67.38874}, {67.38873, \
-30.29055, 67.38873}, {67.38873,
30.29054, -67.38874}, {67.38873, 30.29054,
67.38873}, {67.38873, 67.38873, -30.29055}, {67.38873,
67.38873, 30.29054}, {70.71067, -70.71068, 0}, {70.71067,
0, -70.71068}, {70.71067, 0, 70.71067}, {70.71067, 70.71067,
0}, {85.47634, -36.70011, -36.70011}, {85.47634, -36.70011,
36.7001}, {85.47634, 36.7001, -36.70011}, {85.47634, 36.7001,
36.7001}, {92.38795, -38.26835, 0}, {92.38795,
0, -38.26835}, {92.38795, 0, 38.26834}, {92.38795, 38.26834,
0}, {100, 0, 0}}

Inc = {{72, 67, 52, 50}, {67, 61, 54, 52}, {88, 83, 67, 72}, {83, 77,
61, 67}, {96, 93, 83, 88}, {93, 85, 77, 83}, {98, 97, 93,
96}, {97, 89, 85, 93}, {85, 69, 61, 77}, {69, 56, 54, 61}, {89,
73, 69, 85}, {73, 57, 56, 69}, {52, 35, 28, 50}, {35, 19, 12,
28}, {54, 41, 35, 52}, {41, 25, 19, 35}, {56, 37, 41, 54}, {37,
21, 25, 41}, {57, 29, 37, 56}, {29, 13, 21, 37}, {21, 9, 19,
25}, {9, 4, 12, 19}, {13, 5, 9, 21}, {5, 1, 4, 9}, {28, 33, 48,
50}, {33, 39, 46, 48}, {12, 17, 33, 28}, {17, 23, 39, 33}, {4, 7,
17, 12}, {7, 15, 23, 17}, {1, 2, 7, 4}, {2, 10, 15, 7}, {15, 31,
39, 23}, {31, 44, 46, 39}, {10, 26, 31, 15}, {26, 42, 44,
31}, {48, 65, 72, 50}, {65, 81, 88, 72}, {46, 59, 65, 48}, {59,
75, 81, 65}, {44, 63, 59, 46}, {63, 79, 75, 59}, {42, 70, 63,
44}, {70, 86, 79, 63}, {79, 91, 81, 75}, {91, 96, 88, 81}, {86,
94, 91, 79}, {94, 98, 96, 91}, {95, 92, 97, 98}, {92, 84, 89,
97}, {87, 82, 92, 95}, {82, 76, 84, 92}, {84, 68, 73, 89}, {68,
55, 57, 73}, {76, 60, 68, 84}, {60, 53, 55, 68}, {71, 66, 82,
87}, {66, 60, 76, 82}, {49, 51, 66, 71}, {51, 53, 60, 66}, {55,
36, 29, 57}, {36, 20, 13, 29}, {53, 40, 36, 55}, {40, 24, 20,
36}, {20, 8, 5, 13}, {8, 3, 1, 5}, {24, 18, 8, 20}, {18, 11, 3,
8}, {51, 34, 40, 53}, {34, 18, 24, 40}, {49, 27, 34, 51}, {27, 11,
18, 34}, {3, 6, 2, 1}, {6, 14, 10, 2}, {11, 16, 6, 3}, {16, 22,
14, 6}, {14, 30, 26, 10}, {30, 43, 42, 26}, {22, 38, 30, 14}, {38,
45, 43, 30}, {27, 32, 16, 11}, {32, 38, 22, 16}, {49, 47, 32,
27}, {47, 45, 38, 32}, {43, 62, 70, 42}, {62, 78, 86, 70}, {45,
58, 62, 43}, {58, 74, 78, 62}, {78, 90, 94, 86}, {90, 95, 98,
94}, {74, 80, 90, 78}, {80, 87, 95, 90}, {47, 64, 58, 45}, {64,
80, 74, 58}, {49, 71, 64, 47}, {71, 87, 80, 64}};

bmesh = ToBoundaryMesh["Coordinates" -> bnodes,
bmesh["Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm[LightBlue],
PlotRange -> {Automatic, {0, 100}, Automatic}]]


That I have to fill with a solid body, which I created using.

ball = AddMeshMarkers[BallMesh[{0.0, 0.0, 0.0}, 100, 4],
"MeshElementsMarker" -> 1];


The problem is that the boundary mesh has to cover the solid mesh exactly and as you can see in the very first picture, (red nodes: boundary mesh, blue: solid mesh) there are some nodes that do not fit the solid mesh.

So I was wandering if there is a simple way to create the solid mesh, that you can see in the third picture, starting from the given nodes and mesh connectivities from the boundary mesh. So that both meshes fit in the end perfectly.

I hope the explanation is now a bit more precise. :)

• I'm not sure exactly what conditions you need on the output, and I don't think I understand yet what problem the inner/outer distinction and the misalignment of vertices introduces, but how far does ConvexHullMesh go towards solving your problem? Commented May 11, 2021 at 22:50

# Update

In view of the updated question, the following workflow may be appropriate. If we assume that the boundary mesh is close to the boundary mesh produced by BallMesh, then we can use the Nearest function to match the nodes that are close and simply replace them. For example:

nf = Nearest[ball["Coordinates"] -> "Index"];
newcrd = ball["Coordinates"];
newcrd[[Flatten@nf[bnodes]]] = bnodes;
ball2 = ToElementMesh["Coordinates" -> newcrd,
"MeshElements" -> ball["MeshElements"]];
ball2["Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm[LightBlue],
PlotRange -> {Automatic, {0, 100}, Automatic}]]


First, let me note a specific function within the MeshTools package called SphericalShellMesh that will give you control over the number of layers of a spherical shell hexahedron element mesh. The following will show you how to create a layer of hexahedron elements given a quad mesh.

First, generate a Quad mesh of a spherical shell.

(*Import required FEM package*)
Needs["NDSolveFEM"];
(*Install MeshTools*)
(*Uncomment if not installed*)
(*ResourceFunction["GitHubInstall"]["c3m-labs","MeshTools"]*)
Needs["MeshTools"]
bmesh = SphereMesh[6];
bmesh["Wireframe"[
"MeshElementStyle" -> FaceForm[LightBlue]]
]
bmesh["Wireframe"[
"MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> FaceForm[LightBlue],
PlotRange -> {Automatic, {0, 1}, Automatic}]
]


Second, use this code to extrude a single layer of hexahedron elements at 90% of the initial radius.

crd = bmesh["Coordinates"];
inc = ElementIncidents[bmesh["BoundaryElements"]][[1]];
ncrd = Length@crd;
crd = Join[crd, ScalingTransform[0.9 {1, 1, 1}]@crd];
hexinc = inc /. {{i_, j_, k_, l_} -> {i, j, k, l, i + ncrd, j + ncrd,
k + ncrd, l + ncrd}};
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {HexahedronElement[hexinc]}];
mesh["Wireframe"[
"MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm[LightBlue],
PlotRange -> {Automatic, {0, 1}, Automatic}]
]


• Thats not really the answer of my question, but I think it was hard to guess what I really wanted. So I added some more information to the main question.
– Max
Commented May 12, 2021 at 8:33
• @Max I think the only straightforward way to convert a general quad surface mesh to a volume mesh in Mathematica is to convert quads into tris and create a Tet mesh. Higher-end proprietary meshers may allow you to create a hex-only volume mesh, but the underlying Mathematica mesher is TetGen. Therefore, it is designed for creating tetrahedral meshes and not hex meshes. Commented May 12, 2021 at 14:09
• @Max I updated the answer assuming that the surface mesh produced by BallMesh  is similar to bmesh`. The approach will not work on an arbitrary unstructured Quad mesh. Commented May 12, 2021 at 17:09
• Thank you very much for this answers!
– Max
Commented May 17, 2021 at 12:02
• @Max Thank you for accepting my answer! I hope you will find the technique useful. Commented May 19, 2021 at 3:25