13
$\begingroup$

as I'm not satisfied with the way Mathematica interpolates data on unstructured grids, I would like to implement my own interpolation. I will aim for Crouzeix - Raviart type of (quadratic) triangular elements. Let $v_1$, $v_2$ and $v_3$ be nodes of a triangle element $f$. Let $\lambda_i$ be linear interpolant functions for element $f$, thus satisfying

$$ \lambda^{(f)}_i (v^{(f)}_j) = \delta_{ij} $$ for $v_j$ that belongs to the element $e$ and zero for any other points in mesh.

The whole 2D (3D) mesh is constructed from triangular (tetrahedral) elements and each element gives rise to three (four) linear interpolants. So far so good, when we take a set of points in 2D or 3D and perform Delaunay triangulation upon it, we can approximate any function in respective domain by its $\lambda$-interpolation (linear). Nodal values are exact values of the said function in points $v_i$, values in any other points are obtained by the following expression:

$$ f(x, y) = \sum_f \sum_{i = 1}^3 f (v^{(f)}_i) \, \lambda^{(f)}_i (x, y) $$ (sum over $f$ is sum over triangular elements)

Basically, this is what I want to implement.

Now let's take it one step further - the Crouzeix - Raviart elements I'm interested in.

In addition to nodal points of each triangle, there are also edge midpoints and triangle centers - together seven points. Let the edge midpoint opposite to nodal point $i$ be labeled as $v_{i+3}$ (so the point between $v_2$ and $v_3$ is $v_4$ etc). Triangle center is labeled as $v_7 = (1/3) (v_1 + v_2 + v_3)$.

Now the appropriate interpolants are as follows

$$ \begin{aligned} \varphi_1 &= \lambda_1 (2 \lambda_1 - 1) + 3 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_2 &= \lambda_2 (2 \lambda_2 - 1) + 3 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_3 &= \lambda_3 (2 \lambda_3 - 1) + 3 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_4 &= 4 \lambda_2 \lambda_3 - 12 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_5 &= 4 \lambda_3 \lambda_1 - 12 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_6 &= 4 \lambda_1 \lambda_2 - 12 \lambda_1 \lambda_2 \lambda_3 \\ \varphi_7 &= 27 \lambda_1 \lambda_2 \lambda_3 \end{aligned} $$

And, again, the expression for interpolation of any function defined in nodal points + edge midpoints + triangle centers is as follows: $$ f(x, y) = \sum_f \sum_{i = 1}^7 f (v^{(f)}_i) \, \varphi^{(f)}_i (x, y) $$

I have no clue how to start on this one.

In Mathematica, one can define linear interpolands as follows:

ClearAll[a, b, c, x1, x2, x3, y1, y2, y3]
lambda1 = 
  First@Solve[{a x1 + b y1 + c == 1, a x2 + b y2 + c == 0, 
     a x3 + b y3 + c == 0}, {a, b, c}];
lambda2 = 
  First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 1, 
     a x3 + b y3 + c == 0}, {a, b, c}];
lambda3 = 
  First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 0, 
     a x3 + b y3 + c == 1}, {a, b, c}];
λ1 = (a x + b y + c /. lambda1);
λ2 = (a x + b y + c /. lambda2);
λ3 = (a x + b y + c /. lambda3);

And quadratic interpolands just by the definition:

φ1 = λ1 (2 λ1 - 1) + 
    3 λ1 λ2 λ3 // Simplify;
φ2 = λ2 (2 λ2 - 1) + 
    3 λ1 λ2 λ3 // Simplify;
φ3 = λ3 (2 λ3 - 1) + 
    3 λ1 λ2 λ3 // Simplify;
φ4 = 
  4 λ2 λ3 - 12 λ1 λ2 λ3 // 
   Simplify;
φ5 = 
  4 λ3 λ1 - 12 λ1 λ2 λ3 // 
   Simplify;
φ6 = 
  4 λ1 λ2 - 12 λ1 λ2 λ3 // 
   Simplify;
φ7 = 27 λ1 λ2 λ3 // Simplify;

You can work with the following simple mesh:

ImportString["g Domain2D

v -0.45 -1
v 0.45 -1
v 1 -0.45
v 1 0.45
v 0.45 1
v -0.45 1
v -1 0.45
v -1 -0.45
v 0 0

f 1 2 9
f 2 3 9
f 3 4 9
f 4 5 9
f 5 6 9
f 6 7 9
f 7 8 9
f 8 1 9", "OBJ"]

(just save it to some file with obj extension and import it to Mathematica to see it)

and the following data to be interpolated (defined in triangle nodes, edge midpoints and centers):

data = {{-0.45, -1., -0.134552}, {0.45, -1., -0.783103}, {1., -0.45, \
-0.191574}, {1., 0.45, 3.39192}, {0.45, 1., 0.185477}, {-0.45, 
  1., -1.10313}, {-1., 0.45, -1.56601}, {-1., -0.45, 2.01748}, {0., 
  0., 1.}, {0., -1., -0.416147}, {-0.225, -0.5, 
  0.577764}, {0., -0.666667, 
  0.235238}, {0.725, -0.725, -0.568926}, {0.225, -0.5, 
  0.530302}, {0.483333, -0.483333, 0.520129}, {1., 0., 
  2.31016}, {0.5, -0.225, 1.20265}, {0.666667, 0., 1.84406}, {0.725, 
  0.725, 1.65813}, {0.5, 0.225, 1.73818}, {0.483333, 0.483333, 
  1.48638}, {0., 1., -0.416147}, {0.225, 0.5, 0.925275}, {0., 
  0.666667, 0.235238}, {-0.725, 0.725, -1.353}, {-0.225, 0.5, 
  0.182791}, {-0.483333, 0.483333, -0.215702}, {-1., 0., 
  0.627223}, {-0.5, 0.225, 0.291137}, {-0.666667, 0., 
  0.607318}, {-0.725, -0.725, 0.874062}, {-0.5, -0.225, 
  0.826662}, {-0.483333, -0.483333, 0.750545}}
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Encoding your example OBJ file like this is probably more convenient for others. $\endgroup$ – Michael E2 Feb 10 '18 at 16:27
  • $\begingroup$ That's great, I was thinking how to make the code so it works immediately. $\endgroup$ – user16320 Feb 10 '18 at 20:18
  • $\begingroup$ I was wondering if you could elaborate a bit on why you are unsatisfied with Mathematica's way of interpolating on unstructured grids. Have you done a comparison with the method given in the answer? What were your findings? $\endgroup$ – user21 Dec 20 '18 at 6:43
10
$\begingroup$

First, I saved the OP's mesh data in chat:

foo = Import[  (* fetches the chat message *)
   "https://chat.stackexchange.com/transcript/message/42789966#42789966",
   {"HTML", "XMLObject"}];
mreg = First@First@ToExpression@   (* runs ImportString[] - see chat *)
    Cases[foo, XMLElement["pre", _, t_] :> t, Infinity]

The undocumented function Region`Mesh`MeshMemberCellIndex can be used to fetch the cell a point lands in:

cellI = Region`Mesh`MeshMemberCellIndex[mreg, {-0.1, 0.01}]
(*  {2, 7}  *)

The following constructs the data structures we will need for the interpolation. These should be put in a function that returns a crInterpolatingFunction (see further down):

nodes = Join[MeshCoordinates[mreg], 
  MeshCells[mreg, 1] /. Line[p_] :> Mean@MeshCoordinates[mreg][[p]], 
  MeshCells[mreg, 2] /. 
   Polygon[p_] :> Mean@MeshCoordinates[mreg][[p]]
  ];
nodesNF = Nearest[nodes -> "Index"];
valsOnNodes = Flatten[Nearest[data[[All, {1, 2}]] -> data[[All, 3]], nodes], 1];
crElements = MeshCells[mreg, 2] /. Polygon[p_] :> Flatten[{
     p,
     RotateLeft@nodesNF@MovingAverage[
        MeshCoordinates[mreg][[Append[#, First[#]] &@p]],
        2],
     nodesNF@Mean@MeshCoordinates[mreg][[p]]
     }]
(*
  {{1, 2, 3, 11, 12, 10, 26}, {2, 4, 3, 14, 11, 13, 27}, {4, 5, 3, 16, 14, 15, 28},
   {5, 6, 3, 18, 16, 17, 29}, {6, 7, 3, 20, 18, 19, 30}, {7, 8, 3, 22, 20, 21, 31},
   {8, 9, 3, 24, 22, 23, 32}, {9, 1, 3, 12, 24, 25, 33}}
*)

Notes:

  • nodes joins three arrays of points, the vertices, the edge midpoints, and the centroids.

  • nodesNF is for mapping coordinates of nodes in a CR element to their indices in the array nodes.

  • valsOnNodes maps each node $v$ in nodes to the value data[[All, 3]] corresponding to the element in data[[All, {1,2}]] which matches (or is "nearest") the node $v$. It is important that the order of valsOnNodes correspond to the order of nodes. Nearest is an easy and efficient way to avoid making sure you get the index arithmetic right if you manually piece together the vertex/edge/centroid components.

  • crElements is a list whose elements represent Crouzeix-Raviart triangle elements. Each element has is a list of indices into nodes, {v1, v2, v3, e1, e2, e3, c1}, where v1, v2, v3 are the indices of the vertices, e1, e2, e3 are the indices of the midpoint of the edges, and c1 is the index of the centroid.

We can inspect the work so far, to make sure that the elements are in the right order:

Show[
 mreg,
 Graphics[{MapIndexed[Text[#2, #1] &, nodes ~Part~ First@crElements]}]
 ]

Mathematica graphics

One suggestion if the interpolation is going to be used at machine-precision is to use Compile to compute the phi's. It needs to be a function of the eight variables {x, y, x1, y1, x2, y2, x3, y3}:

phiv = {φ1, φ2, φ3, φ4, φ5, φ6, φ7};
xv = {x1, x2, x3};
yv = {y1, y2, y3};

Join[{x, y}, Riffle[xv, yv]]
(*  {x, y, x1, y1, x2, y2, x3, y3}  *)

phiC = With[{
    vars = Join[{x, y}, Riffle[xv, yv]],
    body = phiv},
   Compile @@ Hold[vars, body]
   ];

Here is the basic sequence for interpolation: Get the cell index; from that get the CR mesh element; get the data values corresponding to the nodes in that element; and finally compute the sum of the values times phi with a dot product.

Block[{x = -0.01, y = 0.01},
 cellI = Last@Region`Mesh`MeshMemberCellIndex[mreg, {x, y}];
 elem = crElements[[cellI]];
 fv = valsOnNodes[[elem]];
 fv.(phiC @@ Flatten[{x, y, nodes ~Part~ Take[elem, 3]}])
 ]
(*  0.989778  *)

We can put it in a function that has as a head the following data structure:

crInterpolatingFunction[
 mreg,    (* the mesh, a MeshRegion[] *)
 nodes,   (* all the nodes: the vertices, edge points, and centers *)
 vals,    (* the function values on the nodes in the same order *)
 elems    (* CR elements: each a list of indices of vertices, edge points, center *)
]

Then the interpolation is evaluated via a subvalue:

ClearAll[crInterpolatingFunction];
crInterpolatingFunction[
   mreg_MeshRegion, nodes_?ArrayQ, vals_?VectorQ, elems_?ArrayQ
   ][x_?NumericQ, y_?NumericQ] := With[
  {el = elems ~Part~ Last@Region`Mesh`MeshMemberCellIndex[mreg, {x, y}]},
  vals[[el]].(phiC @@ Flatten[{x, y, nodes ~Part~ Take[el, 3]}])
  ]

OP's example:

crIF = crInterpolatingFunction[mreg, nodes, valsOnNodes, crElements];

Show[
 Plot3D[crIF[x, y], {x, y} ∈ mreg],
 Graphics3D[{Darker@Red, PointSize@Large, Point@data}]
 ]

Mathematica graphics

I'm unfamiliar with C-R interpolation, so the only evaluation of the method I can make is to observe that it goes through the data points.

Update: Differentiable interpolating functions

The built-in InterpolatingFunctions are differentiated by keeping track of the derivative order internally, so we take the same approach below (see the argument dorder). To compute the derivative the derivative of the CR interpolating code is compiled (and memoized) when needed (see phiC).

The interface is slightly different than above since the full nodes are not needed in the CR interpolating code, just MeshCoordinates[mreg]. Some checking and extrapolation handling have also been added.

ClearAll[crInterpolatingFunction];
crInterpolatingFunction::dmval = InterpolatingFunction::dmval;
crInterpolatingFunction[
    mreg_MeshRegion, vals_?VectorQ, elems_?ArrayQ, dorder_: {0, 0}, extrap_: Automatic
    ][x_?NumericQ, y_?NumericQ] := Module[{cellI, el},
   cellI = Last@Region`Mesh`MeshMemberCellIndex[mreg, {x, y}];
   If[cellI == 0,
    Message[crInterpolatingFunction::dmval, {x, y}];
    If [extrap === Automatic,
     cellI = Last@Region`Mesh`MeshMemberCellIndex[mreg, RegionNearest[mreg, {x, y}]],
     Return[extrap@{x, y}, Module]
     ]
    ];
   el = elems[[cellI]];
   vals[[el]].(phiC[dorder] @@ Flatten[{x, y, MeshCoordinates[mreg] ~Part~ Take[el, 3]}])
   ];

phiv = {φ1, φ2, φ3, φ4, φ5, φ6, φ7};
xv = {x1, x2, x3};
yv = {y1, y2, y3};

Clear[phiC];
mem : phiC[dorder: {m_, n_}] := mem = Block[{x, y, x1, x2, x3, y1, y2, y3},
    With[{vars = Join[{x, y}, Riffle[xv, yv]],
      body = D[phiv, {x, m}, {y, n}]},
     Compile @@ Hold[vars, body]
     ]];

Derivative[m_, n_][
   crInterpolatingFunction[mreg_, vals_, elems_, dorder_: {0, 0}, extrap_: Automatic]] :=
  crInterpolatingFunction[mreg, vals, elems, dorder + {m, n}, extrap];

Examples:

crIF = crInterpolatingFunction[mreg, valsOnNodes, crElements];
D[crIF[x, y], y]
(*  crInterpolatingFunction[..., {0, 1}, Automatic][x, y]  *)

Show[
 Plot3D[Evaluate@{crIF[x, y], D[crIF[x, y], x]}, {x, y} ∈ mreg],
 Graphics3D[{Darker@Red, PointSize@Large, Point@data}],
 AxesLabel -> {x, y, z}
 ]
Improve this answer
$\endgroup$
10
  • $\begingroup$ I should point out that this is a minimal implementation. There is no checking if {x, y} is in the mesh region (cellI will be 0 in such a case), or extrapolation handling. I worked with the OP's setup, but I wonder if a more efficient interpolation routine phiC can be derived if all nodes are used, instead of writing it in terms of the vertices. I also noticed that nodes in crInterpolatingFunction can be omitted from the data structure and replaced with MeshCoordinates[mreg], which would save space for a large mesh. $\endgroup$ – Michael E2 Feb 10 '18 at 21:30
  • $\begingroup$ Hm, that valsonnodes line throws an error: The dimension of {{-0.45,-1.},...} and {-1.,0.967,-0.973,...} does not match. $\endgroup$ – user16320 Feb 10 '18 at 21:37
  • $\begingroup$ Well, the valsonnodes expression didn't change...still the same error for me. Isn't the problem with data[[All, {1,2}]] -> data[[All,3]]? Those two things have different dimensions.The first one is a list of doublets, the second is just a list of numbers. What is the purpose of this object? I Think it should be the list of values closest to nodal points (in this case, as the points are directly on nodes, it finds the appropriate value)...however the syntax might be erroneous. $\endgroup$ – user16320 Feb 10 '18 at 22:28
  • $\begingroup$ @user16320 I think before V10, you might have to use Map, something like: valsOnNodes = With[{nf = Nearest[data[[All, {1, 2}]] -> data[[All, 3]]]}, Flatten[nf /@ nodes, 1]];, but I can't check pre-V11.2. $\endgroup$ – Michael E2 Feb 10 '18 at 22:39
  • $\begingroup$ Oh yes, now it works...however, two things to consider: it is very slow (I tried on a bigger set of data). The second thing, more importantly, it doesn't support derivatives, if I do D[crIF[x, y], x], there is just a bunch of errors (mainly "The expression ... cannot be used as a part specification" and "Cannot take positions 1 through 3 in ..."). Is there a way to compute partial derivatives of such function? $\endgroup$ – user16320 Feb 10 '18 at 22:44

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.