4
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ Post your code in copyable format. $\endgroup$ – David G. Stork Feb 8 '16 at 4:16
  • $\begingroup$ Yes, I just tried to put it in "code" format. $\endgroup$ – Majid Hasan Feb 8 '16 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 Feb 8 '16 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 Feb 8 '16 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 Feb 8 '16 at 13:47
3
$\begingroup$

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
$\endgroup$
1
  • $\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$ – Majid Hasan Feb 13 '16 at 4:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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