# Trouble meshing a Corbino disc

I am trying to simulate current flow using FEM through this particular device (a Corbino disk) and I have difficulties generating the mesh for it. Here is the region:

Needs["NDSolveFEM"]

Ω = ImplicitRegion[
0   <= x^2 + y^2 <= 8^2 && 1>= z>= 0||
0^2 <= x^2 + y^2 <= 2^2 && 3>= z>= 1||
3^2 <= x^2 + y^2 <= 4^2 && 3>= z>= 1||
5^2 <= x^2 + y^2 <= 6^2 && 3>= z>= 1||
7^2 <= x^2 + y^2 <= 8^2 && 3>= z>= 1,
{{x, -8, 8}, {y, -8, 8}, {z, 0, 3}}
];

Show[RegionPlot3D[Ω, PlotPoints -> 100], AspectRatio -> Automatic] When I try to mesh I am left with this:

    bmesh = ToElementMesh[Ω];
bmesh["Wireframe"]


Throw: Uncaught Throw[\$Failed,RegionMeshRegionException[]] returned to top level.

I am able to mesh a quarter of the device and a half of the device with no problems. A quarter of the device and the mesh are shown below:

 Ω = ImplicitRegion[
0   <= x^2 + y^2 <= 8^2 && x - y >= 0 && x >= 0 && y >= 0 && 1 >= z >= 0 ||
0^2 <= x^2 + y^2 <= 2^2 && x - y >= 0 && x >= 0 && y >= 0 && 3 >= z >= 1 ||
3^2 <= x^2 + y^2 <= 4^2 && x - y >= 0 && x >= 0 && y >= 0 && 3 >= z >= 1 ||
5^2 <= x^2 + y^2 <= 6^2 && x - y >= 0 && x >= 0 && y >= 0 && 3 >= z >= 1 ||
7^2 <= x^2 + y^2 <= 8^2 && x - y >= 0 && x >= 0 && y >= 0 && 3 >= z >= 1,
{{x, 0, 8}, {y, 0, 8}, {z, 0, 3}}];

Show[RegionPlot3D[Ω], AspectRatio -> Automatic]

bmesh = ToElementMesh[Ω];
bmesh["Wireframe"]  I have also tried building the device piece by piece and then simply using RegionUnion, but that also does not seem to do the trick.

I now suspect that there may be some sort of problem with intersecting boundaries, but I have had that issue before and it always returns an error from ElementMesh, and not from Throw.

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["NDSolveFEM"]
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. 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[]]]] This mesh has 610K tets. You'd probably need to find something in between.

• I'm able to run the first mesh but when I try and improve with a better boundary discretization I'm given the error: ElementMesh:The input has or generated an intersecting boundary and cannot be processed. I've just copied and pasted your code, and am not sure what I did wrong. Any ideas? Jun 27 '18 at 8:44
• @A.May what version of M- are you using? Jun 27 '18 at 9:03
• Version 11.0.0.0 Jun 27 '18 at 9:14
• It works fine for me in Version 11.3. Can you update? Jun 27 '18 at 9:20
• Will do, and I'll inform you once I've run it there. Jun 27 '18 at 10:24

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["NDSolveFEM"]
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]
},
], 2],
Join[radii[[-1]] circle, ConstantArray[h[], {n, 1}], 2],
Join[radii[[-1]] circle, ConstantArray[h[], {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],
];

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: 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}]]
}] 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]]] You get quite pretty 3D unstructured mesh of hexahedron elements.

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