# PeriodicInterpolation does not work in ElementMeshInterpolation

I want to use ElementMeshInterpolation to generate interpolation function with periodic boundary condition.

I use below data as an example

data=Flatten[Table[{i,j,Sin[i+j]},{i,0,2\[Pi],2\[Pi]/50},{j,0,2\[Pi],2\[Pi]/50}],1];
ListContourPlot[data]


which gives This data is periodic along x and y direction.

Using Interpolation

f = Interpolation[data, PeriodicInterpolation -> True];
{ContourPlot[f[x, y], {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}],
ContourPlot[f[x, y], {x, 0, 4 \[Pi]}, {y, 0, 4 \[Pi]}]}


gives We can see the Interpolation function is fine with periodic condition as wanted.

using ElementMeshInterpolation

Though Interpolation works fine for this data set. But Interpolation has problem that it frequently run into "femimq" problem. So ElementMeshInterpolation on a refined mesh is necessary sometimes.

mesh = ToElementMesh[data[[;; , 1 ;; 2]]];
f = ElementMeshInterpolation[{mesh}, data[[;; , -1]],
PeriodicInterpolation -> {True, True}];
{ContourPlot[f[x, y], {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}],
ContourPlot[f[x, y], {x, 0, 4 \[Pi]}, {y, 0, 4 \[Pi]}]}


this gives You see the generated Interpolation function has no periodicity.

using ListInterpolation

mesh can also be used in ListInterpolation

mesh = ToElementMesh[data[[;; , 1 ;; 2]]];
f = ListInterpolation[data[[;; , -1]], mesh,
PeriodicInterpolation -> {True, True}];
{ContourPlot[f[x, y], {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}],
ContourPlot[f[x, y], {x, 0, 4 \[Pi]}, {y, 0, 4 \[Pi]}]}


but this gives the same result as ElementMeshInterpolation.

So the question is how to correctly make periodic interpolation function using ElementMeshInterpolation.

You can not really. The fact that Interpolation can do this hinges on the data being structured. In other works what I am going to show next is not easily generally possible for meshes that represent a non rectangullar domain; which is the common case for FEM meshes.

You can hack it by using the ExtrapolationHandler option.

    Needs["NDSolveFEM"]
data = Flatten[
Table[{i, j, Sin[i + j]}, {i, 0, 2 \[Pi], 2 \[Pi]/50}, {j, 0,
2 \[Pi], 2 \[Pi]/50}], 1];
mesh = ToElementMesh[data[[;; , 1 ;; 2]]];
f = ElementMeshInterpolation[{mesh}, data[[;; , -1]]];


Now, we can use f as a function in the extrapolation handler and map the coordinates outside the domain back onto f. This mapping back to the original domain is tricky to do generally. Here we use Mod.

    f2 = ElementMeshInterpolation[{mesh}, data[[;; , -1]],
"ExtrapolationHandler" -> {Function[{x, y},
f[Mod[x, 2 \[Pi]], Mod[y, 2 \[Pi]]]]}]

{ContourPlot[f2[x, y], {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}],
ContourPlot[f2[x, y], {x, 0, 4 \[Pi]}, {y, 0, 4 \[Pi]}]} If you want to switch of the warning message you can do so with:

f2 = ElementMeshInterpolation[{mesh}, data[[;; , -1]],
"ExtrapolationHandler" -> {Function[{x, y},
f[Mod[x, 2 \[Pi]], Mod[y, 2 \[Pi]]]], "WarningMessage" -> False}]


Perhaps an idea for a future implementation.

• Thank you so much. I am interested in the overhead of "ExtrapolationHandler". I found a peculiar thing. The ContourPlot of f2 is even faster than f. However, Table[Quiet@f2[x,y],{x,100.,150.,1},{y,100.,150.,1}];//AbsoluteTiming is much slower than Table[f[x,y],{x,100.,150.,1},{y,100.,150.,1}];//AbsoluteTiming. Why is that? Feb 4, 2021 at 8:22
• And what do you mean by "Interpolation can do this hinges on the data being structured." Feb 4, 2021 at 8:23
• @matheorem, for your second question see update Feb 4, 2021 at 8:29
• @matheorem, here is a guess, If you are querying f outside it will have to generate an extrapolation value and it has to generate the Message. That might be the cost. If f2 is queried outside then you have the expense of the finding that the value is outside the domain plus the evaluation of f to get the value. Feb 4, 2021 at 8:32
• @matheorem, I have added this example to the documentation. I hope you do not mind. Feb 4, 2021 at 8:33

Since currently ElementMeshInterpolation does not support PeriodicInterpolation and Interpolation only support PeriodicInterpolation on rectangular grid. Apart from user21's workaround, I developed a workaround for arbitrary parallel or parallelipiped grid periodic interpolation.

The idea is naive, just to pull back points outside region by base vectors. Below is helper function.

pullBack2Dcom=Compile[{x1,x2,y1,y2,x,y},Mod[{(-x2 y+x y2)/(-x2 y1+x1 y2),(x1 y-x y1)/(-x2 y1+x1 y2)},1].{{x1,y1},{x2,y2}}];
pullBack3Dcom=Compile[{x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z},Mod[{(x3 y2 z-x2 y3 z-x3 y z2+x y3 z2+x2 y z3-x y2 z3)/(x3 y2 z1-x2 y3 z1-x3 y1 z2+x1 y3 z2+x2 y1 z3-x1 y2 z3),(x3 y1 z-x1 y3 z-x3 y z1+x y3 z1+x1 y z3-x y1 z3)/(-x3 y2 z1+x2 y3 z1+x3 y1 z2-x1 y3 z2-x2 y1 z3+x1 y2 z3),(x2 y1 z-x1 y2 z-x2 y z1+x y2 z1+x1 y z2-x y1 z2)/(x3 y2 z1-x2 y3 z1-x3 y1 z2+x1 y3 z2+x2 y1 z3-x1 y2 z3)},1].{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3}}];
pullBack2D[{{x1_,y1_},{x2_,y2_}},{x_,y_}]:=pullBack2Dcom[x1,x2,y1,y2,x,y];
pullBack3D[{{x1_,y1_,z1_},{x2_,y2_,z2_},{x3_,y3_,z3_}},{x_,y_,z_}]:=pullBack3Dcom[x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z];


The expressions above seems complicated, but they are just solution of LinearSolve plus using Mod function.

Now I prepare a parallel grid data set which has periodic boundary condition.

data=N@Flatten[Table[Append[i{1,0}+j*{1,1},Sin[Total[i{1,0}+j*{1,1}]*2*\[Pi]]],{i,0,1.,1/100},{j,0,1.,1/100}],1];
ListContourPlot[data,AspectRatio->Automatic]


which gives and then do mesh interpolation like below

mesh=ToElementMesh[data[[;;,1;;2]]];
f=ElementMeshInterpolation[{mesh},data[[;;,-1]]];
ContourPlot[f[x,y],{x,y}\[Element]ConvexHullMesh[data[[;;,1;;2]]],AspectRatio->Automatic]


gives which is good.

Now we can use pullBack2D to generate a periodic function using base vector bvecs

bvecs = {{1, 0}, {1, 1}}
g[x_?NumericQ, y_?NumericQ] := f @@ pullBack2D[bvecs, {x, y}]


plot it using

ContourPlot[g[x, y], {x, 0, 2}, {y, 0, 2}, AspectRatio -> Automatic]


we get • You should be able to use a technique like the one here if you have the basis vectors of your "fundamental parallelogram". Feb 4, 2021 at 13:15
• @J.M. Thank you for the link. I see you use LinearSolve and that is exactly what I use except that I use the explict symbolic solution of LinearSolve to use Compile. I found Compile is much faster :) Feb 4, 2021 at 13:19