I have wrapped the interpolation method using the Fenics Project into a template class. The class is template on the interpolation basis function so that the user may choose the order of interpolation. The class also contains an array of y values so that it may serve as a table of values. Part of this is a map to the variable name for convenience when calling for tabled values.
Since all the work goes into extrapolating to a higher basis function when the Constructor is called, lookups are cheap to retrieve. The cost should be about the level of evaluating a polynomial of the order of the interpolation level (since the weights are known).
An example of it being used:
boost::shared_ptr<Vector> xs (new Vector(n));
boost::shared_ptr<Vector> ys1 (new Vector(n));
boost::shared_ptr<Vector> ys2 (new Vector(n));
std::vector<boost::shared_ptr<Vector> > yss;
std::vector<string> names;
// Populate test data
for(int i=0; i<n; i++)
{
xs->setitem(i,i+i*0.1);
ys1->setitem(i,(*xs)[i]*(*xs)[i]);
ys2->setitem(i,(*xs)[i]*(*xs)[i]*(*xs)[i]);
}
yss.push_back(ys1);
names.push_back("squared");
yss.push_back(ys2);
names.push_back("cubed");
InterpolationTable<Interpolate3::FunctionSpace> interpTable (xs, yss, names);
std::cout << interpTable.eval(0,2.12) << "\t" << interpTable.eval("squared",2.12) << "\n";
std::cout << interpTable.eval(1,2.12) << "\t" << interpTable.eval("cubed",2.12) << "\n";
The code listing:
// InterpolationTable.h
// Class for interpolating tabled values using the Fenics Project.
// Charles R. Cook v1.0 29 Nov 2011
// Copyright (c) 2011 Charles R. Cook
// http://www.charlesrcook.com
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#ifndef _INTERPOLATIONTABLE_H
#define _INTERPOLATIONTABLE_H
#include <dolfin.h>
#include "Interpolate1.h"
#include "Interpolate2.h"
#include "Interpolate3.h"
#include "Interpolate4.h"
using namespace dolfin;
template <class TFunctionSpace>
class InterpolationTable
{
public:
InterpolationTable(boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in);
InterpolationTable(boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in,
std::vector<std::string> names_in);
double eval(int index, double x);
double eval(std::string name, double x);
// These vectors contain the raw data
boost::shared_ptr<Vector> xs;
boost::shared_ptr<std::vector<boost::shared_ptr<Vector> > > ys;
// The mesh (defined as nodes at x locations)
boost::shared_ptr<Mesh> mesh;
boost::shared_ptr<Interpolate1::FunctionSpace> V1;
boost::shared_ptr<TFunctionSpace> V2;
boost::shared_ptr<std::vector<Function> > fy1;
boost::shared_ptr<std::vector<Function> > fy2;
boost::shared_ptr<std::vector<std::string> > Names;
private:
void Setup(boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in);
boost::shared_ptr<std::map<std::string,int> > mapping;
};
template <class TFunctionSpace>
double InterpolationTable<TFunctionSpace>::eval(
std::string name, double x)
{
if((*Names).size()==0)
throw ("names not defined");
// check key exists in map
if ((*mapping).find(name) == (*mapping).end() )
throw ("name was not found");
int index = (*mapping)[name];
return eval(index, x);
}
template <class TFunctionSpace>
double InterpolationTable<TFunctionSpace>::eval(
int index, double x)
{
if(index >= (*ys).size())
throw ("index out of range");
if(x < mesh->coordinates()[0] ||
x > mesh->coordinates()[(*mesh).num_cells()])
throw ("x location out of range");
return (*fy2)[index](x);
}
template <class TFunctionSpace>
void InterpolationTable<TFunctionSpace>::Setup(
boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in)
{
// copy the x locations in
xs = xs_in;
// copy the data in
for (uint i=0; i<ys_in.size(); i++)
{
(*ys)[i] = ys_in[i];
}
// Setup the mesh to interpolate on
// Note the count is the number of elements, not nodes
// when the interval is constructed
for(uint i=0; i<xs->size(); i++)
mesh->coordinates()[i] = (*xs)[i];
// Setup the Function space(s)
for (uint i=0; i<ys->size(); i++)
{
Function fv1(V1,(*ys)[i]);
fy1->push_back(fv1);
Function fv2(V2);
fv2.extrapolate(fv1);
fy2->push_back(fv2);
}
}
template <class TFunctionSpace>
InterpolationTable<TFunctionSpace>::InterpolationTable(boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in,
std::vector<std::string> names_in) :
xs(new Vector(xs_in->size())),
ys(new std::vector<boost::shared_ptr<Vector> > (ys_in.size())),
mesh(new Interval(xs_in->size()-1,0,1)),
V1(new Interpolate1::FunctionSpace(mesh)),
V2(new TFunctionSpace(mesh)),
fy1(new std::vector<Function>),
fy2(new std::vector<Function>),
Names(new std::vector<std::string>),
mapping(new std::map<std::string, int>)
{
Setup(xs_in, ys_in);
(*Names) = names_in;
for(uint i=0; i < (*Names).size(); i++)
(*mapping)[(*Names)[i]] = i;
}
template <class TFunctionSpace>
InterpolationTable<TFunctionSpace>::InterpolationTable(boost::shared_ptr<Vector> xs_in,
std::vector<boost::shared_ptr<Vector> > ys_in) :
xs(new Vector(xs_in->size())),
ys(new std::vector<boost::shared_ptr<Vector> > (ys_in.size())),
mesh(new Interval(xs_in->size()-1,0,1)),
V1(new Interpolate1::FunctionSpace(mesh)),
V2(new TFunctionSpace(mesh)),
fy1(new std::vector<Function>),
fy2(new std::vector<Function>),
Names(new std::vector<std::string>),
mapping(new std::map<std::string, int>)
{
Setup(xs_in, ys_in);
}
#endif /* _INTERPOLATIONTABLE_H */