CurveDerivative.hpp

#ifndef __FIBER_GRID_TYPES_CURVEDERIVATIVE_HPP
#define __FIBER_GRID_TYPES_CURVEDERIVATIVE_HPP

#include "gridtypesDllApi.h"
#include <grid/Grid.hpp>
#include <eagle/PhysicalSpace.hpp>

#include <grid/EuclideanMetric.hpp>
#include "LineSet.hpp"

//#define VERBOSE

namespace Fiber
{

struct  IdentityOperator
{
        template <class ValueType, class PointType>
 static const ValueType&compute(const ValueType&Value, index_t, const PointType&)
        {
                return Value;
        }
};

struct  ComputeUnitVector
{
        template <class ValueType, class PointType>
 static const ValueType&compute(const ValueType&Value, index_t, const PointType&)
        {
                return unit(Value);
        }
};


struct DerivativeOperatorData
                {
                        WeakPtr<Grid>             theGrid;
                        const RefPtr<FragmentID>& fragV;
                        string                    fieldName;
                        
                        DerivativeOperatorData(WeakPtr<Grid> theGridP, const string& fieldNameP, const RefPtr<FragmentID>& fragVP )
                        : theGrid( theGridP )
                        , fragV( fragVP )
                        , fieldName(fieldNameP)
                        {}
                };

template <class Type, class Metric = EuclideanMetric, class ValueOperator = IdentityOperator>
struct  CurveDerivative : public Metric, public ValueOperator
{
        typedef typename Metric::point  point;
        using Metric::distance;

        using ValueOperator::compute;

        RefPtr<MemArray<1, Type> >      Quantity;
        RefPtr<MemArray<1, Type> >      DerivedQuantity;

        typedef MemArray<1, Type>                       ResultArray_t;
//      typedef std::pair<WeakPtr<Grid>, string>        Constructor_t;
        typedef DerivativeOperatorData                  Constructor_t;

        RefPtr<ResultArray_t> result() const
        {
                return DerivedQuantity;
        }

        typedef std::string Parameters;

        CurveDerivative(const Constructor_t&C, const MemBase::Creator_t&Crec )
        {
#ifdef VERBOSE
                puts("CurveDerivative::CurveDerivative() enter");fflush(stdout);
#endif
                if (!C.theGrid)
                {
                        puts("CurveDerivative::CurveDerivative() No Grid!");fflush(stdout);
                        return; 
                } 

                LineSet LS = C.theGrid; 

                if (!LS)
                {
                        puts("CurveDerivative::CurveDerivative() No LineSet!");fflush(stdout);
                        return; 
                } 

                        
                if (!LS.CartesianVertices)
                {
                        puts("CurveDerivative::CurveDerivative() No Cartesian Vertices!");fflush(stdout);
                        return; 
                }

                if (!LS.getCoords( C.fragV ) )
                {
                        puts("CurveDerivative::CurveDerivative() No Coords!");fflush(stdout);
                        return;
                } 

#ifdef VERBOSE
                puts("CurveDerivative::CurveDerivative() progressing");fflush(stdout);
#endif

        string  QuantityFieldName = C.fieldName; 

#ifdef VERBOSE 
        puts("CurveDerivative::CurveDerivative() progressing 2");fflush(stdout);
        std::cout << "CurveDerivative::CurveDerivative() doing: " << QuantityFieldName << std::endl;fflush(stdout);
#endif          

        RefPtr<Field> QuantityVertexField = (*LS.CartesianVertices)[ QuantityFieldName ]; 
                if (!QuantityVertexField)
                {
                        puts("CurveDerivative::CurveDerivative() No Quantity Field found!");fflush(stdout);
                        return; 
                }
                if (Quantity = QuantityVertexField->getData( C.fragV ) )
                {
                        DerivedQuantity = new ResultArray_t( LS.getCoords( C.fragV )->nElements(), Crec ); 

#ifdef VERBOSE 
                puts("CurveDerivative::CurveDerivative() iterating"); fflush(stdout);
#endif 
                LS.ApplyOnVertexFragment( *this,  C.fragV );

#ifdef VERBOSE
                puts("CurveDerivative::CurveDerivative() post iteration"); fflush(stdout);
#endif

                } 
                else
                {
                        puts("CurveDerivative::CurveDerivative() No Data retrieved!");fflush(stdout);
                }
        }

        void    apply(const CreativeIterator<Eagle::PhysicalSpace::point>&Vertices,
                      const LineSet::LineIndices_t&Line)
        {
#ifdef VERBOSE
                puts("CurveDerivative::apply()");fflush(stdout); 
#endif
        MultiArray<1, Type>&Values          = *Quantity;
        MultiArray<1, Type>&ValueDerivative = *DerivedQuantity; 

                if (Line.size()<1) // no data
                {
#ifdef VERBOSE
                        puts("CurveDerivative::apply() No Data"); 
#endif
                        return; 
                }
                if (Line.size()<2) // one element - derivative is zero
                {
#ifdef VERBOSE
                        puts("CurveDerivative::apply() One Element"); 
#endif
                const point&P0 = Vertices[ Line[0] ];
                const Type&Q = ValueOperator::compute( Values  [ Line[0] ], Line[0], P0 ); 

                        ValueDerivative[ Line[0] ] = Q-Q;
#ifdef VERBOSE 
                        std::cout << "First" << ValueDerivative[ Line[0] ] << std::endl;
#endif
                        return;
                } 
                // at least two elements, 0,1, here

                // First point: forward difference
                {
                const point&P0 = Vertices[ Line[0] ],
                           &P1 = Vertices[ Line[1] ]; 

/*
                const Type&Q0 = ValueOperator::compute( Values  [ Line[0] ], Line[0], P0 );
                const Type&Q1 = ValueOperator::compute( Values  [ Line[1] ], Line[1], P1 );

                double  ds = Metric::distance(Line[0], Line[1], P0, P1);

                        ValueDerivative[ Line[0] ] = (Q1 - Q0) / ds;
*/ 

                const Type&Q0 = Values  [ Line[0] ]; 
                const Type&Q1 = Values  [ Line[1] ]; 

                double  ds = Eagle::norm(P1 - P0);

                        ValueDerivative[ Line[0] ] = (Q1 - Q0) / ds; 
#ifdef VERBOSE 
                        std::cout << "Q0= " << Q0 << ", Q1= " << Q1 << std::endl;
                        std::cout << "First " << ValueDerivative[ Line[0] ] << std::endl;
#endif  
                } 

                // 
                // TODO: Reuse already computed values and point coordinates. 
                //

        index_t I = Line.size();
                for(unsigned i=1; i<I-1; i++)
                {
                const point&Pm = Vertices[ Line[i-1] ],
                           &P  = Vertices[ Line[i  ] ],
                           &Pp = Vertices[ Line[i+1] ]; 

/*              const Type&Qm = ValueOperator::compute( Values  [ Line[i-1] ], Line[i-1], Pm ),
                           Qp = ValueOperator::compute( Values  [ Line[i+1] ], Line[i+1], Pp ); 

                double  ds_m = Metric::distance(Line[i-1], Line[i  ], Pm, P),
                        ds_p = Metric::distance(Line[i  ], Line[i+1], P , Pp);

                        ValueDerivative[ Line[i] ] = (Qp - Qm) / ( ds_m + ds_p);
*/ 

                const Type&Qm = Values  [ Line[i-1] ],
                          &Q  = Values  [ Line[i  ] ],
                          &Qp = Values  [ Line[i+1] ]; 

                double  ds2_m = 2*Eagle::norm( P - Pm ),
                        ds2_p = 2*Eagle::norm( Pp - P ); 

                        ValueDerivative[ Line[i] ] =  (Q - Qm ) / ds2_m
                                                      +
                                                      (Qp - Q ) / ds2_p
                                                      ; 

#ifdef VERBOSE 
                std::cout << "i: " << i << " " << ValueDerivative[ Line[i] ] << std::endl;
#endif                          
                }

                // Last point: backward difference 
                if (I>2)
                {
                index_t I = Line.size();

                const point&P0 = Vertices[ Line[I-2] ],
                           &P1 = Vertices[ Line[I-1] ]; 

/*              const Type&Q0 = ValueOperator::compute( Values  [ Line[I-2] ], Line[I-2], P0 );
                const Type&Q1 = ValueOperator::compute( Values  [ Line[I-1] ], Line[I-1], P1 );

                double  ds = Metric::distance(Line[I-2], Line[I-1], P0, P1);
*/ 

                const Type&Q0 = Values  [ Line[I-2] ]; 
                const Type&Q1 = Values  [ Line[I-1] ]; 

                double  ds = Eagle::norm(P1-P0);

                        ValueDerivative[ Line[I-1] ] = (Q1 - Q0) / ds; 
#ifdef VERBOSE 
                std::cout << "Last: " << ValueDerivative[ Line[I-1] ] << std::endl;
#endif
                }

        }
};

} // namespace Fiber

#endif // __FIBER_GRID_TYPES_CURVEDERIVATIVE_HPP