derivatives.h

Go to the documentation of this file.

#pragma once
 
#include "grid.h"
#include "dx.h"
 
namespace dg{
 
    // Check derivativesT.h for more overloads
 
namespace create{
namespace detail
{
template<class real_type, size_t Nd>
void update_left_right( unsigned coord, EllSparseBlockMat<real_type, thrust::host_vector>& mat, const aRealTopology<real_type,Nd>& g)
{
    // Also used by fast_interpolation ...
    unsigned right_size = 1, left_size = 1;
    for( unsigned u=0; u<coord; u++)
        right_size*= g.shape(u);
    for( unsigned u=coord+1; u<Nd; u++)
        left_size *= g.shape(u);
    mat.set_right_size( right_size);
    mat.set_left_size( left_size);
}
}
 
 
template<class real_type, size_t Nd>
EllSparseBlockMat<real_type, thrust::host_vector> derivative( unsigned coord,
    const aRealTopology<real_type, Nd>& g, dg::bc bc, direction dir = centered)
{
    if( coord >= Nd)
        throw Error( Message(_ping_)<<"coord>=Nd not allowed! You typed: "<<coord<<" while Nd is "<<Nd);
    EllSparseBlockMat<real_type, thrust::host_vector> dd = dx_normed(
            g.n(coord), g.N(coord), g.h(coord), bc, dir);
    detail::update_left_right( coord, dd, g);
    return dd;
}
 
template<class real_type, size_t Nd>
EllSparseBlockMat<real_type, thrust::host_vector> jump( unsigned coord,
    const aRealTopology<real_type, Nd>& g, dg::bc bc)
{
    if( coord >= Nd)
        throw Error( Message(_ping_)<<"coord>=Nd not allowed! You typed: "<<coord<<" while Nd is "<<Nd);
    EllSparseBlockMat<real_type, thrust::host_vector> dd = jump(
            g.n(coord), g.N(coord), g.h(coord), bc);
    detail::update_left_right( coord, dd, g);
    return dd;
}
 
 
 
} //namespace create
 
} //namespace dg