# How to create an ElementMesh of a Sphere

I want to create an ElementMesh of a "sphere". To reduce the computational effort I want to use HexahedronElement; the resulting mesh can have a hole at the center.

The first steps is to create a boundary ElementMesh of a spherical surface using QuadElement. Then we can project this mesh toward the center of the sphere building successive spherical shells.

What I did is to create a mesh of a circumscribed cube and then project the coordinates on the sperical surface.

Some preliminary definitions:

deepLookup[assoc_,listOfLists_,levelspec_:{-2}]:=Map[Lookup[assoc,#]&,listOfLists,levelspec];
partitionOp[spec__][list_]:=Partition[list,spec]
flattenOp[spec_][list_]:=Flatten[list,spec]
partOp[spec__][list_]:=Part[list,spec]


The surface density of the mesh (in the real case n can be large: 4000, 8000 and more for example):

n = 5;
m = 2 n + 1;
cl = Range[-1, 1, 1/n] // N


{-1., -0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8, 1.}

Update: Even better:

cl = Tan@Range[-\[Pi]/4, \[Pi]/4, \[Pi]/(4 n)] // N;


{-1., -0.726543, -0.509525, -0.32492, -0.158384, 0., 0.158384, 0.32492, 0.509525, 0.726543, 1.}

The "structured" list of coordinates on the cube:

pts = Table[
Tuples[MapAt[#[[{f}]] &, {cl, cl, cl}, {d}]], {d, 3}, {f, {-1, 1}}];
pts = pts // flattenOp // partitionOp[m] // partitionOp[m];
Dimensions[pts]


{6, 11, 11, 3}

Building the quads on each of the six faces of the cubes:

coords = DeleteDuplicates@Flatten[pts, 2];
incidents = deepLookup[assoc, pts, {-3}];
faceIncs =
incidents // Map[partitionOp[{2, 2}, {1, 1}]] //
flattenOp[{{1}, {2}, {3}, {4, 5}}] //
partOp[All, All, All, {1, 3, 4, 2}];


Now project the coordinates on the sperical surface:

coords = Normalize /@ coords;


and make a flat list of quads incidents:

quadIncs = Flatten[faceIncs, 2];


A plot of the result:

Graphics3D[
GraphicsComplex[
coords, {FaceForm[Opacity[.8, LightGray]],
Axes -> True, AxesLabel -> {"x", "y", "z"}, Lighting -> "Neutral"]


Equidistance Update: Equiangular I can also create successfully the mesh:

mesh = ToBoundaryMesh[
"Coordinates" -> coords,
];
mesh["Wireframe"] The (possible) problem is that the vertices of each QuadElement are not coplanar.

assoc2 = AssociationThread[Values@assoc, coords];
coplanarMeasure[
pts_ /; MatrixQ[pts, NumericQ] && Dimensions[pts] == {4, 3}] :=
Det[Append /@ pts]
coplanarMeasure /@ deepLookup[assoc2, quadIncs] // Histogram


Equidistance Update: Equiangular [Ancillary question. How this situation is handled by Finite Element framework? How this situation can affect the quality of the solution returned by NDSolve?]

Suppose I want to adjust the coordinates of the mesh so that the vertices of QuadElement are coplanar and the resulting mesh is the best approximation of the spherical surphace (in some sense).

How this can be done?

• there is no inherent requirement of the finite element method that quad element faces must be exactly planar. – george2079 Nov 16 '16 at 13:10

Perhaps this is what you want: The only trick is scaling the coordinates of each point to lie on the sphere inscribed in the cube symmetric about the origin circumscribing the point.

cube = ToElementMesh[Cuboid[{-1, -1, -1}, {1, 1, 1}],
MaxCellMeasure -> {"Volume" -> 0.005}];

rescale[v_] := Max[Abs@v] Normalize[v]

emesh = ToElementMesh[
"Coordinates" -> rescale /@ cube["Coordinates"],
"MeshElements" -> cube["MeshElements"]];

emesh["Wireframe"[
"MeshElementStyle" -> FaceForm[Directive[Opacity[0.7], LightGray]]]] • This is a fast and interesting approach to get the result I called "Equidistance" for a filled sphere; with some change I think I can use for my needs instead of my more complex approach. But the real problem is what I wrote at the end: the vertices of each QuadElement or the faces of HexahedronElement are not coplanar. Is this a problem for FEM? How can I adjust their positions to fix? – unlikely Jan 6 '16 at 18:19
• @unlikely It doesn't seem to be a problem for NDSolve[]. Is there an example where it is? – Michael E2 Oct 16 '16 at 23:30
• @unlikely. No, that should not be a problem. – user21 Nov 16 '16 at 1:09