How to mesh a region using adaptive cubic elements

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

• It seems to me that you are trying to re-implement the generation of a 3D mesh, which may be complicated. Do the elements have to be cubes? Because otherwise you could try DiscretizeRegion[ ImplicitRegion[{1 <= x^2 + y^2 + z^2 <= 2, x >= 0, y >= 0}, {x, y, z}], MaxCellMeasure -> 10], which gives you this. If you want more control, you may want to look at ToElementMesh in the NDSolveFEM  package. Unfortunately, I have been unable to force either method to use cubic elements (i.e. HexahedronElement). – MarcoB Mar 26 '16 at 13:16
• For reference, these questions are relevant to the (limitations of) Mathematica's FEM approach, and why it doesn't really solve your problem yet: How to create an ElementMesh of a Sphere; this answer to FEM: how to choose FEM element type specifies that HexahedronElement is only used for rectangular meshes. – MarcoB Mar 26 '16 at 13:26
• It would be lovely to be able to 1) generate a fine-grained list of your cuboids; 2) convert them to a mesh; 3) have the internal FEM machinery refine that mesh by merging elements. Unfortunately, I am not familiar enough with the internal FEM machinery to even know if it's doable, let alone do it myself :-) – MarcoB Mar 26 '16 at 13:29
• @MarcoB - thx for your hints. I will dig a bit deeper into this. Unfortunately yes, the elements have to be cubes. I'm focusing an FDTD simulation package. As far as I know, the grid is restricted to cartesian coordinates and cuboid shapes. – Kay Mar 26 '16 at 13:36
• I'm glad I could help. Thank you for the accept and bounty! – MarcoB Apr 1 '16 at 12:33

2 Answers

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


• @MaroB - First of all: nicely done, really! At a primary glimpse, it looks impressive. To my shame, I figured out, that I missed one fact in the question, but I hope it is easy to handle though. In my example (quarter of a sphere), I used an analytically defined function. The analytical function I have in mind will be treaded with a convolution first. Give me some time to figure out, if that causes me any problems (I cannot remember if the convolution was analytic or numeric). I'm on vacations at the moment, but will be back. – Kay Mar 28 '16 at 19:18
• I have one question regarding the last figure: Why is only 1 red cube visible? I mean, the adjacent cubes to the top or to the right (2x green, 4x yellow and 8x blue) are equal 1x red. Why happened a splitting over there? – Kay Mar 28 '16 at 19:32
• @Kay I'll have to look more carefully, but I suspect that's because of the slight curvature of the back side of the region, which prevents a full larger cube from fitting. Also, a numerical region should not be a problem, as long as one can generate a reasonable function that determines whether a point is or isn't within the region. – MarcoB Mar 28 '16 at 19:51
• Sounds great! I'll have more Information within the next 2-3 days (sorry for the delay). Is it possible to have some private chat or something similar, in case I have more precise questions to your solution? I mean, the comment option might be a bit inadequat... – Kay Mar 28 '16 at 19:55
• @Kay I'm not overly familiar with chat, but I think we should be able to open a chat room on the SE system to chat about this at length. I'll be looking into my proposal a little more because I think it is still relatively inefficient to refine the whole mesh at every pass, but let me know when you have more questions. – MarcoB Mar 28 '16 at 22:13

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


• This is essentially how a cartesian mesher works, isn't it? And if the OP accepted non cuboid shapes, we would then proceed to identifying to which pattern the intercepted cells belong to (by the analysis of their vertices position: inside or outside), and then substitute these cuboids by their corresponding non-hexahedral elements, fitted to best match the boundary. I wonder what are WRI plans for a (more) complete meshing framework? I would love to have a framework where I could then define a prismatic boundary-layer, etc. – P. Fonseca Mar 27 '16 at 13:01
• wow, Looks like some heavy stuff :D Before going through it, just a question regarding your prior comment.I don't get this: We start with the largest cuboids and break down the border cuboids into four smaller cuboids. In my question/code above, the cuboids are all unity, noone is larger than any other. And second, I wanted to go the other way around: merging some of the cuboids...? – Kay Mar 28 '16 at 12:03
• Kay, you can start the process with larger cube, say radius 4 or 2 and take it from there – Zviovich Mar 28 '16 at 12:31
• @Zviovich - True. Indeed, I realized that shortly after starting to study your code. Originally, I was thinking the other way around, as I'd like to define a minimum feature size first and allow than for arbitrary dimensioned cubes and cuboids (obtained by merging). Am I right: your code is restricted to, 1. I need a good guess for a starting cube-dimension, 2. only cubes/no cuboids are allowed, 3. the smaller neighbor cube is always half in size? Well 1. and 3. are some sort of the same issue, as I could guess the start dimension as 2, 4, 8, 16 and so on... I'll think about this. – Kay Mar 28 '16 at 12:55
• @Zviovich - Btw, what do the transparent cubewires represent? Outside the sphereVolume? – Kay Mar 28 '16 at 13:00