posts - 81, comments - 262, trackbacks - 0

A C++ Interpolation Class for Tabled Data


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 */

 


Print | posted on Tuesday, November 29, 2011 12:27 PM | Filed Under [ Technical Computing ]

Powered by:
Powered By Subtext Powered By ASP.NET