FieldInterpolator.hpp

#ifndef _FieldInterpolator_HPP
#define _FieldInterpolator_HPP
 
#include <fish/fiber/vector/Interpolate.hpp>
#include <fish/fiber/vector/LinearIpol.hpp>
#include <fish/fiber/vector/CubicIpol.hpp>
#include <aerie/AnalyticFunction.hpp>
#include <fish/fiber/baseop/RectilinearInterpolation.hpp>
#include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
 
namespace Fiber
{
typedef MemArray<3, Eagle::tvector3> VectorArray_t;
class   FieldInterpolatorBase
{
protected:
        RefPtr<Field> TheField;
 
public:
        // supported interpolation types, used as template parameter in the FielInterpolator class
        enum InterpolType{LINEAR, CUBIC, RECTILINEAR, ANALYTIC, EIGENLINEAR, EIGENCUBIC};
 
        FieldInterpolatorBase( RefPtr<Field>& TheFieldP )
        :TheField(TheFieldP)
        {}
 
//      typedef pair<Eagle::PhysicalSpace::point, RefPtr<FragmentID> > FragmentedPoint_t;
        typedef LocalPoint FragmentedPoint_t;
 
        // convenience function
        RefPtr< MemBase > getArray(const FragmentedPoint_t& localPoint, const double time )
        {
//      const string&FragName = localPoint.second; 
                //const Eagle::PhysicalSpace::point &float_index = localPoint.first; 
 
//      RefPtr<FragmentID> FragID; 
//
//              if(FragName.compare("") != 0 )
//                      FragID = new FragmentID( FragName ); 
                
                return TheField->getCreator( localPoint.frag_id )->create();
        }
};
 
template<typename FieldType,  typename CoordType, int InterpolType>
class FieldInterpolator : public FieldInterpolatorBase
{
public:
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP)
        {}
 
        FieldType doIt( const  FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
                puts("template<FieldType>::Interpolator::doIt() no valid Interpolation type specialization"); 
 
        FieldType data; 
                return data;
        }  
};
 
template<typename FieldType, typename CoordType>
class FieldInterpolator<FieldType, CoordType, FieldInterpolatorBase::LINEAR> : public FieldInterpolatorBase
{
        typedef Interpolate<CoordType::Dims, FieldType, LinearIpol< FieldType> > LinearInterpolator_t;
public:
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) {}
 
        FieldType doIt( const  FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        RefPtr<MemArray<CoordType::Dims, FieldType> > ToIntArr = getArray(localPoint, time); 
        assert( ToIntArr );
 
        FieldType data; 
        LinearInterpolator_t MyInterpolator(*ToIntArr, localPoint.local_point); 
 
                data = MyInterpolator.eval(); 
 
                return data;
        }
};
 
template<typename FieldType, typename CoordType>
class FieldInterpolator<FieldType, CoordType, FieldInterpolatorBase::RECTILINEAR> : public FieldInterpolatorBase
{
public:
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) {}
 
        FieldType doIt( const  FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        assert(CoordType::Dims == 3);
        RefPtr<MemArray<3, FieldType> > ToIntArr = getArray(localPoint, time); 
        assert( ToIntArr );
        
        FieldType data;
        RectilinearInterpolation::rectilinear_interpolation(localPoint.local_point, TheField, data);
        return data;
        }
};
 
template<typename FieldType, typename CoordType>
class FieldInterpolator<FieldType, CoordType, FieldInterpolatorBase::CUBIC> : public FieldInterpolatorBase
{
        typedef Interpolate<CoordType::Dims, FieldType, CubicIpol< FieldType> > CubicInterpolator_t;
public:
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) {}
 
        FieldType doIt( const FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        RefPtr<MemArray<CoordType::Dims, FieldType> > ToIntArr =getArray(localPoint, time); 
        assert( ToIntArr ); 
 
        FieldType data; 
 
        CubicInterpolator_t MyInterpolator(*ToIntArr, localPoint.local_point); 
                data = MyInterpolator.eval(); 
 
                return data;
        }       
};
 
template<typename FieldType, typename CoordType>
class FieldInterpolator<FieldType, CoordType, FieldInterpolatorBase::ANALYTIC> : public FieldInterpolatorBase
{
typedef Eagle::PhysicalSpace::AnalyticFunction<FieldType> Eqn_t;
 
public:
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) {}
 
        FieldType doIt( const  FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        FieldType data; 
 
                if ( MemCore::InterfacePtr<Eqn_t> AnalyticField = *TheField )
                {
                Eqn_t&AF = *AnalyticField; 
                        data = AF(time, worldlocation);
                } 
                else
                {
                        puts("Interpolator<FieldType, Dims, InterpolatorBase::ANALYTIC>::doIt() ERROR: No analytic equation!");
                } 
 
                return data;
        }
};
 
 
class EigenVectorFlipper
{
public:
        tvector unique_orientation;
 
        EigenVectorFlipper( tvector& unique_orientationP )
        : unique_orientation ( unique_orientationP )
        {}
 
        tvector limit( const tvector& T) const
        {
                //std:: cout << "EigenVectorFlipper::limit() orientation: " << unique_orientation << std::endl;
 
                if( dot( T, unique_orientation ) < 0)
                        return -T;
 
                return T;
        }
};
 
template<typename CoordType>
class FieldInterpolator<tvector, CoordType, FieldInterpolatorBase::EIGENLINEAR> : public FieldInterpolatorBase
{
        typedef Interpolate<CoordType::Dims, tvector, LinearIpol< tvector>, double, EigenVectorFlipper > LinearInterpolator_t;
 
        tvector unique_orientation;
 
public:
 
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) 
        , unique_orientation(1.0, 0.0, 0.0)
{}
 
        void setOrientation(tvector& unique_orientationP )
        {
                unique_orientation =  unique_orientationP;      
        }
 
        tvector doIt( const FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        RefPtr<MemArray<CoordType::Dims, tvector> > ToIntArr = getArray(localPoint, time); 
        assert( ToIntArr ); 
 
        EigenVectorFlipper EigenDelimiter( unique_orientation );
 
        tvector data; 
        LinearInterpolator_t MyInterpolator(*ToIntArr, localPoint.local_point, EigenDelimiter); 
 
                data = MyInterpolator.eval(); 
 
                return data;
        }
};
 
template<typename CoordType>
class FieldInterpolator<tvector, CoordType, FieldInterpolatorBase::EIGENCUBIC> : public FieldInterpolatorBase
{
        typedef Interpolate<CoordType::Dims, tvector, CubicIpol< tvector>, double, EigenVectorFlipper > LinearInterpolator_t;
 
        tvector unique_orientation;
 
public:
 
        FieldInterpolator( RefPtr<Field>& TheFieldP )
        :FieldInterpolatorBase(TheFieldP) 
        , unique_orientation(1.0, 0.0, 0.0)
{}
 
        void setOrientation(tvector& unique_orientationP )
        {
                unique_orientation =  unique_orientationP;      
        }
 
        tvector doIt( const FragmentedPoint_t& localPoint, const CoordType& worldlocation, const double time )
        {
        RefPtr<MemArray<CoordType::Dims, tvector> > ToIntArr = getArray(localPoint, time); 
        assert( ToIntArr ); 
 
        EigenVectorFlipper EigenDelimiter( unique_orientation );
 
        tvector data; 
        LinearInterpolator_t MyInterpolator(*ToIntArr, localPoint.local_point, EigenDelimiter); 
 
                data = MyInterpolator.eval(); 
 
                return data;
        }
};
 
} // namespace Fiber
 
#endif //_FieldInterpolator_HPP