# 3D Inclusion with structured mesh and coarse and arbitrary matrix

I am wondering if there is a simple way to define a structured 3D mesh (inclusion) like e.g. that

and surround it with a corse and unstructured Matrix. It should be possible to refine both more or less independently (of course when the inclusion is refined the outer matrix will also get finer at the direct inclusin-matrix-boundary).

My attempts so fare always results in an unstructured mesh for the inclusion as well as for the matrix like this:

I have used the following code to generate the last mesh...

xI = 200; yI = 200; zI = 20;
InclusionRegion =
Region[Hexahedron[{{-xI/2, -yI/2, -zI/2}, {xI/2, -yI/2, -zI/2}, {xI/
2, yI/2, -zI/2}, {-xI/2, yI/2, -zI/2}
, {-xI/2, -yI/2, zI/2}, {xI/2, -yI/2, zI/2}, {xI/2, yI/2,
zI/2}, {-xI/2, yI/2, zI/2}}], Axes -> True];

xM = xI*2; yM = yI*2; zM = zI*2;
MatrixRegion =
Region[Hexahedron[{{-xM/2, -yM/2, -zM/2}, {xM/2, -yM/2, -zM/2}, {xM/
2, yM/2, -zM/2}, {-xM/2, yM/2, -zM/2}
, {-xM/2, -yM/2, zM/2}, {xM/2, -yM/2, zM/2}, {xM/2, yM/2,
zM/2}, {-xM/2, yM/2, zM/2}}], Axes -> True];

mesh = ToElementMesh[
DiscretizeGraphics[
RegionDifference[MatrixRegion, InclusionRegion]]
, "RegionMarker" -> {{{0., 0., 0.}, 1, 10000}, {{xM/2, yM/2, zM/2},
2, 1000}}
, MaxCellMeasure -> {"Volume" -> 10000}
, "MeshOrder" -> 1];


I am grateful for any help or suggestions. Thanks in advances

Max

• Perhaps "MeshTools-1.0.0.paclet" is what you're looking for. Nov 16, 2020 at 11:58
• Nov 16, 2020 at 12:13

This answer extends @user21's to include different mesh densities of the inclusion along the X, Y, and Z directions.

It is important to note that the current mesher (version 12.1.1) likes to produce an isotropic mesh. One can accomplish the different mesh densities by creating a parameterized (I, J, K) structured mesh that ranges between zero and the number of elements in each direction. Then, one can rescale the coordinates from I, J, K space to user scaled coordinates.

First, let's create an isotropic structured mesh:

nx = 10; ny = 40; nz = 5;
isoMesh =
ToElementMesh[Cuboid[{0, 0, 0}, {nx, ny, nz}],
"MeshOrder" -> 1, MaxCellMeasure -> 1,
"RegionMarker" -> {{{nx, ny, nz}/2, 1}},
"MeshElementType" -> TetrahedronElement];
isoMesh["Wireframe"]


Second, let's create a rescaling transformation function from I, J, K space to user scaled coordinates:

scaledToUser =
RescalingTransform[{{0, nx}, {0, ny}, {0, nz}}, {{-xI/2,
xI/2}, {-yI/2, yI/2}, {-zI/2, zI/2}}];


Now, we can create the inner mesh by simply rescaling the coordinates like so:

innerMesh =
ToElementMesh[
"Coordinates" -> scaledToUser /@ isoMesh["Coordinates"],
"MeshElements" -> isoMesh["MeshElements"]];
innerMesh["Wireframe"]


Now, simply follow@user21's workflow with the new definition of innermesh and you can achieve different mesh densities along the X, Y, Z directions.

finalMesh[
"Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> (Directive[FaceForm[#],
EdgeForm[]] &  /@ {Orange, Blue}),
PlotRange -> {All, All, {-zM, zI/2}}]]
finalMesh[
"Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> (Directive[FaceForm[#],
EdgeForm[]] &  /@ {Orange, Blue}),
PlotRange -> {All, {0, yI/2}, {-zM, zI/2}}]]


# Structured hex mesh approach

As I alluded to in the comments, if you want to use a structured hex mesh for the inclusion, you probably want to propagate that through the entire mesh because the current version of Mathematica does not support pyramid and wedge type elements in 3D.

Depending on the nature of the physics you're trying to solve, there often can be sharp gradients at the interface regions. In this case, your solutions can often benefit by having a boundary layer mesh (or anisotropic mesh) where you have fine element layers at the interface that grow exponentially into the domain. These types of meshes can be quite economical in terms of element count.

## Workflow

### Helper functions

First, we will define some helper functions to create an anisotropic mesh.

(*Import required FEM package*)
Needs["NDSolveFEM"];
(* Define Some Helper Functions For Structured Quad Mesh*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=
Quiet@Abs@
FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 100000},
Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]
]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


### Construct a tensor product grid using RegionProduct tensor product mesh

Now, we can glue a bunch of segments that have different refinement strategies together in the horizontal, vertical, and depth directions as shown in the following workflow.

(*Define parameters*)
(*Lengths*)
h = 100;(*Horizontal*)
v = 10;(*Vertical*)
d = h;(*Depth*)
(*Number of elements per segment*)
nh = 10;
nv = 10;
nd = 10;
(*Association for Clearer Region Assignment*)
reg = <|"main" -> 1, "incl" -> 2|>;
(*Create mesh segments*)
(*Horizontal segments*)
(* left segment *)
(*First element is 1/50th of seg length*)
sh = rightSegmentGrowth[h, nh, h/50];
Print["Horizontal and depth segments"]
rh = pointsToMesh@(reflectRight@reflectRight[sh] - 2 h)
(*Vertical segment*)
sv = rightSegmentGrowth[v, nv, v/50];
Print["Vertical segment"]
rv = pointsToMesh@(reflectRight@reflectRight[sv] - 2 v)
(*Create tensor product grid with RegionProduct*)
rp = RegionProduct[rh, rv, rh];
(*Show the mesh*)
Print["Highlighted RegionProduct mesh"]
HighlightMesh[rp, Style[1, Orange]]


### Convert MeshRegion to ElementMesh with region markers

(*Extract Coords from RegionProduct*)
crd = MeshCoordinates[rp];
(*grab hexa element incidents RegionProduct mesh*)
inc = Delete[0] /@ MeshCells[rp, 3];
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {HexahedronElement[inc]}];
(*Extract bmesh*)
bmesh = ToBoundaryMesh[mesh];
(*Inclusion RegionMember Function*)
Ω3Dinclusion = Cuboid[{-h, -v, -h}, {h, v, h}];
rmf = RegionMember[Ω3Dinclusion];
regmarkerfn = If[rmf[#], reg["main"], reg["incl"]] &;
(*Get mean coordinate of each hexa for region marker assignment*)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@
ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
(*Create and view element mesh*)
Print["Converted Hexa Element Mesh Cutaway Drawing"]
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"],
"MeshElements" -> {HexahedronElement[inc, regmarkers]}];
mesh[
"Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> (Directive[Opacity[0.5], FaceForm[#](*,
EdgeForm[]*)] &  /@ {Blue, Orange}),
ViewPoint -> {-1.5, 0.8, -3}, ViewVertical -> {0, 1, 0},
PlotRange -> {{0, 2 h}, {0, 2 v}, {0, 2 h}}]]


Using a fully structured hex mesh, we created a fairly economical mesh (46656 hex elements) with very fine refinement at the interface.

• Top! That is exactly what I was looking for. Thank you.
– Max
Nov 17, 2020 at 14:25
• For comparison reasons I have to create the same structure, but instead of Tetrahedron elements with Hexahedron elements. That means a structured inclusion and a coarse unstructured matrix. I tried to apply the same procedure as you introduced for the Tets, but it does not work out so well. I can create the inclusion and rescale it. That is almost similar except of "MeshElementType" -> HexahedronElement. Then I want to transform it into bmesh1 = ToBoundaryMesh[ Cuboid[{-xM/2, -yM/2, -zM/2}, {xM/2, yM/2, zM/2}]] bmesh2 = ToBoundaryMesh[innerMesh]
– Max
Nov 23, 2020 at 16:33
• But the creation of the outerMesh does not really seem to work out that way.So I am curious if the generation of the Hexahedron mesh as I described above can be done analogously or not...
– Max
Nov 23, 2020 at 16:55
• @Max pure hex meshing can be quite difficult. There are open-source alternatives to produce a hex dominant mesh, but that requires the support for different element shapes such as pyramids and wedges (which the current version of Mathematica does not support). If your geometry is simple, your best bet to produce a pure hex mesh is to make the geometry completely structured. It still can be economical if you use anisotropic meshing as shown [here] (mathematica.stackexchange.com/a/232763/61809). Nov 23, 2020 at 17:10
• @Max I added a fully structured hex mesh approach in case you're still interested. Dec 4, 2020 at 4:38

There is an example of something similar in the PDEModel collection in the Acoustic Cloak Model. Here is a 3D version.

Some setup:

Needs["NDSolveFEM"]
xI = 200; yI = 200; zI = 20;
xM = xI*2; yM = yI*2; zM = zI*2;


We start by creating the inner mesh:

innerMesh =
ToElementMesh[Cuboid[{-xI/2, -yI/2, -zI/2}, {xI/2, yI/2, zI/2}],
"MeshOrder" -> 1, MaxCellMeasure -> 10000,
"RegionMarker" -> {{{0., 0., 0.}, 1}},
"MeshElementType" -> TetrahedronElement]

innerMesh["Wireframe"]


Check that the marker is there:

innerMesh["MeshElementMarkerUnion"]
{1}


Next we create a boundary mesh for the outer shape:

bmesh1 = ToBoundaryMesh[
Cuboid[{-xM/2, -yM/2, -zM/2}, {xM/2, yM/2, zM/2}]]


and extract the boundary mesh from the inner mesh:

bmesh2 = ToBoundaryMesh[innerMesh]


With the FEMAddOns you can combine them:

ResourceFunction["FEMAddOnsInstall"][]

Needs["FEMAddOns"]
bmesh = BoundaryElementMeshJoin[bmesh1, bmesh2]

bmesh["Wireframe"]


Now comes the key point. When we generate the full outer mesh we make sure that no new nodes are inserted on the boundary. That is done with setting "SteinerPoints"->False.

outerMesh = ToElementMesh[bmesh,
"SteinerPoints" -> False,
"RegionHoles" -> {{0, 0, 0}},
"RegionMarker" -> {{{xM/2, yM/2, zM/2}, 2, 1000}},
MaxCellMeasure -> {"Volume" -> 10000}, "MeshOrder" -> 1]


Now, that we have an inner and outer mesh that align at the inner material region we can make the final full mesh:

innerCoordinates = innerMesh["Coordinates"];
outerCoordinates = outerMesh["Coordinates"];
finalMesh =
ToElementMesh[
"Coordinates" -> Join[outerCoordinates, innerCoordinates],
"MeshElements" ->
Flatten[{outerMesh["MeshElements"],
TetrahedronElement, {ElementIncidents[
innerMesh["MeshElements"]] + Length[outerCoordinates],
ElementMarkers[innerMesh["MeshElements"]]}]}]]


Check that the markers are there:

finalMesh["MeshElementMarkerUnion"]
{1,2}


And Visualize:

finalMesh[
"Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> (Directive[FaceForm[#],
EdgeForm[]] & /@ {Orange, Blue}),
PlotRange -> {All, All, {-zM, zI/2}}]]


If you want to generate a second order mesh you can do so by

MeshOrderAlteration[finalMesh, 2]

• Thank's so fare. Is there also a way to define the mesh density of the inclusion via x-, y- and z- directions. So that you define hard the numbers of elements in each direction e.g with {nx,ny,nz} and not via the maximal cell measure?
– Max
Nov 16, 2020 at 16:40
• It is nice to have a 3D example of BoundaryElementMeshJoin`. I extended your answer to include different mesh densities along the X, Y, and Z directions. Let me know if you want to include it in your answer, and I will delete mine since you did all the heavy lifting. Nov 16, 2020 at 23:25
• @Tim Laska, sure modify as you like. Nov 17, 2020 at 5:42