# Finite Element Method and ElementMeshInterpolation issue in v10.3

Bug introduced in 10.3

I recently installed Mathematica v10.3 and checked a project I was working on few months ago under Mathematica v10.1.

In this project I solve a PDE with NDSolve and Finite Element Method over a 2D anulus region and I and ContourPlot that solution.

Under v10.1 (and v10.2) the solution appear like this:

Under v10.3 the solution appear like this:

The source code is exactly the same.

At the beginning I was convinced the problem was with NDSolve (and this question reflected that, I later edited). After more investigation I discovered the solution in the two releases is nearly the same.

You can reproduce the problem running:

Import["http://bit.ly/cp103issue"]


I searched on the documentation, news section, but I didn't find anything that could explain this or restore the previous behavior.

Any chance to understand the reason of the changed behavior and eventually restore the old?

# UPDATE

The problem is not related to ContourPlot but to the InterpolatingFunction returned by NDSolve (thanks to JasonB). Probably the InterpolatingFunction is built under the cover by a call to ElementMeshInterpolation. Indeed the following code in Mathematica 10.3 returns result completely different from the results in Mathematica 10.2.

<<NDSolveFEM
{Ri,Ro}={3000000.,6371000.};
rgn=RegionDifference[Disk[{0.,0.},Ro],Disk[{0.,0.},Ri]];
mesh=ToElementMesh[rgn];
grid=mesh["Coordinates"];
vog=(Exp[-Norm[#]^2/Ri^2])&/@grid;
if=ElementMeshInterpolation[{mesh},vog];
ContourPlot[if[x,y],{x,y}\[Element]mesh,PlotLegends->Automatic,Epilog->{Point[grid]}]
TakeLargest[(if@@@grid-vog),20]


In 10.3:

In 10.2 (after som InterpolatingFunction::dmval message):

As you can see, in 10.2, the largest differences in "interpolation errors" are negligible, while in 10.3 they are even greather than the actual expected value.

The problem appear more specifically related with the "scale" of the mesh. Indeed if you change the sizes:

{Ri,Ro}={3000000.,6371000.}/10^6;


all apear fine:

• Unless you can post some kind of example it is hard to figure it out, sorry – Jason B. Feb 18 '16 at 10:51
• @JasonB Of course, just asking if someone is aware of some change in the FEM framework between 10.2 and 10.3... I'll try to prepare a reduced example if possible... – unlikely Feb 18 '16 at 11:11
• What loading/deformation? Plane stress 2 point, 4 point forces deformation/ stress plots? Solves an isotropic 2D stress function pde? – Narasimham Feb 18 '16 at 14:16
• @Narasimham Yes, the plots shoud represents the modulus and horizontal component of the displacement vector; there is a load on the volume, horizontal outward on the x axis and vertical inward on the y axis, that gradually canges in-between; there is an essential homogeneus boundary condition on the inner border, and natural homogeneus on the outer. The material is composed of many isothropic anulus with differents Lamé coefficients. – unlikely Feb 18 '16 at 16:09
• The stress/deformation results are essentially the same. Extra features introduced elsewhere perhaps reduce the display/appearance .. that you can ignore until WRI would correct later on. – Narasimham Feb 18 '16 at 16:21

Too long to make as a comment, so I'll post this non-answer as an answer.

I just want to point out that this definitely is not a ContourPlot problem, and it is absolutely a problem with the interpolating function.

So let's evaluate the following code in versions 10.2 and 10.31

Import["http://bit.ly/cp103issue"];
{sol1, sol2} = testSolution;
grid = testSolution[[1]]["Grid"];
vals1 = testSolution[[1]]["ValuesOnGrid"];
vals2 = testSolution[[2]]["ValuesOnGrid"];
norm = Norm /@ Transpose[{vals1, vals2}];
normgrid = Transpose[Append[Transpose[grid], norm]];


So I have extracted the grid points from the interpolating functions, as well as the values on those grid points. I have checked, and they are the same across both versions to within an accuracy of 10^-17. So we should be dealing with the same interpolating functions. But we aren't.

Here is the nail in the coffin,

ListPlot[sol1 @@@ grid - vals1]
ListPlot[sol2 @@@ grid - vals2]


Evaluating this code in version 10.2, we get

which looks great. In 10.3.1 we get

which is just plain horrible. So what we have here is some backsliding in the interpolating function resulting from an NDSolve call using a nonrectangular region.

Now here is something interesting, if you take those grid points and the values on the grid, and you make an interpolation function using Interpolation, then you get basically the same plot in both versions.

ContourPlot[
Interpolation[Transpose[{grid, norm}], InterpolationOrder -> 1][x,
y], {x, y} ∈ testMesh]


• Paging @user21 - here is some more interpolation madness! – Jason B. Feb 19 '16 at 12:40
• Thanks for your investigation. I was near to confirm the same results. So practically speaking, with ListPlot[sol1 @@@ grid - vals1] test you are proving that the InterpolatingFunction doesn't even interpolate on the grid points? – unlikely Feb 19 '16 at 13:43
• @unlikely - yes!! The most basic thing you can expect from an interpolating function. – Jason B. Feb 19 '16 at 13:47
• Thanks Jason; I think paging only works if I had a contribution on this particular site already; in any case I did not get a notification. But I certainly will have a look what goes wrong. – user21 Feb 21 '16 at 3:01

Update: This is fixed in Version 11.0:

<< NDSolveFEM
{Ri, Ro} = {3000000., 6371000.};
rgn = RegionDifference[Disk[{0., 0.}, Ro], Disk[{0., 0.}, Ri]];
mesh = ToElementMesh[rgn];
grid = mesh["Coordinates"];
vog = (Exp[-Norm[#]^2/Ri^2]) & /@ grid;
grid]},PlotLegends\[Rule]Automatic]*)

if = ElementMeshInterpolation[{mesh}, vog];
(*ContourPlot[if[x,y],{x,y}\[Element]mesh,PlotLegends\[Rule]Automatic,\
Epilog\[Rule]{Point[grid]}]*)
TakeLargest[(if @@@ grid - vog), 20]

{4.996*10^-16, 4.85723*10^-16, 4.44089*10^-16, 3.88578*10^-16,
3.33067*10^-16, 3.33067*10^-16, 2.498*10^-16, 2.498*10^-16,
2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16,
2.22045*10^-16, 1.94289*10^-16, 1.66533*10^-16, 1.66533*10^-16,
1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16}


Old message:

This is a bug; you can workaround that by setting

SetSystemOptions[ "FiniteElementOptions" -> {"CacheInterpolationElements" -> False}]


Update:

It seems the default tolerance is not high enough for mesh over a large space compared to the default bounding box:

SetSystemOptions[
"FiniteElementOptions" -> {"InterpolationToleranceFactor" -> 10^-8}]


This would allow you to still make use of cached elements which is faster for structured queries (which *Plot is not)

• Sorry? Wait, is this your fault? Grab the pitchforks. :-) – Jason B. Feb 21 '16 at 18:19
• @JasonB, sometimes things just get out of control.... – user21 Feb 22 '16 at 20:28
• Both workaround works to fix the aspect of my plots. But no one, apparently, fix the numerical problem exposed in the @JasonB answer (InterpolatingFunction that doesn't interpolate), at least in my specific case: ifn=Import["http://bit.ly/1ThEi3G"];ifn["ValuesOnGrid"]-ifn@@@ifn["ElementMesh"]["Coordinates"]//Abs//Max is not 0... It's 10^-2 w.r.t. the order of magnitude of ifn values on grid... – unlikely May 17 '16 at 9:57
• @unlikely, what do you get 0?, what OS, and \$Version? – user21 May 17 '16 at 10:11
• @user21 - I'm using version 10.4 Linux, and here are the results. You have to use the "InterpolationToleranceFactor" option from this answer to get the correct result. – Jason B. May 17 '16 at 11:27