# Problem when using the MeshRefinementFunction option of FEMToElementMesh

I'm trying to solve Stokes equation in the following domain

The border of the inner half circle where I define a certain boundary condition different from zero and where I am interested in obtaining a smooth solution for the pressure. To obtain a smooth result I need to use a fairly small grid on the centre of the domain an my plan was to use position-dependent grid to speed up the computation.

While checking different functions to use for the MeshRefinementFunction I found the following:

<< NDSolveFEM
Table[
z0 = 0;
r2 =
RegionIntersection[
RegionDifference[
Cuboid[{-L, -L}, {L, L}],
RegionUnion[Disk[{0, -z0}, 1], Disk[{0, z0}, 1]]],
Disk[{0, 0}, L]];
mesh =
ToElementMesh[r2,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.00125 (1 + Norm[Abs[Mean[vertices]] - {0, 0}])^2]];
Show[mesh["Wireframe"], PlotRange -> {{0, 2}, {-1, 1}}],
{L, 10., 20, 5}]


See how the resolution of the inner circle is lost. Something similar happened here

However, if I define the MaxCellMeasure the problem disappears

mesh = ToElementMesh[ r2, "MaxCellMeasure" -> .005]


The problem seems to be the ratio between system size and grid resolution since the same happens in the first case if I decrease the refinement constant for a given system size.

Using MaxBoundaryCellMeasure is not a solution btw, since it introduces other problems of its own. Using such a small grid uniformly makes the computation very time consuming and I don't really need such resolution on the outer half circle.

So, is this is a bug? Is there an easy way to avoid this?

Update

Try using AccuracyGoal

Table[

z0 = 0;
r2 = RegionIntersection[
RegionDifference[Cuboid[{-L, -L}, {L, L}],
RegionUnion[Disk[{0, -z0}, 1], Disk[{0, z0}, 1]]],
Disk[{0, 0}, L]];

mesh = ToElementMesh[r2,
MeshRefinementFunction -> Function[{vertices, area},
area > 0.00125 (1 + Norm[Abs[Mean[vertices]] - {0, 0}])^2],
AccuracyGoal -> 2];


The issue seems to be related to the automatic improvement of boundary position.

Table[

z0 = 0;
r2 = RegionIntersection[
RegionDifference[Cuboid[{-L, -L}, {L, L}],
RegionUnion[Disk[{0, -z0}, 1], Disk[{0, z0}, 1]]],
Disk[{0, 0}, L]];

mesh = ToElementMesh[r2,
MeshRefinementFunction -> Function[{vertices, area},
area > 0.00125 (1 + Norm[Abs[Mean[vertices]] - {0, 0}])^2],
"ImproveBoundaryPosition" -> False];

(*mesh["Wireframe"]*)

Show[mesh["Wireframe"], PlotRange -> {{0, 2}, {-1, 1}}], {L, 10., 20,
5}]


As a solution, I've often used DiscretizeRegion (which also eliminates those annoying extrapolation warnings)

Table[

z0 = 0;
r2 = RegionIntersection[
RegionDifference[Cuboid[{-L, -L}, {L, L}],
RegionUnion[Disk[{0, -z0}, 1], Disk[{0, z0}, 1]]],
Disk[{0, 0}, L]];

mesh = DiscretizeRegion[r2,
MeshRefinementFunction -> Function[{vertices, area},
area > 0.00125 (1 + Norm[Abs[Mean[vertices]] - {0, 0}])^2],
Method -> "Continuation", AccuracyGoal -> 5, PrecisionGoal -> 5];

(*mesh["Wireframe"]*)

Show[mesh, PlotRange -> {{0, 2}, {-1, 1}}], {L, 10., 20, 5}]


• your first method, if I add a small value for "MaxBoundaryCellMeasure" seems to solve the problem quite nicely. As for the annoying warnings, I'm using "ExtrapolationHandler" -> {# &, "WarningMessage" -> False} and seems to work quite well. Thanks for the answer! Aug 12, 2016 at 19:43
• when setting AccuracyGoal the mesh is created nicely but the result is completely unphysical. It also throws a warning about comparing real and imaginary numbers. Any clues why? I can update the code to show that if you are interested. Aug 12, 2016 at 19:56
• @tsuresuregusa sure, I'm interested Aug 12, 2016 at 20:58
• took me a while to find a minimal working example and while doing so discovered that the error with AccuracyGoal appears only when I define the interpolation function as I did in my previous comment. If I delete that, the accuracyGoal solution is better than the one with MaxBoundaryCellMeasure. Both have ugly discontinuities when looking at the inner circle though, so some method of extrapolation is necessary. Tomorrow I will check your discretizeRegion method. Aug 12, 2016 at 23:11
• @tsuresuregusa using DiscretizeRegion is not the same as using ToElementMesh careful here. The difference is that DiscretizeRegion will only produce a first order accurate mesh while ToElementMesh` will by default generate a second order mesh which in turn will result in a better FEM solution. See the ElementMesh generation tutorial for more details. Aug 16, 2016 at 13:17