Trouble meshing a Corbino disc

I am a friend of constructing meshes by hand. This might not be as convenient as using automated discretization tools but the outcome is a bit more controllable. Here is is a code for that:

Needs["NDSolve`FEM`"]
n = 120;
radii = N@{2, 3, 4, 5, 6, 7, 8};
h = N@{0, 1, 3};
circle = CirclePoints[{1., 0.}, n];
p = Join[
   Flatten[
    Table[
     {
      Join[radii[[i]] circle, ConstantArray[h[[3 - UnitStep[(-1)^i]]], {n, 1}], 2],
      Join[radii[[i]] circle, ConstantArray[h[[2 + UnitStep[(-1)^i]]], {n, 1}], 2]
      },
     {i, 1, Length[radii] - 1}
     ], 2],
   Join[radii[[-1]] circle, ConstantArray[h[[3]], {n, 1}], 2],
   Join[radii[[-1]] circle, ConstantArray[h[[1]], {n, 1}], 2]
   ];
polys = Join[
   {Range[1, n]},
   Flatten[Table[
     Join[
      Partition[Range[i n + 1, (i + 1) n], 2, 1, 1],
      Reverse /@ Partition[Range[(i + 1) n + 1, (i + 2) n], 2, 1, 1],
      2
      ],
     {i, 0, 2 Length[radii] - 2}], 1],
   {Range[(2 Length[radii] - 1) n + 1, 2 n Length[radii]]}
   ];

R = BoundaryMeshRegion[p, Polygon[polys]]; // AbsoluteTiming // First
S = ToElementMesh[R] // AbsoluteTiming

1.53396

{1.49476, ElementMesh[{{-8., 8.}, {-8., 8.}, {0., 3.}}, {TetrahedronElement[ "<" 34739 ">"]}]}

This is how R looks like:

enter image description here

You can make it "rounder" by increasing n in the code above.


I will show another approach for creating hexahedral mesh on OP's domain that allows some more control over element size and distribution. It involves using open source package MeshTools (I am developing it.).

One should first install the package and then load it.

Needs["MeshTools`"]

Create 2D region of Corbino disc top view and convert it to boundary mesh. bmesh here contains just discretized circles.

reg = RegionUnion[
   Disk[{0, 0}, 2],
   Annulus[{0, 0}, {3, 4}],
   Annulus[{0, 0}, {5, 6}],
   Annulus[{0, 0}, {7, 8}]
 ];
bmesh = ToBoundaryMesh[reg]
(* ElementMesh[{{-8., 8.}, {-8., 8.}}, Automatic] *)

Make 2D mesh of triangles over boundary mesh bmesh and assign different markers (integer numbers) to parts of Corbino disk that will have different height.

meshTri = ToElementMesh[bmesh,
  "MeshOrder" -> 1,
  "RegionHoles" -> None,
  "RegionMarker" -> {{{0,0},1},{{2.5,0},2},{{3.5,0},1},{{4.5,0},2},{{5.5,0},1},{{6.5,0},2},{{7.5,0},1}},
   MeshQualityGoal -> 1
]
(* ElementMesh[{{-8., 8.}, {-8., 8.}}, {TriangleElement["<" 570 ">"]}] *)

Use TriangleToQuadMesh and SmoothenMesh functions to convert triangular mesh to nice quadrilateral mesh and visualize both meshes. We get the best quality if we smoothen triangular mesh, convert it to quadrilaterals and smoothen it again.

meshQuad = SmoothenMesh@TriangleToQuadMesh@SmoothenMesh[meshTri]
(* ElementMesh[{{-8., 8.}, {-8., 8.}}, {QuadElement["<" 1200 ">"]}] *)

Row[{
  meshTri["Wireframe"["MeshElementStyle" -> FaceForm /@ {Red, Green}]],
  meshQuad["Wireframe"["MeshElementStyle" -> FaceForm /@ {Red, Green}]]
 }]

mesh2D

Then we split the meshQuad by element marker to two different meshes, extrude them for different thickness (number of element layers should be chosen accordingly) and merge them together in one ElementMesh.

meshHigh = SelectElements[meshQuad, ElementMarker == 1]
meshLow = SelectElements[meshQuad, ElementMarker == 2]

mesh = MergeMesh[{
   ExtrudeMesh[meshHigh, 3, 6],
   ExtrudeMesh[meshLow, 1, 2]
 }]
(* ElementMesh[{{-8., 8.}, {-8., 8.}, {0., 3.}}, {HexahedronElement["<" 5104">"]}] *)

mesh["Wireframe"["MeshElementStyle" -> FaceForm[LightBlue]]]

mesh3D

You get quite pretty 3D unstructured mesh of hexahedron elements.

Through[{Min, Mean, Max}@Flatten[mesh["Quality"]]]
(* {0.363397, 0.855278, 0.999059} *)

Here is a way to do it: First we create a symbolic region by

w = Fold[RegionDifference, Cylinder[{{0, 0, 0}, {0, 0, 3}}, 8], 
   RegionDifference[Cylinder[{{0, 0, 1}, {0, 0, 4}}, #], 
      Cylinder[{{0, 0, 0}, {0, 0, 4}}, # - 1]] & /@ {7, 5, 3}];

When we mesh that:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[w]
mesh["Wireframe"[
  "MeshElementStyle" -> Directive[FaceForm[LightBlue], EdgeForm[]]]]

we get a mesh with 26K tets and some internal edges are not as sharp as they could be.

enter image description here

We can improve that with a better boundary discretization like:

mesh = ToElementMesh[\[CapitalOmega], 
  "BoundaryMeshGenerator" -> {"BoundaryDiscretizeRegion", 
    Method -> {"MarchingCubes", PlotPoints -> 33}}]
mesh["Wireframe"[
  "MeshElementStyle" -> Directive[FaceForm[LightBlue], EdgeForm[]]]]

enter image description here

This mesh has 610K tets. You'd probably need to find something in between.