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?

• Perhaps my answer here 231273 will give you some ideas how to create a mapped mesh. Oct 12 '20 at 4:10
• @TimLaska Thanks. You have indeed obtained a graded mesh in your answer. The code is too complicated for me to follow. Is it possible to show how it can be simply done for the cube in my question?
– Deep
Oct 12 '20 at 4:23
• In my opinion, the 3D unstructured meshing (being more complicated), is not as developed as the 2D meshing as of v12.1.1. So, some functionality may not yet be robust. I can show an example of a graded 3D mesh, but I am operating at a lower level and essentially building the mesh by hand. Therefore, it probably won't meet your simple criterion. Also note that the MaxCellMeasure in your RegionMarker spec is on a volume basis and not length. I don't think it will help in this case (at least for me). Oct 12 '20 at 15:05
– Deep
Oct 13 '20 at 3:46

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["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];
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 /@ 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] • Awesome. Thanks.
– Deep
Oct 15 '20 at 2:49
• @Deep, Tim made a proposal to add anisotropic meshing capabilities to the FEM mesh generation here. If you think this is worthwhile, consider giving this proposal an upvote Feb 22 at 10:36
• @user21 Done. I believe it is a useful feature.
– Deep
Feb 24 at 9:42