5
$\begingroup$

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!

$\endgroup$

1 Answer 1

6
$\begingroup$

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

enter image description here

$\endgroup$
2
  • $\begingroup$ Wow! Thanks a lot. Can something similar be done for axisymmetric objects with parameteric representation of {r(t),z(t)} that is rotated say around the z axis? $\endgroup$
    – bchakra
    Jul 8, 2022 at 13:18
  • $\begingroup$ If you mean the second case in your example, then yes, that can be done in the same way. $\endgroup$
    – user21
    Jul 8, 2022 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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