16
$\begingroup$

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

enter image description here

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:

enter image description here

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

Improve this question
$\endgroup$
2

2 Answers 2

Reset to default
12
$\begingroup$

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

Structured mesh

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

User scaled mesh

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

Mesh with inclusion

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["NDSolve`FEM`"];
(* 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]]

Region product mesh

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

Structured hex mesh

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

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Top! That is exactly what I was looking for. Thank you. $\endgroup$
    – Max
    Nov 17, 2020 at 14:25
  • $\begingroup$ 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] $\endgroup$
    – Max
    Nov 23, 2020 at 16:33
  • $\begingroup$ 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... $\endgroup$
    – Max
    Nov 23, 2020 at 16:55
  • $\begingroup$ @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). $\endgroup$
    – Tim Laska
    Nov 23, 2020 at 17:10
  • $\begingroup$ @Max I added a fully structured hex mesh approach in case you're still interested. $\endgroup$
    – Tim Laska
    Dec 4, 2020 at 4:38
15
$\begingroup$

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

Some setup:

Needs["NDSolve`FEM`"]
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"]

enter image description here

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

enter image description here

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"], 
     MapThread[
      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}}]]

enter image description here

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

MeshOrderAlteration[finalMesh, 2]
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – Max
    Nov 16, 2020 at 16:40
  • $\begingroup$ 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. $\endgroup$
    – Tim Laska
    Nov 16, 2020 at 23:25
  • $\begingroup$ @Tim Laska, sure modify as you like. $\endgroup$
    – user21
    Nov 17, 2020 at 5:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.