FiniteDifferences.hpp

#ifndef __FIBER_FINITE_DIFFERENCING_HPP
#define __FIBER_FINITE_DIFFERENCING_HPP
 
#include <eagle/Vector.hpp>
#include "MultiArray.hpp"
#include "MultiOperate.hpp"
 
namespace Fiber
{
 
template <class ResultType>
struct  FiniteDifferenceTrait
{
        template <Dims_t  Dims, class DifferenceType>
static  void assign( MultiArray<Dims, ResultType>&result,
                     const MultiIndex<Dims>&Index,
                     int EvalDim,
                     const DifferenceType&DT)
        {
                result [ Index ][ EvalDim ] = DT;
        }
};
 
template <>
struct  FiniteDifferenceTrait<double>
{
        template <Dims_t Dims, class DifferenceType>
static  void assign( MultiArray<Dims, double>&result,
                     const MultiIndex<Dims>&Index,
                     int EvalDim,
                     const DifferenceType&DT)
        {
                result [ Index ] = DT;
        }
};
 
 
/*
template <class ValueType>
struct  FiniteDifferenceDifferentionTrait
{
static  auto difference(const ValueType&L, const ValueType&R)
        {
                return L - R;
        }
};
*/
 
 
#if     0
 
/* @ingroup marray
   Compute Finite Differences along one line within
   a multidimensional array in the given index.
   The evaluation dimension is given at runtime.
   Not sure if we need this.
   @note Internal use only.
 */
template <Dims_t Dims, typename Value, typename ResultType>
const MultiArray<Dims, ResultType>&
FiniteDifference( MultiArray<Dims, ResultType>&result,
                  const MultiArray<Dims, Value>&ValueField,
                  const MultiIndex<Dims>&Index,
                  int   EvalDimension)
{
const MultiIndex<Dims>&Sz = ValueField.Size();
 
        if (Sz[EvalDimension] < 2)
        {
//              throw something ?
                return result;
        }
 
MultiIndex<Dims> StartIndex = Index,
                MiddleIndex = Index,
                   EndIndex = Index;
 
        {
                StartIndex[ EvalDimension ] = 0;
                EndIndex  [ EvalDimension ] = 1;
 
                FiniteDifferenceTrait<ResultType>::assign( result,
                                                           StartIndex, EvalDimension,
                                                           ValueField[ EndIndex ] - ValueField[ StartIndex ] );
 
//              result[ StartIndex ][EvalDimension] = ValueField[ EndIndex ] - ValueField[ StartIndex ];
        }
 
        for(index_t i=1; i<Sz[EvalDimension]-1; i++)
        {
                StartIndex [ EvalDimension ] = i-1;
                MiddleIndex[ EvalDimension ] = i  ;
                EndIndex   [ EvalDimension ] = i+1;
 
                FiniteDifferenceTrait<ResultType>::assign( result,
                                                           MiddleIndex, EvalDimension,
                                                           0.5*(ValueField[ EndIndex ] - ValueField[ StartIndex ]) );
 
 
//              result[ MiddleIndex ][EvalDimension] = 0.5*(ValueField[ EndIndex ] - ValueField[ StartIndex ] );
        }
 
        {
                StartIndex[ EvalDimension ] = Sz[EvalDimension]-2;
                EndIndex  [ EvalDimension ] = Sz[EvalDimension]-1;
 
                FiniteDifferenceTrait<ResultType>::assign( result,
                                                           EndIndex, EvalDimension,
                                                           ValueField[ EndIndex ] - ValueField[ StartIndex ] );
 
//              result[ EndIndex ][EvalDimension] = ValueField[ EndIndex ] - ValueField[ StartIndex ];
        }
 
 
        return result;
}
 
 
template <Dims_t Dims, typename Value, typename ResultType>
struct  FiniteDifferenceOperation
{
        MultiArray<Dims, ResultType>*   result;
        const MultiArray<Dims, Value>*  ValueField;
        int                             EvalDimension;
 
        FiniteDifferenceOperation()
        {}
 
        void    apply(const MultiIndex<Dims>&Index)
        {
                FiniteDifference( *result, *ValueField, Index, EvalDimension);
        }
};
 
template <Dims_t Dims, typename Value, typename ResultType, int Diff>
const MultiArray<Dims, ResultType>&
FiniteDifference( MultiArray<Dims, ResultType>&result,
                  const MultiArray<Dims, Value>&ValueField, const MultiIndex<Diff>&)
{
typedef FiniteDifferenceOperation<Dims, Value, ResultType>
        FiniteDifferenceOperation_t;
 
typedef MultiOperate<Dims, FiniteDifferenceOperation_t, Diff> FD0_t;
 
FD0_t   FD0;
        FD0.result     = &result;
        FD0.ValueField = &ValueField;
        FD0.EvalDimension = Diff;
 
        FD0.iterate( ValueField.Size() );
 
        return result;
}
#endif
 
 
 
template <Dims_t EvalDimension>
struct  ComputePartialDerivative
{
 
        template <Dims_t Dims, typename Value, typename ResultType>
static  const MultiArray<Dims, ResultType>&
        FiniteDifference( MultiArray<Dims, ResultType>&result,
                          const MultiArray<Dims, Value>&ValueField,
                          const MultiIndex<Dims>&Index)
        {
        const MultiIndex<Dims>&Sz = ValueField.Size();
 
                if (Sz[EvalDimension] < 2)
                {
//              throw something ?
                        return result;
                }
 
        MultiIndex<Dims> StartIndex = Index,
                        MiddleIndex = Index,
                           EndIndex = Index;
 
                {
                        StartIndex[ EvalDimension ] = 0;
                        EndIndex  [ EvalDimension ] = 1;
 
                        FiniteDifferenceTrait<ResultType>::assign( result,
                                                                   StartIndex, EvalDimension,
                                                                   ValueField[ EndIndex ] - ValueField[ StartIndex ] );
 
        //              result[ StartIndex ][EvalDimension] = ValueField[ EndIndex ] - ValueField[ StartIndex ];
                }
 
                for(index_t i=1; i<Sz[EvalDimension]-1; i++)
                {
                        StartIndex [ EvalDimension ] = i-1;
                        MiddleIndex[ EvalDimension ] = i  ;
                        EndIndex   [ EvalDimension ] = i+1;
 
                        FiniteDifferenceTrait<ResultType>::assign( result,
                                                                   MiddleIndex, EvalDimension,
                                                                   0.5*(ValueField[ EndIndex ] - ValueField[ StartIndex ]) );
 
 
        //              result[ MiddleIndex ][EvalDimension] = 0.5*(ValueField[ EndIndex ] - ValueField[ StartIndex ] );
                }
 
                {
                        StartIndex[ EvalDimension ] = Sz[EvalDimension]-2;
                        EndIndex  [ EvalDimension ] = Sz[EvalDimension]-1;
 
                        FiniteDifferenceTrait<ResultType>::assign( result,
                                                                   EndIndex, EvalDimension,
                                                                   ValueField[ EndIndex ] - ValueField[ StartIndex ] );
 
        //              result[ EndIndex ][EvalDimension] = ValueField[ EndIndex ] - ValueField[ StartIndex ];
                }
                return result;
        }
 
 
template <Dims_t Dims, typename Value, typename ResultType>
        struct  FiniteDifferenceOperation
        {
                MultiArray<Dims, ResultType>*   result;
                const MultiArray<Dims, Value>*  ValueField;
 
                FiniteDifferenceOperation()
                {}
 
                void    apply(const MultiIndex<Dims>&Index)
                {
                        FiniteDifference( *result, *ValueField, Index);
                }
        };
 
 
template <Dims_t Dims, typename Value, typename ResultType>
static  MultiArray<Dims, ResultType>&
        FiniteDifference( MultiArray<Dims, ResultType>&result,
                          const MultiArray<Dims, Value>&ValueField)
        {
        typedef FiniteDifferenceOperation<Dims, Value, ResultType>
                FiniteDifferenceOperation_t;
 
        typedef MultiOperate<Dims, FiniteDifferenceOperation_t, EvalDimension> FD_t;
 
        FD_t    FD;
                FD.result     = &result;
                FD.ValueField = &ValueField;
 
                FD.iterate( ValueField.Size() );
 
                return result;
        }
 
template <Dims_t Dims, typename Value, typename ResultType>
static  MultiArray<Dims, ResultType>&
        compute( MultiArray<Dims, ResultType>&result,
                 const MultiArray<Dims, Value>&ValueField)
        {
                return ComputePartialDerivative<EvalDimension-1>::compute(
                        FiniteDifference( result, ValueField), ValueField
                                        );
        }
};
 
 
template <>
struct  ComputePartialDerivative<Dims_t(-1)>
{
template <Dims_t Dims, typename Value, typename ResultType>
static  MultiArray<Dims, ResultType>&
        compute( MultiArray<Dims, ResultType>&result,
                 const MultiArray<Dims, Value>&ValueField)
        {
                return result;
        }
 
};
 
 
template <Dims_t Dims, typename Value, typename ResultType>
const MultiArray<Dims, ResultType>&
        ComputeDerivative( MultiArray<Dims, ResultType>&result,
                           const MultiArray<Dims, Value>&ValueField)
{
        return ComputePartialDerivative<Dims-1>::compute( result, ValueField);
}
 
 
 
} /* namespace Fiber */
 
#endif /*  __FIBER_FINITE_DIFFERENCING_HPP */