# MeshRefinementFunction missing points

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["NDSolveFEM"]
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.

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Feb 5 '16 at 16:47
• It seems that there is a problem with your code. I have tried It, but after I launch the cell myMesh=ToElementMesh..., Mathematica is hanging. No error/warning message. I have aborded the calculations after 15 minutes. (my machine : I7, Mathematica 10.0) – andre314 Feb 5 '16 at 17:29
• I am running 10.3.0.0, so maybe it's a version issue. Can you get it to work without the mesh refinement? – Commander Ellen Tigh Feb 5 '16 at 17:37
• It doesn't work even without the mesh refinement (the modified code is myMesh = ToElementMesh[DiscretizeRegion@myRegion]; myMesh["Wireframe"["MeshElementStyle" -> FaceForm[White]]]). I must leave. be back in 1 hour. – andre314 Feb 5 '16 at 17:49
• I have tried your code on Mathematica 10.3.1 (on the Web, on "Wolfram Developement Platform") without the mesh refinement and it works. If I add the mesh refinement, the result is exactly the same. This new problem is due to an error in your code (this time). The correct syntax is myMesh = ToElementMesh[DiscretizeRegion[myRegion, MeshRefinementFunction -> meshRefine], "ImproveBoundaryPosition" -> True]; myMesh["Wireframe"["MeshElementStyle" -> FaceForm[White]]]. Now I see the problem that you describes wich is the object of you question. – andre314 Feb 5 '16 at 19:48

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.