Meshing a thin tube

Question

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]

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?

Answer 1

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"]

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"]

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

Histogram[mesh["Quality"]]