My recent research has been with the Fenics Project, which is an amazing finite element project. In general it provides the tools needed to solve differential equations with the finite element method. For some examples check out their website and view the demos and applications (my interest is in CFD).
In my particular code I need to interpolate tabulated properties (steam tables) with some reasonable level of accuracy. The solution in MATLAB is straight forward, call interp1:
y_i = interp1(xs, ys, x_i, 'cubic');
Where xs and ys are vectors of known x and y values and y_i is the interpolated value at x_i.
The code I am working on is already using the Fenics Project to solve a set of PDEs, so I wanted to see if I could use it for the interpolation as well through the finite element's basis functions. As the title of this post indicates, it can be done and here is how through a code snippet.
In C++…
// Setup data vectors
Vector xs(n);
Vector ys(n);
// Populate test data (this would be the tabulated data)
for(int i=0; i<n; i++)
{
xs.setitem(i,i+i*0.1);
ys.setitem(i,xs[i]*xs[i]);
}
// Setup the mesh to interpolate on
// Note the count is the number of elements, not nodes.
// there are general n-1 elements than nodes.
Interval mesh (n-1,0,1);
for(int i=0; i<n; i++)
mesh.coordinates()[i] = xs[i];
// create the function space of order one.
// this is done so that dof maps directly to
// data points. This provides a linear interpolation.
Interpolate1::FunctionSpace V (mesh);
// copy the data into the function space
Function fy (V, ys);
// Create quadratic function space (quadratic interpolation)
Interpolate2::FunctionSpace V2 (mesh);
Function fy2 (V2);
// 'fit' to the quadratic basis functions through least squares
fy2.extrapolate(fy);
// The fitted function is now available in fy2 as a quadratic fit
std::cout << "O(1): " << fy(2)
<< "\tO(2): " << fy2(2)
<< "\n";
Note that two UFL files (form files) are needed for this snippet.
Interpolate1.ufl
# First order (linear) lagrange elements (polynomials)
element = FiniteElement("Lagrange", interval, 1)
Interpolate2.ufl
# Second order lagrange elements (polynomials)
element = FiniteElement("Lagrange", interval, 2)
To create higher order fitting you can change the element to be of the desired order (instead of 2).
The very nice feature of interpolation through this method is that it will work in two dimensions and three dimensions as well. To do so, the elements would be updated to 2-d or 3-d element types instead of 'interval' (say triangle or tetrahedral). Too cool!