9
$\begingroup$

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["NDSolve`FEM`"]

Ω = 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]

enter image description here

When I try to mesh I am left with this:

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

Throw: Uncaught Throw[$Failed,Region`Mesh`RegionException[]] 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"]

enter image description here

enter image description here

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.

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – A. May Jun 27 '18 at 8:44
  • $\begingroup$ @A.May what version of M- are you using? $\endgroup$ – user21 Jun 27 '18 at 9:03
  • $\begingroup$ Version 11.0.0.0 $\endgroup$ – A. May Jun 27 '18 at 9:14
  • $\begingroup$ It works fine for me in Version 11.3. Can you update? $\endgroup$ – user21 Jun 27 '18 at 9:20
  • $\begingroup$ Will do, and I'll inform you once I've run it there. $\endgroup$ – A. May Jun 27 '18 at 10:24
7
$\begingroup$

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.

$\endgroup$
5
+100
$\begingroup$

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 packages, FEMAddOns by WRI and MeshTools with my set of utilities.

One should first install the packages and load them.

<<FEMAddOns`
<<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}}
]
(* ElementMesh[{{-8., 8.}, {-8., 8.}}, {TriangleElement["<" 522 ">"]}] *)

Use ToQuadMesh and ElementMeshSmoothing functions from FEMAddOns package to convert triangular mesh to nice quadrilateral mesh and visualize both meshes.

meshQuad = ElementMeshSmoothing@ToQuadMesh[meshTri]
(* ElementMesh[{{-8., 8.}, {-8., 8.}}, {QuadElement["<" 1090 ">"]}] *)

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

mesh2D

Using functions from MeshTools package, 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 = SelectElementsByMarker[meshQuad, 1];
meshLow = SelectElementsByMarker[meshQuad, 2];

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

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

mesh3D

You get quite pretty 3D unstructured mesh of hexahedron elements. Then you can convert it to tetrahedron elements with ToTetrahedronMesh or increase mesh order with MeshOrderAlteration[mesh,2]. All these methods give you plenty of possibilities to test in your FEM analysis and compare their performance.

ToTetrahedronMesh[mesh]
(* ElementMesh[{{-8., 8.}, {-8., 8.}, {0., 3.}}, {TetrahedronElement["<" 22860 ">"]}] *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.