Finite element mesh coordinate lies outside solution interpolation function

Question

I solve a finite element problem and then want to look at values at mesh nodes for further calculations. Unfortunately the solution interpolation function identifies some mesh coordinates as being outside the solution interpolation function. This is rather similar to this problem but that was associated with a refinement function which is not used here. It is also similar to this question but in a comment to the question where the problem was first noticed user21 is skeptical that it is the same.

To get going we need a stress operator a mesh and a solution. Here they are

Needs["NDSolve`FEM`"];
planeStress[
   Y_, ν_] := {Inactive[
      Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(2 \
(1 - ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
    Inactive[
      Div][({{-(Y/(1 - ν^2)), 
         0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[
         Grad][u[x, y], {x, y}]), {x, y}], 
   Inactive[
      Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y \
ν)/(1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
     Inactive[
      Div][({{-((Y (1 - ν))/(2 (1 - ν^2))), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        v[x, y], {x, y}]), {x, y}]};

L = 2/10;(*Bracket side length*)
d = 2/100;(*End edge lenght*)
r = L - d;(*radius of curved edge*)
r1 = d;(*radius of hole*)
L2 = 4/1000;(*nail thickness*)
L3 = 4/100;(*nail length*)
L4 = L/8;(*location of nail from bottom*)
Y = 10^3;(*modulus of elasticity*)
ν = 33/100;(*Poisson ratio*)

reg = RegionDifference[
  RegionUnion[
   Rectangle[{0, 0}, {L, L}],
   Rectangle[{-L3, L4 - L2/2}, {0, L4 + L2/2}],
   Rectangle[{-L3, (L - L4) - L2/2}, {0, (L - L4) + L2/2}]
   ],
  RegionUnion[
   Disk[{L, 0}, r],
   Disk[{2 d, L - 2 d}, r1]]
  ];
mesh = ToElementMesh[reg, {{-L3, L}, {0, L}}, 
   "MaxBoundaryCellMeasure" -> 0.001, MaxCellMeasure -> 1.1*^-6, 
   AccuracyGoal -> 8, MeshQualityGoal -> 1, 
   "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 200}];
s1 = -0.5;
x =.; y =.;
{uif, vif} = 
  NDSolveValue[{planeStress[Y, ν] == {0, NeumannValue[s1, y == L]},
    (* top nail *)
    DirichletCondition[
     u[x, y] == 0, -L3 < x < -L3/3 && y == L4 + L2/2], 
    DirichletCondition[
     u[x, y] == 0, -L3 < x < -L3/3 && y == L4 - L2/2],
    DirichletCondition[v[x, y] == 0, 
     x == -L3 && L4 - L2/2 < y < L4 + L2/2],

    (* bottom nail *)
    DirichletCondition[
     u[x, y] == 0, -L3 < x < -L3/3 && y == (L - L4) + L2/2], 
    DirichletCondition[
     u[x, y] == 0, -L3 < x < -L3/3 && y == (L - L4) - L2/2],
    DirichletCondition[v[x, y] == 0, 
     x == -L3 && (L - L4) - L2/2 < y < (L - L4) + L2/2]

    }, {u, v}, {x, y} ∈ mesh];

Now we extract the coordinates of the mesh and put them into the solution. Errors are identified.

coords = mesh["Coordinates"];
uu = uif[#[[1]], #[[2]]] & /@ coords;

(*
InterpolatingFunction::femdmval: Input value {2.10684*10^-17,0.177219} lies outside the range of data in the interpolating function.

InterpolatingFunction::femdmval: Input value {7.45819*10^-18,0.177018} lies outside the range of data in the interpolating function.

*)

I now find these points that failed and plot the mesh and failed coordinates.

nn = Flatten@Position[NumberQ[#] & /@ uu, False];
cc = coords[[nn]];
eps = 0.002;
{x1, y1} = cc[[1]];
Show[mesh["Wireframe"], 
 PlotRange -> {{x1 - eps, x1 + eps}, {y1 - eps, y1 + eps}},
 Epilog -> {Red, PointSize[0.02], Point[cc]}]

The points lie in a corner which the mesh has not followed. However I am not worried by that approximation.

Interestingly if the option "ImproveBoundaryPosition" -> False is added to ToElementMesh then the problem goes away. However, user21 is skeptical about that solution.

What is happening?

Answer 1

Thanks very much for posting this, so I can have a look at this. It seems you were right. The mid side nodes get moved too much in some cases and the algorithm well need to become better for that. Part of the reason is that the mesh becomes fine and mid side nodes get moved more than the elements radius for example.

Concerning workarounds: Since you use a fine grid anyways the benefit of curved elements is not that high so using "ImproveBoundaryPosition" -> False seem OK.

As a second alternative you can use (Method 2) from the other post you mention:

mesh = ToElementMesh[reg, RegionBounds[reg], 
   "BoundaryMeshGenerator" -> {"BoundaryDiscretizeRegion"}, 
   "MaxBoundaryCellMeasure" -> 0.001, MaxCellMeasure -> 1.1*^-6];

This behaves better. I hope this helps.