# Does Interpolation support multi-dimension input and output?

I have an interpolation problem with both multi-dim input and output. Here is an example of the error:

Interpolation[{{{316.,
228.}, {0.9587445571260307, -0.019395144189562147}}, {{161.,
424.}, {1.070025266546878, 0.040399144303199186}}, {{191.,
298.}, {1.0690628597419136, 0.01710950051961435}}}]


Instead of giving me an interpolation function, the output simply repeats the input:

Interpolation[{{{316., 228.}, {0.958745, -0.0193951}}, {{161.,
424.}, {1.07003, 0.0403991}}, {{191., 298.}, {1.06906,
0.0171095}}}]


When I use a single dim output with a similar input format, the Interpolation function works fine:

InterpolatingFunction[{{82., 353.}, {342., 439.}}, {
5, 4225, 0, {3, 3}, {2, 2}, 0, 0, 0, 0, Automatic, {}, {}, False}, {
NDSolveFEMElementMesh[{{84., 421.}, {82., 342.}, {353., 439.}}, {
NDSolveFEMTriangleElement[{{2, 3, 1}}]}, {
NDSolveFEMLineElement[{{3, 1}, {1, 2}, {2,
3}}]}]}, {0.9760543173343315, 0.964633671206072, \
0.9573634991631621}, {Automatic}]


Does that mean the Interpolation function only works for either multi-dim input or output, not both?

No, interpolating vectors is possible, but I think that when the data points form an unstructured input grid, then interpolation is restricted to scalar. At least I didn't succeed.

However on a rectangular grid, it works:

func = Interpolation@
Flatten[Table[{{x, y}, {x^2, y^3}}, {x, 5}, {y, 5}], 1];

func[3, 2]

(*  {9, 8}  *)


Interpolation on an unstructured grid is done through ElementMeshInterpolation. Here is some evidence that interpolation of vectors is unsupported:

Needs["NDSolveFEM"];
mesh = ToElementMesh[Cuboid[{0, 0}], MaxCellMeasure -> 1];
values = Function[{x, y}, x^2 + y^2] @@@
mesh["Coordinates"];

ElementMeshInterpolation[{mesh}, values]         (* succeeds *)
ElementMeshInterpolation[{mesh}, List /@ values] (* succeeds *)
ElementMeshInterpolation[{mesh}, Transpose@{values, values}] (* fails *)
(*
InterpolatingFunction[{{0., 1.}, {0., 1.}}, "<>"]

InterpolatingFunction[{{0., 1.}, {0., 1.}}, "<>"]

ElementMeshInterpolation[{ElementMesh[{{0., 1.}, {0.,
1.}}, {TriangleElement["<" 2 ">"]}]}, {{0., 0.}, {1., 1.}, {2.,
2.}, {1., 1.}, {0.25, 0.25}, {0.25, 0.25}, {0.5, 0.5}, {1.25,
1.25}, {1.25, 1.25}}]
*)


There is no error message with the last one, which leaves open the possibility that it should have worked. Or they forgot to code a message for this case. I'm not sure if the second example with 1-vectors was intended to work or not; the result is same as the first example: scalar output.