7
$\begingroup$

I'd like to mesh a thin tube. The following gives a strange result when thickness becomes too small:

height = 3;
radius = 1;
thickness = .06;

ir = ImplicitRegion[(radius - thickness)^2 <= x^2 + y^2 <= radius && 
    0 <= z <= height, {x, y, z}];
mesh = ToElementMesh[
  ir, {{-radius, radius}, {-radius, radius}, {0, height}}, 
  "MeshOrder" -> 1, MaxCellMeasure -> 0.0003]

enter image description here

As you can see it is not quite a cylinder...

Ideally, I'd like to mesh cylinders with thicknesses as small as 0.005. How can it be done?

$\endgroup$
9
$\begingroup$

This can sometimes happen if the region bounds intersect with the region it self. If you look at the output of the ElementMesh you will see that it's bounds do not go from -radius to radius.

mesh["Bounds"]
{{-0.995306, 0.995306}, {-0.995306, 0.995306}, {0., 3.}}

The fix is easy: just enlarge the bounding box a bit:

height = 3;
radius = 1;
thickness = .06;
rt = radius + 0.1;

Needs["NDSolve`FEM`"]
ir = ImplicitRegion[(radius - thickness)^2 <= x^2 + y^2 <= radius && 
    0 <= z <= height, {x, y, z}];
mesh = ToElementMesh[ir, {{-rt, rt}, {-rt, rt}, {0, height}}, 
  "MeshOrder" -> 1, MaxCellMeasure -> 0.0003];

Look at the bounds:

mesh["Bounds"]
{{-1., 1.}, {-1., 1.}, {0., 3.}}

Look at the mesh:

mesh["Wireframe"]

enter image description here

Update:

Another option is to manually generate the a hex element based mesh for this case:

Needs["NDSolve`FEM`"]

nx = 100; ny = 5; nz = 100;
coordinates = 
  Flatten[ Table[{r Cos[\[Theta]], r Sin[\[Theta]], h}, {h, 0., 
     3., (3 - 0)/(nz - 1)}, {r, 1. - 0.05, 
     1., (1. - (1 - 0.05))/(ny - 1)}, {\[Theta], 0., 
     2 Pi, (2 Pi - 0.)/(nx - 1)}], 2];

mkIncidents = 
  Compile[{{nx, _Integer, 0}, {ny, _Integer, 0}, {nz, _Integer, 0}},
   Flatten[
    Table[Block[{p1 = (j - 1)*nx + i, p2 = j*nx + i, p3 = p2 + 1, 
       p4 = p1 + 1, p5, p6, p7, p8},
      {p5, p6, p7, p8} = {p1, p2, p3, p4} + k*nx*ny;
      {p1, p2, p3, p4} += (k - 1)*nx*ny;
      {p1, p2, p3, p4, p5, p6, p7, p8}], {i, 1, nx - 1}, {j, 1, 
      ny - 1}, {k, 1, nz - 1}], 2]
   ];

incidents = mkIncidents[nx, ny, nz];

mesh = 
 ToElementMesh["Coordinates" -> coordinates, 
  "MeshElements" -> {HexahedronElement[incidents]}]

ElementMesh[{{-0.999497, 1.}, {-0.999874, 0.999874}, {0., 
   3.}}, {HexahedronElement["<" 39204 ">"]}]

mesh["Wireframe"]

enter image description here

Playing with the nx,ny and nz allows you to influence the the quality which you can inspect with:

Histogram[mesh["Quality"]]
$\endgroup$
  • $\begingroup$ Do you manage to mesh with thickness = 0.01? I either get an error, or I can reduce MaxCellMeasure but it never ends. $\endgroup$ – anderstood Feb 2 '18 at 18:03
  • $\begingroup$ @anderstood, yes with MaxCellMeasure -> 0.0000015 it gives me 1.3M tets. I am a bit surprised that I need to specify the MaxCellMeasure, but I assume because of the high aspect ratio the automatic MaxCellMeasure computation fails. This does not look quite correct. $\endgroup$ – user21 Feb 2 '18 at 19:16
  • $\begingroup$ OK, thanks. I stopped the computation after 30 min... Anyway I think the mesh is not quite optimal for the FEM (elements are pretty distorted). I need to read more about meshing with Mathematica. Thanks again for your help. $\endgroup$ – anderstood Feb 2 '18 at 19:25
  • 1
    $\begingroup$ @anderstood with V11.2 on my laptop it took 70sec. so 30 Min is way too long. Another idea: See update. $\endgroup$ – user21 Feb 2 '18 at 20:02

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.