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?