4

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@
 Pick[First@
   Thread[MeshCells[ΩSubset, 3], Tetrahedron], 
  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.

Community
  • 1
Majid Hasan
  • 41
  • 4
  • Post your code in copyable format. – David G. Stork Feb 08 '16 at 04:16
  • Yes, I just tried to put it in "code" format. – Majid Hasan Feb 08 '16 at 04:18
  • 1
    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. – MarcoB Feb 08 '16 at 05:15
  • Yes, i know, it's just that sometimes i need to flatten the grid, so I leave it inside Flatten. – Majid Hasan Feb 08 '16 at 05:21
  • 1
    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. – MarcoB Feb 08 '16 at 05:36
  • In the first code i2min is not defined, so the Iterator {i1,1,101-i2min} fails. In the second code section there are syntax problems. Please [edit] your question. – rhermans Feb 08 '16 at 10:03
  • 2
    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 http://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementProgramming.html – Michael E2 Feb 08 '16 at 13:47

1 Answers1

3

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
bill s
  • 68,936
  • 4
  • 101
  • 191
  • 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. – Majid Hasan Feb 13 '16 at 04:46