How to mesh a region using adaptive cubic elements

Question

Assume a list of {x,y,z} coordinates. The coordinates are arranged in a regular fashion through a fixed raster size (called raster in the code below). For simplicity I visualized the list by applying Cuboid while {x,y,z} is now represented by the center of the cuboid. I found that a quarter of a sphere is a perfect shape to test and applied some boundary conditions to the list to fulfill that purpose. This quarter-sphere now consists of single cuboids with unit dimension (figure under the code below).

For my purpose, the surface of this quarter-sphere is fairly rough. To solve this problem, one could decrease the raster, but only at the expense of an increase in the total number of cuboids (raster = 1 -> $\approx 10^3$ cuboids, raster = .2 -> $\approx 10^5$ cuboids), which causes higher loading times, less handling ability, and so on.

However, allowing cuboids with non-unity dimension (here defined by cubeDim) could fix this problem; this would be equivalent to merging several of the smaller cuboids to a single cuboid with larger dimension.

Unfortunately, I do not have a clue how to do that with the results of my code below.

range = 15;
radiusMax = 10;
radiusMin = 6;
raster = 2;

(* -- define sphere volume -- *)
sphereVolume[x_, y_, z_] := 
  radiusMin <= Sqrt[x^2 + y^2 + z^2] <= radiusMax;

(* -- table one's (sphereVolume) and zero's (outside of sphereVolume) -- *)  
data = Table[
   If[sphereVolume[x, y, z], 1, 0], {x, 0, range, raster}, {y, 0, 
    range, raster}, {z, -range, range, raster}];

(* -- get position of one's, serving as position for cube center (for visualization) -- *)
cubeCenter = Position[data, 1];
cubeDim = 1;

(* -- list of all cubes -- *)
cubeList = 
  Table[Cuboid[{cubeCenter[[a]] - cubeDim/2, 
     cubeCenter[[a]] + cubeDim/2}], {a, 1, Length@cubeCenter}];

Graphics3D[cubeList, Boxed -> False]

My Question :

How to enforce Mathematica to find a solution merging some of the cubes? It isn't necessary to find the ultimate best solution, but drastically reducing the number of elements (cuboids) will greatly help me.

Remember that the cubes just serve for visualization. In the end, I want to make a split from {x, y, z} coordinates (with unit distances) to {x, y, z} coordinates (representing the center of a cuboid) now linked to the dimension of the particular non-uniform cube.

I hope my question is clear. Just give me a hint, if something is unclear. In addition, if you have a better suggestion for the question title, suggestions are appreciated.

Edit:

I realized the need to be more precise regarding my problem starting point, as two answers (some kind of similar in their approach) have been uploaded. In my example above (the quarter-sphere) I first defined an analytical function. Then I transferred the analytical function into a numerical list by applying Table. This list or array (if you so will), has to be the starting point. There is no option for me, to shift this transfer, as the numerical list will undergo a convolution, which is only numerically solvable.

Therefore, could you please be especially focusing on the point that I have to start with either an array of zero’s and one’s or a list of {x, y, z}-coordinates with a fixed distance dx, dy, dz (being the unrefined meshgrid somehow). Thanks!!

Answer 1

I discarded my previous approach to generate cubes, then fuse them together, since it seems to do a lot of wasted work. Instead, I propose here my version of a cartesian mesher.

The approach is conceptually the same as the one delineated by Zviovich, but I wasn't entirely satisfied with his results, as it seems to me that his process still leads to significant over-refinement of the cells.

---

OP asked for a method that employed larger cubical cells in the core of the region to minimize the number of cells generated.

The choice to refine a cell is made based on these considerations:

• is each of the vertices of the cell already entirely within the region to be meshed? --> no further refinement is necessary; the cell is left as is.
• is the cell entirely outside of the region? --> then that cell is dropped from the mesh.
• if the cell is partially within the region, its volume is calculated; if it is under a user-defined cutoff, then no further refinement is necessary; otherwise, the cells is split into eights.

The refinement process is repeated until no further changes are possible to the mesh according to these rules (using FixedPoint).

I first defined a smaller, but still interesting sub-region of the one in the OP to work on:

region = ImplicitRegion[
    {6 <= Norm[{x, y, z}] <= 13, 0 <= x <= 8, 0 <= y <= 4, 0 <= z <= 8}, 
    {x, y, z}
  ];

DiscretizeRegion[region]

Here are some helper functions:

rmf = RegionMember[region];

(* A function to split a cube cell into eights *)
Clear[cubeSplit]
cubeSplit[{p0_, p1_}] :=
 Module[
   {center, edge, newbasecube, baselayer, toplayer},
   center = Mean[{p0, p1}];
   edge = EuclideanDistance[p1, p0]/Sqrt[3];
   newbasecube = {p0, center};
   baselayer = 
     SortBy[Norm] /@ 
       Table[RotationTransform[n Pi/2, {0, 0, 1}, center] /@ newbasecube, {n, 0, 3}];
   toplayer = 
     baselayer /. {x_, y_, z_} -> ({x, y, z} + {0, 0, Sign[(p1 - p0)[[3]]] edge/2});
   Flatten[{baselayer, toplayer}, 1]
 ]

Below is the cell refinement function which implements the logic described above. The function takes the coordinates of the opposite vertices of a cube as descriptors of the cell, a function object returning whether a point is within the region to be meshed (True) or not (False), and a user-defined mesh quality parameter, here using the maximum admissible cell volume for those cells that are at the region boundary:

Clear[refineCube]
refineCube[{p0_, p1_}, regionmemberfun_, maxvol_] := Module[
  {vertices, verteval},

  (*generate all vertices*)
  vertices = 
    Table[RotationTransform[n Pi/2, {0, 0, 1}, Mean[{p0, p1}]] /@ {p0, p1}, {n, 0, 3}];
  verteval = Flatten@regionmemberfun@vertices;

  Which[
    (*all False: cube must be completely out of region; drop it*)
    Nor @@ verteval, Nothing,

    (*all True: all vert are part of region, cube is already refined*)
    And @@ verteval, {p0, p1},

    (*at least one vertex in, but not all: consider area and split*)
    True,
    If[
      Abs[Times @@ (p1 - p0)] > maxvol,(*vol of cube > maxvol*)
      Sequence @@ cubeSplit[{p0, p1}],
      {p0, p1}
    ]
  ]
]

Here is the refinement process:

maxCellVolume = 0.05;

cubemeshpoints = FixedPoint[
     Function[{cube}, Map[refineCube[N@#, rmf, maxCellVolume] &, cube]],
     cubeSplit@{{0, 0, 0}, {13, 13, 13}}
   ];

Graphics3D[
  {Red, Cuboid @@@ cubemeshpoints},
  Axes -> True, Lighting -> "Neutral"
]

Notice that some internal cells have been allowed to remain larger, so as to minimize the number of cells in the mesh while at the same time adaptively reducing the cell size at the boundaries to improve the fit there.

This property can also be visualized as a histogram distribution of the cell volumes in the mesh:

Histogram[
  Abs[Times @@ (#2 - #1)] & @@@ cubemeshpoints, "Log", "LogCount",
  FrameStyle -> Directive[Black, 16],
  Frame -> {True, True, False, False}, 
  FrameLabel -> (Style[#, 20] & /@ {"cell volume", "no. of cells"}),
  ChartStyle -> {Opacity[0.8], ColorData[97, "ColorList"]}, 
  ImageSize -> Large
]

Here is a similarly color-coded representation of this distribution on the calculated mesh:

Graphics3D[
  MapThread[
    {#1, Cuboid @@@ #2} &,
    {
      ColorData[97, "ColorList"][[1 ;; 4]],
      GatherBy[cubemeshpoints, Abs[Times @@ (#[[2]] - #[[1]])] &]
    }
  ],
  Lighting -> "Neutral", ImageSize -> Large
]

---

Here is the method applied to a slightly different region, a solid arch, to highlight the difference in scale among the cubes comprising the mesh:

Clear[region, rmf]
region = ImplicitRegion[
    {6 <= Norm[{x, z}] <= 13, 0 <= x <= 8, 0 <= y <= 4, 0 <= z <= 8},
    {x, y, z}
  ];
rmf = RegionMember[region];

maxCellVolume = 0.05;
cubemeshpoints = FixedPoint[
     Function[{cube}, Map[refineCube[#, rmf, maxCellVolume] &, cube]],
     cubeSplit@{{0., 0., 0.}, {8, 8, 8}}
   ];

Graphics3D[
 {Opacity[0.5], EdgeForm[Opacity[0.15]],
  MapThread[
   {#1, Cuboid @@@ #2} &,
   {
    ColorData[97, "ColorList"][[1 ;; 4]],
    GatherBy[cubemeshpoints, Abs[Times @@ (#[[2]] - #[[1]])] &]
   }
  ]},
 Lighting -> "Neutral", ImageSize -> Large
]

---

Finally, we can apply this method to the region originally mentioned in the OP:

region = ImplicitRegion[{6 <= Norm[{x, y, z}] <= 10, x >= 0, y >= 0}, {x, y, z}];
rmf = RegionMember[region];

maxCellVolume = 0.2;

cubemeshpoints = FixedPoint[
   Function[{cube}, Map[refineCube[#, rmf, maxCellVolume] &, cube]],
   cubeSplit@{{-10., -10., -10.}, {10., 10., 10.}}
 ];
Graphics3D[{Red, Cuboid @@@ cubemeshpoints}, Axes -> True, Lighting -> "Neutral"]

In this case it is difficult to visualize the larger interior cubes since they are surrounded by the smaller exterior ones. Instead, let's slice through the region with a horizontal plane:

frames = ParallelTable[
   Graphics3D[
    MapThread[
     {#1, Cuboid @@@ #2} &,
     {
      ColorData[97, "ColorList"][[1 ;; 3]],
      GatherBy[cubemeshpoints, Round[Abs[Times @@ (#[[2]] - #[[1]])], 1*^-4] &]
      }
    ],
    Lighting -> "Neutral", ImageSize -> Large,
    ClipPlanes -> {{0, 0, -1, zmax}}
   ],
   {zmax, 10.25, 0, -0.25}
  ];

Answer 2

My take on the problem. We start with the largest cuboids and break down the border cuboids into four smaller cuboids. The function can be nested as many time as needed.

rMin = 6;
rMax = 10;

(*--define sphere volume--*)
sphereVolume[{x_, y_, z_}] := rMin <= Sqrt[x^2 + y^2 + z^2] <= rMax;

(*Our Data Structure will look as follows
element={{x0,y0,z0},distanceToWall,volumeBoundQ}
*)

(*--Function to see if any vertex of the cuboid touches the region \
function *)
vertexBound[element_, regionFunction_] := 
 Module[{x, y, z, d}, {x, y, z} = element[[1]]; d = element[[2]]; 
  Fold[Or, False, 
   regionFunction[#] & /@ 
    Flatten[Table[{i, j, k}, {i, x - d, x + d, 2 d}, {j, y - d, y + d,
        2 d}, {k, z - d, z + d, 2 d}], 2]]]

(*--Function to see if cuboid is inside the region function *)
completeBound[element_, regionFunction_] := 
 Module[{x, y, z, d}, {x, y, z} = element[[1]]; d = element[[2]]; 
  Fold[And, True, 
   regionFunction[#] & /@ 
    Flatten[Table[{i, j, k}, {i, x - d, x + d, 2 d}, {j, y - d, y + d,
        2 d}, {k, z - d, z + d, 2 d}], 2]]]

(*--Break Cuboids touching the region function into smaller cuboids--*)


explode[element_, regionFunction_] := 
 Module[{x, y, z, distance = element[[2]]/2}, {x, y, z} = 
   element[[1]]; 
  Select[Flatten[
    Table[{{i, j, k}, distance, 0}, {i, x - distance, x + distance, 
      2 distance}, {j, y - distance, y + distance, 2 distance}, {k, 
      z - distance, z + distance, 2 distance}], 2], 
   vertexBound[#, regionFunction] &]]

explodeList[list_, regionFunction_] := 
 Flatten[explode[#, regionFunction] & /@ list, 1]

show3D[validCentroids_] := 
 Graphics3D[{Opacity[0], 
   Cuboid[#[[1]] - #[[2]], #[[1]] + #[[2]]] & /@ 
    Select[validCentroids,  ! #[[3]] &], Opacity[1], Red, 
   Cuboid[#[[1]] - #[[2]], #[[1]] + #[[2]]] & /@ 
    Select[validCentroids, #[[3]] &]}]

(*--Process Data Structure--*)
refineMesh[cuboidList_, regionFunction_] := 
 Module[{completeSet, explodeSet, set, unprocessed}, 
  completeSet = Select[cuboidList, #[[3]] &]; 
  explodeSet = 
   explodeList[Select[cuboidList, ! #[[3]] &], regionFunction]; 
  set = Join[completeSet, explodeSet]; 
  unprocessed = Flatten[Position[set, {___, 0}]]; 
  set[[unprocessed, 3]] = 
   completeBound[#, regionFunction] & /@ set[[unprocessed]]; set]

validCentroids = 
  Select[Flatten[
    Table[{{i, j, k}, 0.5, False}, {i, 1, rMax}, {j, 1, rMax}, {k, 1, 
      rMax}], 2], vertexBound[#, sphereVolume] &];

validCentroids = 
  Nest[refineMesh[#, sphereVolume] &, validCentroids, 3];

show3D[Cases[validCentroids, {{_, a_, _}, ___} /; Abs[a - 4] < 0.5]]