10
$\begingroup$

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.

$\endgroup$
8
  • 1
    $\begingroup$ I guess you could build a list of 2D interpolating functions using ElementMeshInterpolation, and then interpolate between their values in t? $\endgroup$ – user484 May 6 '15 at 0:54
  • $\begingroup$ I tried to look at the documentation for the ElementMeshInterpolation and failed. Could you kindly give a hint, where and how can I get it? $\endgroup$ – Alexei Boulbitch May 6 '15 at 8:12
  • $\begingroup$ @AlexeiBoulbitch You can find some sample here but unfortunately all the samples are only for stationary problems (search for "post-process"). A time-dependent case is missing. $\endgroup$ – unlikely May 6 '15 at 8:16
  • $\begingroup$ Is there a reason why you don't want to build a 3D mesh of the whole region directly at the start and then interpolate over that? Otherwise, doing that seems like the simplest option. $\endgroup$ – Virgil May 6 '15 at 14:49
  • $\begingroup$ @Virgil What if you want to extend to a f[t, x, y, z]? $\endgroup$ – unlikely May 6 '15 at 15:55
6
$\begingroup$

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.}]

Mathematica graphics

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.

$\endgroup$
5
$\begingroup$

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}]

Mathematica graphics

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]
$\endgroup$
13
  • $\begingroup$ Clearly something like this is done for solutions to transient PDEs. I kept thinking I could hit on the right combinations of arguments to ElementMeshInterpolation to achieve it; but I'm thinking now it is done internally on such a low level as you've done, and there is no FEM` function that we can use. $\endgroup$ – Michael E2 May 6 '15 at 15:35
  • $\begingroup$ @MichaelE2, I added an example to the FiniteElementProgramming tutorial that shows how to do it. That will be available in the next release. If you do have FEM related questions, feel free to @ me. $\endgroup$ – user21 May 8 '15 at 9:16
  • $\begingroup$ @user21 Good news! $\endgroup$ – unlikely May 8 '15 at 9:46
  • $\begingroup$ @user21 Thank you very much. $\endgroup$ – Michael E2 May 8 '15 at 9:57
  • 2
    $\begingroup$ @unlikely, ElementMeshInterpolation[{tRange, mesh}, data] will work, Where Length[mesh["Coordinates"]] is A and Length[tRange] B and Dimensions[data] is {B,1,A} $\endgroup$ – user21 May 18 '15 at 15:33

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.