# Interpolation of data on a 4D grid

I have a 4D grid

l=50; Ω = Table[{(i1 - 9/10)/(2*l), (i2 - 9/10)/(2*l), (i3 - 9/10)/(2*l), (i4 - 9/10)/(2*l),
1 - (i1 - 9/10)/(2*l) - (i2 - 9/10)/(2*l) - (i3 - 9/10)/(2*l) - (i4 - 9/10)/(2*l)}, {i1, 1,  2l}, {i2, 1, 2l - i1 + 1}, {i3, 1, 2l - (i1 + i2)+2}, {i4, 1, 2l - (i1 + i2 + i3)+3}];


and I need to find an interpolating function over this grid. If I use the standard Interpolation function, I get "unstructured grid" error. To work around this problem, I've been using the following code (which I found at Interpolation on a regular square grid spanning a triangular domain) for a 3D grid

Needs["NDSolveFEM"]; Needs["TetGenLink"]; tmpF2T={}; l=50; Ω = Flatten[ Table[{(i1 - 9/10)/(2*l), (i2 - 9/10)/(2*l), (i3 - 9/10)/(2*l),  1 - (i1 - 9/10)/(2*l) - (i2 - 9/10)/(2*l) - (i3 -    9/10)/(2*l)}, {i1, 1, 2l}, {i2, 1, 2l - i1 + 1}, {i3, 1, 2l - (i1 + i2)+2}], 0]; ΩSubset =  DelaunayMesh[Flatten[Ω[[All, All, All, {1, 2, 3}]], 2]]; emesh =  ToElementMesh[ "Coordinates" -> MeshCoordinates[ΩSubset],   "MeshElements" -> {TetrahedronElement@
Pick[First@
Unitize[Chop@
Flatten@PropertyValue[{ΩSubset, 3},
MeshCellMeasure]], 1]}]; AppendTo[tmpF2T,ElementMeshInterpolation[{emesh}, loc[[All, 1]]]];


where loc[[All, 1]] is a list of function values defined over the Ω grid. These function values are largely irrelevant for this question, and can be replaced with, say, Table[50, {i,1,Length[ΩSubset]}]. The code has been working great for 3D. It removes the problematic points that cause issues in interpolation, and then interpolation works just fine.

Now I need to do the same thing for my four dimensional grid. But, DelaunayMesh does not seem to be producing a grid in this case. So my question is, what is the right way to extend this code to four and possibly five dimensions? And should I be using "Tetrahedron" and "TetrahedronElement" for more than 3 dimensions?

Thank you.

• Post your code in copyable format. Commented Feb 8, 2016 at 4:16
• Yes, I just tried to put it in "code" format. Commented Feb 8, 2016 at 4:18
• Majid, what are you trying to achieve by Flatten[..., 0] in your definitions of Omega? That level specification is equivalent to no flattening at all. Commented Feb 8, 2016 at 5:15
• I see. Perhaps you should remove it here, though. It is confusing and it reduces the legibility of your code. More in general, I wonder if you could reduce your problem to a smaller, more manageable size (i.e. a minimal working example). That will be more likely to attract attention and generate good answers. Commented Feb 8, 2016 at 5:36
• I don't think the built-in meshers can go beyond 3D, though the documentation is vague on limitations. Maybe someone will correct me. Example: "DiscretizeRegion::cdim: The region given at position 1 in DiscretizeRegion[Cuboid[{0,0,0,0},{1,1,1,1}]] is in dimension 4. DiscretizeRegion only supports dimensions 1 through 3." One can get a 4D {x, y, z} x {t} interpolation; see "Transient PDEs" in reference.wolfram.com/language/FEMDocumentation/tutorial/… Commented Feb 8, 2016 at 13:47

The problem is your data is not on a simple grid. Let's simplify a bit and we'll see that Interpolation can handle 4-dimensional data.

l = 5;
f[i1_, i2_, i3_, i4_] := i1 + i2 + i3 + i4;
omega = Table[{{i1, i2, i3, i4}, f[i1, i2, i3, i4]},
{i1, l}, {i2, l}, {i3, l}, {i4, l}];


Now if you look at omega, you'll see that the parentheses are all wrong for the form needed to input into Interpolation. This can be fixed by judicious flattening:

g = Interpolation[Flatten[omega, 3]]


Hence, for example,

g[2, 2.1, 2.2, 1.5]

7.8

• But the main issue I am facing with the standard Interpolation function is that it returns me "Interpolation::indim: The coordinates do not lie on a structured tensor product grid." error, when I define omega as: omega = Table[{{i1, i2, i3, i4}, f[i1, i2, i3, i4]}, {i1, l}, {i2, l-(i1)+1}, {i3, l-(i1+i2)+2}, {i4, l-(i1+i2+i3)+3}]; g = Interpolation[Flatten[omega, 3]];` This is why I am using that bit of code shown in my question. Commented Feb 13, 2016 at 4:46