3D Mesh refinement near sphere boundary and on surface revolutions

Question

I am trying to solve an advection-diffusion problem inside surface of revolutions. I have high velocity gradients near the boundary of the bounding surface. Is there a way to refine the mesh near the boundary? I have gone through the documentation of the MeshRefinementFunction but I am unable to figure it out. In particular for general surface of revolution I am using opencascade and it is not clear how to do mesh refinement on the mesh generated by open cascade. Here are two sample codes (one for the sphere and one for a egg shape surface of revolution).

ToElementMesh[Sphere[], "MeshOrder" -> 2, 
  "MaxBoundaryCellMeasure" -> 0.02, 
  MeshRefinementFunction -> 
   Function[{vertices, Vol}, 
    Block[{x, y, z}, {x, y, z} = Mean[vertices];
     If[x^2 + y^2 + z^2 > 0.8, Vol > 0.01, Vol > 0.01]]]];

I will like the mesh in a spherical shell near the boundary to be refined.

ClearAll["Global`*"];
Needs["NDSolve`FEM`"];
Needs["OpenCascadeLink`"];
(*Riemeann map of the unit disc*)
aY = 0.15; aZ = 0.1; area = Pi;
aW = Sqrt[area/Pi + aY^2 + aZ^2];
sigma = Exp[I t]; scale = 5;
P = aW*sigma + aY/sigma + aZ/(Sqrt[2]*sigma^2);
zt = ComplexExpand[Re[P]];
rt = ComplexExpand[Im[P]];

(* Define parametric region *)
\[CapitalOmega] = 
  ParametricRegion[{R rt Cos[\[Theta]], R rt Sin[\[Theta]], 
    R zt}, {{t, 0, \[Pi]}, {\[Theta], 0, 2 \[Pi]}, {R, 0, scale}}];

(* Generate mesh for velocity evaluation *)
npoly = 100;
pts = Table[{0, scale*rt, scale*zt}, {t, 0, Pi, Pi/npoly}];
poly = Line[pts]; 

(* OpenCascade meshing *)
wire = OpenCascadeShape[poly];
axis = {{0, 0, 0}, {0, 0, 1}};
sweep = OpenCascadeShapeRotationalSweep[wire, axis, 2 \[Pi]];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep];
mesh = ToElementMesh[bmesh];
Dimensions[mesh["Coordinates"]]
Show[Graphics3D[{{Red, Thickness -> 0.01, poly}, {Blue, Thick, 
    Arrow[axis]}}], mesh["Wireframe"], Boxed -> False]

Similarly here the mesh needs to be refined around a shell of similar shape near the boundary. Thank you!

Answer 1

With the mesh refinement you can not quite get to the boundary; you'd need an interplay between AccuracyGoal (to refine the actual boundary) and the MeshRefinementFunction. Something like this:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Ball[],
   AccuracyGoal -> 3.1,
   MeshRefinementFunction -> 
    Function[{vertices, Vol}, 
     Block[{x, y, z}, {x, y, z} = Mean[vertices];
      If[x^2 + y^2 + z^2 > 0.8, Vol > 0.00001, Vol > 0.1]]]];
mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "MeshElementStyle" -> {FaceForm[Orange]}, 
  PlotRange -> {{0, 1.1}, {0, 1.1}, {0, 1.1}}]]