core/derivative.hh

#ifndef VTESTBED_CORE_DERIVATIVE_HH
#define VTESTBED_CORE_DERIVATIVE_HH

#include <vtestbed/core/types.hh>

namespace vtb {

    namespace core {

        template <class T, int N>
        void
        derivative(
            Array<T,N> func,
            const T delta,
            int dimension,
            Array<T,N>& result
        );

        template <class T, int N>
        inline Array<T,N>
        derivative(Array<T, N> func, const T delta, int dimension) {
            Array<T,N> result{func.shape()};
            ::vtb::core::derivative<T,N>(func, delta, dimension, result);
            return result;
        }

        template <class T, int N>
        T
        derivative(
            const blitz::TinyVector<int,N>& idx,
            int idx_min,
            int idx_max,
            int dimension,
            const T delta,
            const Array<T,N>& func
        ) {
            typedef blitz::TinyVector<int,N> intN;
            const int idx_d = idx(dimension);
            T result;
            if (idx_d == idx_min) {
                intN idx1(idx);
                ++idx1(dimension);
                intN idx2(idx1);
                ++idx2(dimension);
                result = T(0.5)
                    * (-func(idx2) + T(4)*func(idx1) - T(3)*func(idx))
                    / delta;
            } else if (idx_d == idx_max) {
                intN idx1(idx);
                --idx1(dimension);
                intN idx2(idx1);
                --idx2(dimension);
                result = T(0.5)
                    * (T(3)*func(idx) - T(4)*func(idx1) + func(idx2))
                    / delta;
            } else {
                intN idx0(idx);
                --idx0(dimension);
                intN idx1(idx);
                ++idx1(dimension);
                result = T(0.5)
                    * (func(idx1) - func(idx0))
                    / delta;
            }
            return result;
        }

        template <class T, int N>
        T
        derivative(
            const blitz::TinyVector<int,N>& idx,
            int dimension,
            const T delta,
            const Array<T,N>& func
        ) {
            return ::vtb::core::derivative<T,N>(
                idx, 0, func.extent(dimension)-1, dimension, delta, func
            );
        }

        template <class T, int N>
        blitz::TinyVector<T,N>
        gradient(
            const blitz::TinyVector<int,N>& idx,
            const blitz::TinyVector<T,N>& delta,
            const Array<T,N>& func
        ) {
            blitz::TinyVector<T,N> result;
            for (int i=0; i<N; ++i) {
                result(i) = derivative(idx, i, delta(i), func);
            }
            return result;
        }

        template <class T>
        blitz::TinyVector<T,4>
        gradient(
            const blitz::TinyVector<int,4>& idx,
            const blitz::TinyVector<T,4>& delta,
            const Array<T,4>& func
        ) {
            blitz::TinyVector<T,4> result;
            result(0) = derivative(idx, 0, delta(0), func);
            result(1) = derivative(idx, 1, delta(1), func);
            result(2) = derivative(idx, 2, delta(2), func);
            result(3) = derivative(idx, 3, delta(3), func);
            return result;
        }

    }

}

#endif // vim:filetype=cpp