PartialDerivative.hpp

#ifndef __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP
#define __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP

#include <eagle/Vector.hpp>
#include "MultiArray.hpp"

namespace Fiber
{


template <int IndexChart, typename Value, typename ResultType>
Eagle::Vector<ResultType, IndexChart>&
PartialDerivative(Eagle::Vector<ResultType, IndexChart>&result, const MultiArray<IndexChart, Value>&ValueField, const MultiIndex<IndexChart>&I)
{

const MultiIndex<IndexChart>&Sz = ValueField.Size();

        for(int i=0; i<IndexChart; i++)
        {
        MultiIndex<IndexChart> DiffI = MultiIndex<IndexChart>::Axis(i); 

                if (I[i]<1) // on left boundary, need to do forward-differencing
                {
                        result[i] = ValueField[ I + DiffI ] - ValueField[ I ]; 
                } 
                else if (I[i]>=Sz[i]-1) // on right boundary, need to do backward-differencing
                {
                        result[i] = ValueField[ I ] - ValueField[ I - DiffI ]; 
                }
                else
                {
                        result[i] =  0.5*(ValueField[ I + DiffI ] - ValueField[ I - DiffI ]);
                }
        }

        return result;
}


} /* namespace Fiber */ 

#endif /* __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP */