MeshRefinementFunction missing points

Question

I'm doing some FEM on a rectangular region with a spherical cavity:

myRegion = ImplicitRegion[X^2 + Y^2 + (Z - 1.05)^2 >= 1, {{X, 0, 1.5},{Y,0,1.5}, {Z, 0, 1.05}}];
RegionPlot3D[myRegion, Boxed -> False]

The most important part of the region is that in the bottom right (where the cavity is nearly reaching the bottom), so I'd like to selectively refine in that region. I do this with a function that refines when the mesh points are close to the z-axis:

Needs["NDSolve`FEM`"]
meshRefine = Function[{vertices, volume}, Block[{x, y, z}, {x, y, z} = Mean[vertices]; If[Sqrt[x^2 + y^2] <= 0.5, volume > 0.000005, False]]];

Now for the mesh. I have found that "ImproveBoundaryPosition" is necessary for applying boundary conditions to NDSolve, and that applying DiscretizeRegion before ToElementMesh avoids a "femimq: The element mesh has insufficient quality..." error.

myMesh=ToElementMesh[DiscretizeRegion@myRegion, "ImproveBoundaryPosition" -> True,MeshRefinementFunction -> meshRefine];
myMesh["Wireframe"["MeshElementStyle" -> FaceForm[White]]]

These two images show the result with and without the refinement function (from below):

---

As you can see, the MeshRefinementFunction gets some (but not all) of the points near the thinnest part. It seems that when the initial mesh connects the spherical boundary to the bottom, MeshRefinementFunction does not apply. Any ideas for a better way to refine that region? Thanks in advance.

Answer 1

I suspect that the cells closest to the $z$ axis in your region did not get optimized because their volume was simply too small to trigger the mesh refinement function: the solid is so thin in the vicinity of the $z$ axis, and its thickness is changing so rapidly, that any constant volume requirement small enough not to trigger the generation of thousands of useless cells in the thicker part of your refinement section ends up being too large to trigger refinement of the very thin slice close to the $z$ axis.

You may be able to improve the situation using an "adaptive" refinement function, whose volume cutoff threshold is smaller and smaller the closer you come to the $z$ axis.

Considering that $\sqrt{x^2+y^2}<1$ in your region of interest, one could use:

meshRefine = Function[{vertices, volume},
   Block[{x, y, z},
    {x, y, z} = Mean[vertices];
    volume > 5*^-5 Sqrt@(x^2 + y^2)
    ]
   ];

This generates the following mesh:

DiscretizeRegion[myRegion, MeshRefinementFunction -> meshRefine]

Note that I chose the mesh refinement function's $5\times10^{-5}$ threshold entirely arbitrarily to generate a number of cells that looked appropriate to me; of course you should adjust it for your own purposes.

Answer 2

NOT A ANSWER

Once corrected, your code gives the following with Mathematica 10.3.1 (on the "Wolfram Development PlatForm")

It looks like a little bit better than what you obtain.

EDIT

I can't investigate further : my code crashes on the "Wolfram Development Platform"