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