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["NDSolve`FEM`"]; 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@
   Thread[MeshCells[ΩSubset, 3], Tetrahedron], 
    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.

  • $\begingroup$ Post your code in copyable format. $\endgroup$ Commented Feb 8, 2016 at 4:16
  • $\begingroup$ Yes, I just tried to put it in "code" format. $\endgroup$ Commented Feb 8, 2016 at 4:18
  • 1
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Feb 8, 2016 at 5:15
  • 1
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Feb 8, 2016 at 5:36
  • 2
    $\begingroup$ 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/… $\endgroup$
    – Michael E2
    Commented Feb 8, 2016 at 13:47

1 Answer 1


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]

  • $\begingroup$ 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. $\endgroup$ Commented Feb 13, 2016 at 4:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.