# Finite element mesh coordinate lies outside solution interpolation function

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["NDSolveFEM"];
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))),
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*)
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[#[], #[]] & /@ 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[];
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?

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.

• Thanks. Although BoundaryDiscretizeRegion is in help I do find that there are a bewildering number of options and suboptions for BoundaryMeshGenerator. I am not sure if I know what to use when.
– Hugh
Sep 21 '17 at 9:55
• @Hugh, I understand that this is not optimal and I apologize for that. The "BoundaryDiscretizeRegion" boundary mesh generator takes all options BoundaryDiscretizeRegion does. The RegionPlot boundary mesh generator takes RegionPlot options with the exception of PlotPoints being called SamplePoints. Thanks again for posting, this will give me a change to look this in more detail and hopefully come up with a fix. Sep 21 '17 at 10:08