Build an InterpolatingFunction with ElementMeshInterpolation

Question

Suppose you have a function

f[t_, x_, y_] := t (1 - x^2 - y^2)

defined for {x,y} in Disk[] and an ElementMesh of that Region

<< NDSolve`FEM`
mesh = ToElementMesh[Disk[]]

You can evaluate the function for some t at the vertices and nodes of the mesh:

coords = mesh["Coordinates"];
values = f[0, Sequence@@#]& /@ coords;

and then build an InterpolatingFunction with:

ElementMeshInterpolation[{mesh}, values]

Now suppose you have a list of times:

tl = Range[0, 1, 0.05];

you can evaluate the function for all times and all mesh nodes:

values = Table[f[t, Sequence @@ ##] & /@ coords, {t, tl}];

But it is possible to build a single InterpolatingFunction from these data, accepting x,y,t as arguments?

I know I can do somethiong like

Interpolation@Flatten[Table[{t, Sequence @@ xy, f[t, Sequence @@ xy]}, {t, 
   tl}, {xy, coords}], 1]

ignoring some warning, but I don't think this method can take advantage of a previously built 2nd order ElementMesh, so I'm searching another way, basically the same way used when NDSolve handle a transient equation.

Answer 1

Somewhat manual construction:

<< NDSolve`FEM`
mesh = ToElementMesh[Disk[], MaxCellMeasure -> Infinity];
tl = Range[0, 1, 0.2];
mesh3 = MeshOrderAlteration[
   With[{gc = ElementMeshToGraphicsComplex@MeshOrderAlteration[mesh, 1]},
    gc /. GraphicsComplex[pts_, g_, o___] :>
      ToElementMesh[
       "Coordinates" -> Flatten[
         PadLeft[pts, {Automatic, 3}, #] & /@ tl,
         1],
       "MeshElements" -> {TetrahedronElement[
          Flatten[
           Table[
             Cases[g, 
              Polygon[
                lis_] :> (lis /. 
                 tri : {_Integer, _Integer, _Integer} :>
                  With[{p = Join[tri + (n - 1) Length@pts, tri + n Length@pts]},
                   {p[[{1, 4, 5, 6}]], p[[{1, 5, 2, 6}]], p[[{2, 3, 1, 6}]]}
                   ]),
              Infinity],
             {n, Length[tl] - 1}],
           3]
          ]}
       ]
    ], 2];

mesh3["Wireframe"]

f[t_, x_, y_] := t (1 - x^2 - y^2)
values = f @@ Transpose[mesh3["Coordinates"]];
emifn = ElementMeshInterpolation[{mesh3}, values];

ContourPlot3D[f[t, x, y] == 0.1, {t, 0.`, 1.`}, {x, -1., 1.}, {y, -1., 1.}]

Notes

The ordering for the tetrahedra was deduced from this decomposition of a prism on a triangular base:

foo = ToElementMesh[
   DiscretizeGraphics[
    gc = GraphicsComplex[
      {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}},
      Polygon[{{1, 2, 3}, {6, 5, 4}, {1, 2, 5, 4}, {2, 3, 6, 5}, {3, 1, 4, 6}}]]
    ],
   MaxCellMeasure -> Infinity, "MeshOrder" -> 1];
foo["Coordinates"]
foo["MeshElements"]
(*
  {{0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {0., 1., 1.}, {1., 0., 1.}, {0., 0., 1.}}

  {TetrahedronElement[{{1, 6, 5, 4}, {1, 5, 2, 4}, {2, 3, 1, 4}}]}
*)

I tried RegionProduct but I either got errors or could not make the resulting mesh be the product of the meshes.

Reducing the mesh order to 1, forming the product, then bumping back up to 2 was easier than figuring out how to order the vertices of an order-2 tetrahedron. Surely it can be done, but this seemed easier/lazier.

Answer 2

I like the answer of @Michael E2 but I'm still searching, if available, a simpler way. Experimenting I found this way. My way probably can be applied also for a 3D mesh/region.

The data to interpolate:

<< NDSolve`FEM`
f[t_, x_, y_] := t (1 - x^2 - y^2);
mesh = ToElementMesh[Disk[]];
coords = mesh["Coordinates"];
tl = Range[0, 1, .05];
values = Table[f[t, Sequence @@ ##] & /@ coords, {t, tl}];

Building the interpolating function:

if = NDSolveValue[{D[u[t, x, y], t] == 0, u[0, x, y] == 0}, u, 
     Element[{x, y}, mesh], {t, tl[[1]], tl[[-1]]}, 
     Method -> "FiniteElement"]

if = If[if[[2, 1]] != 5, $Failed, 
  ReplacePart[if, {
    {2, 3} -> 0, (* no derivatives provided *)
    {2, 4, 1} -> Length@tl,(* actual grid points along t *)
    {3, 1} -> tl, (* actual grid along t *)
    {4} -> List /@ values (* values to interpolate *)
    }]]

Apprently this works:

Manipulate[
 Plot3D[if[t, x, y], {x, y} \[Element] mesh, 
  PlotRange -> {{-1, 1}, {-1, 1}, {0, 1}}], {t, 0, 1}]

I'm not completely satisfied with this approach because it relies on the undocumented internal structure of the InterpolatingFunction (partially revealed by this question, thanks!) and I completely ignore if it's "safe" to "tamper" with this internal structure (maybe Interpolation or NDSolve build also some private data accompaining the returned InterpolatingFunction?).

NOTE

I'm more confident with the the time variable after the space variables x,y so I also tried with:

NDSolveValue[{D[u[x, y, t], t] == 0, u[x, y, 0] == 0}, u, 
  Element[{x, y}, mesh], {t, tl[[1]], tl[[-1]]}, 
  Method -> "FiniteElement"]

Unfortunately I get this message

Drop::drop: Cannot drop positions 3 through 3 in {ElementMesh[{{1.,5.26062*10^-22},{0.990147,0.14003},{0.960783,0.2773},{0.91339,0.407085},<<43>>,{0.990147,-0.14003},{-0.601496,-0.540192},{-0.725006,-0.512798},<<951>>},<<2>>,{PointElement[{{1},<<49>>,<<46>>}]}],{<<1>>}}. >>

and soon the Kernel crahes. I think in the past I successfully used this standard...

UPDATE

Thanks to @user21 comment the following approach will be available with the next release.

values = Partition[values, 1];
Dimensions[values] == {Length@tl, 1, Length@mesh["Coordinates"]}
if = ElementMeshInterpolation[{tl, mesh}, values]