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 */