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>
 
 
#define USE_DEPRECATED_LINESET_FUNCTIONS
#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;
        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 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 ) )
                {
                        if (!C.fragV )
                                puts("CurveDerivative::CurveDerivative() No Coords for global fragment!");
                        else
                                printf("CurveDerivative::CurveDerivative() No Coords for fragment [%s]!\n", C.fragV->Name().c_str() );
 
                        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