Create graded mesh

Question

I want to create graded mesh inside a cube, with finer mesh in the left half and a coarser mesh in the right half. My plan was to partition the cube with a 2D plane and then generate boundary mesh with ToBoundaryMesh (it was intended to create a boundary mesh on the partition as well) and then apply ToElementMesh with different MaxCellMeasure on the left and right of the partition. Here is my code:

reg = RegionUnion[Cuboid[], 
  ImplicitRegion[x == 0.5 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}]];
ToBoundaryMesh[reg]

which gives the following error:

BoundaryDiscretizeRegion::brepl: There are components in RegionUnion[Cuboid[{0,0,0}],ImplicitRegion[x==0.5&&0<=y<=1&&0<=z<=1,{x,y,z}]] having dimension lower than the embedding dimension 3 that will not be included in the boundary representation.

How can I achieve a graded mesh? Any help is much appreciated.

Edit

I have tried the following. I created two cuboids with a common boundary and joined them with Or. Then I created a mesh while marking the two cuboid regions (labels 10 and 20) using RegionMarker, in which I also specified the MaxCellMeasure of the two regions to be 1 and 0.01. Here is the code:

reg = ImplicitRegion[Or[0 <= x <= 0.5 && 0 <= y <= 1 && 0 <= z <= 1,
0.5 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1], {x, y, z}];
mesh = ToElementMesh[reg, 
   "RegionMarker" -> {{{0.1, 0.5, 0.5}, 10, 1}, {{0.6, 0.5, 0.5}, 20, 0.01}}];

But it still creates a mesh of uniform size in the two regions. Any ideas on how to overcome the problem?

Answer 1

Here is a method to build an anisotropic hexahedron mesh, although it does not meet the OP's desire for a "simple" solution.

The inspiration comes from the Tensor Product Grid Example. The idea is to create a series of 1D mesh segments along each direction and build the mesh using RegionProduct.

Here is an example of a graded mesh connected to a uniform mesh along the x direction. A practical example of this type of meshing would be conjugate heat transfer where the fluid has a boundary layer or graded mesh to capture thermal gradients at the wall and the solid can be meshed uniformly.

Some helper functions to build graded meshes:

(* Import required 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];
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

A commented workflow to create a graded mesh:

(*Define parameters*)
(*Lengths*)
(*2 horizontal segments*)
h1 = 0.5;
h2 = 0.5;
v = 1;(*Vertical*)
d = 1;(*Depth*)
(*Number of elements per segment*)
nh1 = 15;
nh2 = 10;
nv = 20;
nd = 5;
(*Association for Clearer Region Assignment*)
reg = <|"left" -> 1, "right" -> 2|>;
(*Create mesh segments*)
(*Horizontal segments*)
(* left segment *)
(*First element is 1/50th of seg length*)
(*Flip segment so smallest elm at interface *)
sh1 = flipSegment@meshGrowthByElm0[h1, nh1, h1/50];
(*Make right segment uniform mesh size*)
sh2 = Subdivide[h2, nh2];
(*Glue segments together*)
segh = extendMesh[sh1, sh2];
(*View individual horizontal segments*)
Print["Horizontal segments"]
pointsToMesh /@ {sh1, sh2}
(*View combined segments*)
Print["Combined horizontal segments"]
rh = pointsToMesh@segh
(*Vertical Segment*)
Print["Vertical segment"]
rv = pointsToMesh@Subdivide[v, nv]
(*View Region Product of horiz and vert segs*)
Print["2D Region via RegionProduct"]
RegionProduct[rh, rv]
(*Depth Segment*)
Print["Depth segment"]
rd = pointsToMesh@Subdivide[d, nd]
(*Create a tensor product grid from h,v,and d segments*)
rp = RegionProduct[rh, rv, rd];
(*View mesh*)
Print["Full Extruded 3D Region via RegionProduct"]
HighlightMesh[rp, Style[1, Orange]]
(*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];
(*Iron RegionMember Function*)
Ω3Diron = Cuboid[{0, 0, 0}, {h1, v, d}];
rmf = RegionMember[Ω3Diron];
regmarkerfn = If[rmf[#], reg["left"], reg["right"]] &;
(*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"]
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"], 
   "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
Graphics3D[
 ElementMeshToGraphicsComplex[bmesh, 
  VertexColors -> (ColorData["BrightBands"] /@ 
     Rescale[regmarkerfn /@ bmesh["Coordinates"]])], Boxed -> False]