transform.h

Go to the documentation of this file.

#pragma once
#include "multiply.h"
#include "base_geometry.h"
#include "weights.h"
 
 
namespace dg
{
namespace detail{
template< class Vector, class Functor, class RecursiveVector, size_t ...I>
auto do_pullback( Vector& result, Functor f, const RecursiveVector& map,
        std::index_sequence<I...>)
{
    return dg::blas1::evaluate( result, dg::equals(), f, map[I]...);
}
} //namespace detail
 
template< class Functor, class Geometry>
typename Geometry::host_vector pullback(
        const Functor& f, const Geometry& g)
{
    const std::vector<typename Geometry::host_vector >& map = g.map();
    typename Geometry::host_vector result( map[0]);
    detail::do_pullback( result, f, map,
            std::make_index_sequence<Geometry::ndim()>());
    return result;
}
 
template<class Functor1, class Functor2, class container, class Geometry>
void pushForwardPerp( const Functor1& vR, const Functor2& vZ,
        container& vx, container& vy,
        const Geometry& g)
{
    using host_vec = typename Geometry::host_vector;
    host_vec out1 = pullback( vR, g);
    host_vec out2 = pullback( vZ, g);
    dg::tensor::multiply2d(g.jacobian(), out1, out2, out1, out2);
    dg::assign( out1, vx);
    dg::assign( out2, vy);
}
 
template<class Functor1, class Functor2, class Functor3, class container, class Geometry>
void pushForward( const Functor1& vR, const Functor2& vZ, const Functor3& vPhi,
        container& vx, container& vy, container& vz,
        const Geometry& g)
{
    using host_vec = typename Geometry::host_vector;
    host_vec out1 = pullback( vR, g);
    host_vec out2 = pullback( vZ, g);
    host_vec out3 = pullback( vPhi, g);
    dg::tensor::multiply3d(g.jacobian(), out1, out2, out3, out1, out2, out3);
    dg::assign( out1, vx);
    dg::assign( out2, vy);
    dg::assign( out3, vz);
}
 
template<class FunctorRR, class FunctorRZ, class FunctorZZ, class container, class Geometry>
void pushForwardPerp( const FunctorRR& chiRR, const FunctorRZ& chiRZ, const FunctorZZ& chiZZ,
        SparseTensor<container>& chi,
        const Geometry& g)
{
    using host_vec = typename Geometry::host_vector;
    host_vec chiRR_ = pullback( chiRR, g);
    host_vec chiRZ_ = pullback( chiRZ, g);
    host_vec chiZZ_ = pullback( chiZZ, g);
 
    const dg::SparseTensor<container> jac = g.jacobian();
    std::vector<container> values( 5);
    dg::assign( dg::evaluate( dg::zero,g), values[0]);
    dg::assign( dg::evaluate( dg::one, g), values[1]);
    dg::assign( chiRR_, values[2]);
    dg::assign( chiRZ_, values[3]);
    dg::assign( chiZZ_, values[4]);
    chi.idx(0,0)=2, chi.idx(0,1)=chi.idx(1,0)=3, chi.idx(1,1)=4;
    chi.idx(2,0)=chi.idx(2,1)=chi.idx(0,2)=chi.idx(1,2) = 0;
    chi.idx(2,2)=1;
    chi.values() = values;
    //we do not need 3rd dimension here
    container tmp00(jac.value(0,0)), tmp01(tmp00), tmp10(tmp00), tmp11(tmp00);
    // multiply Chi*J^T -> tmp ( Matrix-Matrix multiplication: "line x column")
    dg::tensor::multiply2d( chi, jac.value(0,0), jac.value(0,1), tmp00, tmp10);
    dg::tensor::multiply2d( chi, jac.value(1,0), jac.value(1,1), tmp01, tmp11);
    // multiply J * tmp -> Chi
    dg::tensor::multiply2d( jac, tmp00, tmp10, chi.values()[2], chi.values()[3]);
    dg::tensor::multiply2d( jac, tmp01, tmp11, chi.values()[3], chi.values()[4]);
}
 
namespace create{
 
 
//Note that for the volume function to work properly all 2d grids must set the g_22 element to 1!!
 
template< class Geometry>
typename Geometry::host_vector volume( const Geometry& g)
{
    using host_vector = typename Geometry::host_vector;
    host_vector vol = dg::tensor::volume(g.metric());
    host_vector weights = dg::create::weights( g);
    dg::blas1::pointwiseDot( weights, vol, vol);
    return vol;
}
 
template< class Geometry>
typename Geometry::host_vector inv_volume( const Geometry& g)
{
    using host_vector = typename Geometry::host_vector;
    using real_type = get_value_type<host_vector>;
    host_vector vol = volume(g);
    dg::blas1::transform( vol, vol, dg::INVERT<real_type>());
    return vol;
}
 
}//namespace create
 
} //namespace dg