I'll give a more general answer independent of implementation.
TL:DR
The point $x$ you seek the value $u$ at lies in a simplex with values at the nodes $u_i$. $x$ has the barycentric coordinates $\xi_i^0$, then:
$u = \sum_{i=1}^{n+1} u_i \xi^0_i$
is the linearly interpolated value at $x$. For a segment, take $\xi_1 = (1-t)$ and $\xi_2 = t$ s.t. $x = (1-t)x_1 + tx_2$ where $x_1,x_2$ are the array bounds.
Recipe for an interpolation scheme
First, the general topic is interpolation: you know a function's values $u_i$ at the points $\{x_i\}_i$ and you want to generate a value at a point $x$ "inside" of the $x_i$. What does inside mean? It means we need a notion of connectivity between the points:
- If $x_i \in \mathbb R$, you can order them as $x_1 \leq x_2 \leq ... \leq x_N$. We can connect the points into segments $[x_i,x_{i+1}]$.
- If $x_i \in \mathbb R^n$, presumably you have some connectivity information such as a mesh. If you don't, you can generate a Delaunay Triangulation of the $\{x_i\}_i$.
They differ from a practical standpoint, but the two cases are identical mathematically: segments, triangles, tetrahedra... are all defined as the convex hull of $2$, $3$, $4$ points.
These "elements" are thus all described by barycentric coordinates. Indeed, since they are convex hulls, they are the set of convex combinations of their vertices. A convex combination is one where the coefficients $\xi_i$ verify $0\leq \xi_i \leq 1$ and $\sum \xi_i = 1$. This combination is also unique. Thus it's very convenient, for a point $x$ within an element $K$, to consider instead its barycentric coordinates $\xi$ in $K$.
Let us now denote $x$ the point we seek a value $u$ at; $x$ lies in a simplex $K$ of vertices $x_1, ..., x_{n+1}$ and has the barycentric coordinates $\xi_i^0, ..., \xi_{n+1}^0$.
There are many ways to assign a value $u$ to $x$.
The general method is to :
- Choose a finite-dimensional family of functions: polynomials of given degree, other polynomial families, solutions to a given PDE indexed by boundary conditions... there is an infinity of possible choices here, not all good or practical. These functions may either be defined over the set of barycentric coordinates, or over the set of physical coordinates.
- Fit the family to the values: that is find the function(s) that best evaluates to the $u_{i}$ at the $x_{i}$ or at the $\xi^0_{i}$
- Evaluate at $x$
The distinction between the $x_i$ and the $\xi_i$ may seem pedantic, but it becomes important if you are not working on a simplex but on a quadrilateral, a rather common element in meshes. Non-parallelogram quadrilaterals are not linear elements (they are quadratic), thus polynomials of the $\xi_i$ are not polynomials of the $x_i$, and polynomials of degree $d$ of the $x_i$ are polynomials of degree $2d$ of the $\xi_i$.
The fitting step may be as simple as solving a set of equations, or it may involve minimization of a chosen norm (like an $L^2$ norm). In some cases, you may have to add additional constraints to ensure a unique solution. For instance, you may want the solution with minimal $H^1$ semi-norm to attenuate oscillations. Anyways, let's take the simplest choice for a simplex:
- Interpolate on degree $1$ polynomials. They either write $f : x \mapsto a_0 + a^Tx$ with $a$ has $n$ entries in physical coordinates or $f:\xi \mapsto a^Tx$ with $a$ has $n+1$ entries in barycentric coordinates.
- We have as many linear equations $f(x_i) = u_i$ as unknowns to define $f$ so to "fit" is simply to solve some linear equations.
- The transformation from $x$ to $\xi$ is linear so we can either consider linear functions of $x$ or of $\xi$ (remain linear under the mapping). I'll use the barycentric coordinates $\xi$.
The barycentric coordinate of the vertex $x_1$ is simply $(1, 0, 0...)$, of the vertex $x_2$, $(0,1,0....0)$ and so on. Indeed, the unique combination $\sum \xi_i x_i = x_1$ is with $\xi_1 = 1$ and all other $\xi_i = 0$. Thus our fitting $f^*$ verifies very simply:
$f^*((1,0,....,0)) = a^T (1,0,...,0) = a_1 = u_1\\
f^*((0,1,0,....,0)) = a^T (0,1,0,...,0) = a_2 = u_2\\ ...$
So our fitting function is very simply
$f^*(\xi) = \sum_{i=1}^{n+1} u_i \xi_i$
and when we evaluate it at $x = \sum \xi^0_i x_i$, we get the interpolated value $u$:
$u = \sum_{i=1}^{n+1} u_i \xi^0_i$.
Generalizations
If the goal had been to arrive at that measly linear combination, all this exposition would have been unnecessary. But I feel it's important to get have a grasp of how to construct your own interpolation scheme for the problems you encounter.
What if my elements are not linear ?
You can consider simplices with non-linear mappings. This requires thinking of elements not as a set of points, but as a mapping defined over the set of barycentric coordinates. For a linear segment $K$ with ends $x_1, x_2$, this is simply:
$F_1 : (\xi_1, \xi_2) \mapsto \xi_1 x_1 + \xi_2 x_2$
A quadratic segment can be describe, for instance, in Bézier form. Introduce $x_{1/2}$ a new point (inside of $[x_1,x_2]$), then
$F_K : (\xi_1,\xi_2) \mapsto \xi_1^2 \xi_1 + 2\xi_1\xi_2 x_{1/2} \xi_2^2 x_2$
is the mapping of the Bézier segment with $x_{1/2}$ as the control point. Under some conditions (invertibility), any point inside of $[x_1,x_2]$ still has unique barycentric coordinates in the element $K$ with mapping $F_K$, and they are generally not the same as if you considered $F_1$ instead. From this point on, the procedure holds: choose a family of functions (whether of $\xi$ or $x$ is another important choice), solve linear equations, evaluate.
What if I don't have elements at all ?
Perhaps you lack connectivity, or you want to interpolate with respect to more points (say $M$) than make up the vertices of a simplex. In that case, you have no choice but to carry out physical interpolation (no barycentric coordinates; technically you could probably consider polytopes in some cases). You still have to choose a family of functions, say $F$ of dimension $N$. With $M\ne N$ in the general case, you cannot expect $f(x_i) = u_i$; this is where fitting comes in: you instead minimize a norm, such as $\sum |f(x_i) - u_i|^2$. You may also incorporate constraints if you have many functions and not many points ($N>M$). Check out least-squares regression. This is not set in stone: no-one prevents you from minimizing $\sum |f(x_i) - u_i|^2 + \sum ||\nabla f(x_i)||^2$ to tone down oscillations.
Some applications
Galerkin schemes (as opposed to Finite Differences) fundamentally define interpolation bases, and then carry out with their analysis.
To give concrete examples, Isoparametric Finite Element methods consider interpolation in the barycentric coordinates by functions of the same degree as the element mappings, Iso-Geometric Analysis considers interpolation in barycentric coordinates by rational functions (the elements are themselves described by rational mappings), physical-frame high-order schemes consider interpolation in the physical coordinates by functions of the same degree as element mappings. On polyhedra, harmonic "barycentric" coordinates are obtained by solving $\Delta \xi = f$ with particular boundary conditions, and then carry on as above. Trefftz schemes consider barycentric coordinates but families of functions that are themselves solutions to simple forms of the problem at hand (PDE). Etc... etc... etc...
