How to create an ElementMesh of a Sphere

Question

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[2] // 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];
assoc = AssociationThread[#, Range@Length@#] &@coords;
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]], 
   EdgeForm[Darker@Gray], Polygon@quadIncs}],
 Axes -> True, AxesLabel -> {"x", "y", "z"}, Lighting -> "Neutral"]

Equidistance

Update: Equiangular

I can also create successfully the mesh:

mesh = ToBoundaryMesh[
   "Coordinates" -> coords,
   "BoundaryElements" -> {QuadElement@quadIncs}
   ];
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[1] /@ 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?

Answer 1

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