Finer mesh for selected subregion of a solid 3D cylinder

Question

I am trying to create a finer mesh in a subregion of a solid cylinder in order to improve the resolution there. I don't want to have too many elements for the entire cylinder because it slows down the computation. From other questions here and the tutorial on Element Mesh Creation https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html, I see that this can be done for 2D shapes and composite shapes (different materials).

1. I'd like a finer mesh along the axis of the cylinder, with a smaller radius (e.g. 1/3 of the radius).

I suppose I can break it in two separate cylinders--the inner solid cylinder and the outer "pipe" and then combine them with RegionUnion. But, then I am not sure how to define the boundary surfaces for the boundary conditions (outer cylindrical surface and the whole end faces). I hope there is a more elegant solution.

2. How to create finer mesh just for the end faces? The boundary surfaces are defined with ElementMarker. I tried defining separate boundary mesh for each face and then turning it ToElementMesh, but it didn't work. Below is the code for the whole solid cylinder.

OpenCascadeShapeType[cyl]; 
cyl = OpenCascadeShape[c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003]];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[cyl]
mesh = ToElementMesh[bmesh]
Show[mesh["Wireframe"], Graphics3D[c1]]

When I visualize the solution with SliceContourPlot3D, the resolution is insufficient to capture small features near the axis.

Answer 1

Version 1

One way to do it could be to generate a 2D disk mesh with an appropriate refinement and the extrude that to 3D.

(m2 = ToElementMesh[Disk[{0, 0}, 0.003], 
    MeshRefinementFunction -> 
     Function[{vertices, area}, 
      area > 0.00000005 (0.001 + 
          500 Norm[Mean[vertices]])]])["Wireframe"]

mesh = ElementMeshRegionProduct[m2, 
  ToElementMesh[Line[{{0}, {0.01}}]]]


mesh["Wireframe"]

Here is a view inside:

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{-0.003, 0}, All, {0, 0.005}}]]

Version 2

A variant of this approach is to use a graded mesh for the element mesh region product, like so:

m1 = ToGradedMesh[{Line[{{0}, {0.01}}], <|
     "Alignment" -> "BothEnds"|>}];
MeshRegion[m1]

mesh = ElementMeshRegionProduct[m2, m1];

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{-0.003, 0}, All, {0, 0.005}}]]

Note, how the mesh is becoming denser towards the bottom (and top not shown) Here is a view of the boundaries:

mesh["Wireframe"[
  "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red], 
    FaceForm[Blue]}, PlotRange -> {{-0.003, 0}, All, {0, 0.01}}]]

Version 3

Here is a version that has a refined boundary but a coarse interior:

Needs["NDSolve`FEM`"]
c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
mesh = ToElementMesh[c1, "MaxBoundaryCellMeasure" -> 0.1, 
   "MaxCellMeasure" -> Infinity];
Show[mesh["Wireframe"], Graphics3D[c1]]

Cross section:

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{0, 0.005}, {0, 0.003}, {0, 0.003}}]]

Version 4

Yet another way to do it would be by using OpenCascade and create an inner cylinder that is refined to a different level.

Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]
c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
c2 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003/2];

s1 = OpenCascadeShape[c1];
s2 = OpenCascadeShape[c2];

faces = Flatten[OpenCascadeShapeFaces /@ {s1, s2}];
union = OpenCascadeShapeUnion[faces];

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];

bmesh["Wireframe"]

Create to material markers p1 and p2:

p1 = {0.01/2, 0, 0};
p2 = {0.01/2, 0, 0.005/2};
Show[
 Graphics3D[{{Red, PointSize[0.02], Point[p1]}, {Green, 
    PointSize[0.02], Point[p2]}}],
 bmesh["Wireframe"],
 Graphics3D[{Opacity[0.2], c1}]]

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{p1, 1, 10^-12}, {p2, 2, 10^-10}}];
Show[mesh["Wireframe"], Graphics3D[c1]]

Version 5

Here is a version with a MeshRefinementFunction. By testing I found that the mesh changes when a "MaxCellMeasure" around 10^-11 was given. So I create a function that re-scales around this value from 0 to the boundary at 0.003.

Needs["NDSolve`FEM`"]
rescale = Rescale[Sqrt[y^2 + z^2], {0, 0.003}, {10^-14, 10^-11}]

The Sqrt[y^2 + z^2] is just the radial distance from the x-axis. This I then injected into the mesh refinement function:

c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
cf = With[{r = rescale},
   Compile[{{coordinates, _Real, 2}, {vol, _Real, 0}},
    Block[{com, x, y, z, rvol},
     com = Mean[coordinates];
     {x, y, z} = com;
     rvol = r;
     If[vol > rvol, True, False]
     ]
    ]
   ];
mesh = ToElementMesh[c1, MeshRefinementFunction -> cf]
Show[mesh["Wireframe"], Graphics3D[c1]]