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:
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} *)
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.
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.