base_geometry.h

Go to the documentation of this file.

#pragma once
 
#include "grid.h"
#include "tensor.h"
 
namespace dg
{
 
 
template<class real_type>
struct aRealGeometry2d : public aRealTopology<real_type,2>
{
    SparseTensor<thrust::host_vector<real_type> > jacobian()const{
        return do_compute_jacobian();
    }
    SparseTensor<thrust::host_vector<real_type> > metric()const {
        return do_compute_metric();
    }
    std::vector<thrust::host_vector<real_type> > map()const{
        return do_compute_map();
    }
    virtual aRealGeometry2d* clone()const=0;
    virtual ~aRealGeometry2d() = default;
    protected:
    using aRealTopology<real_type,2>::aRealTopology;
    aRealGeometry2d( const aRealGeometry2d& src) = default;
    aRealGeometry2d& operator=( const aRealGeometry2d& src) = default;
    private:
    virtual SparseTensor<thrust::host_vector<real_type> > do_compute_metric()const {
        return SparseTensor<thrust::host_vector<real_type> >(*this);
    }
    virtual SparseTensor<thrust::host_vector<real_type> > do_compute_jacobian()const {
        return SparseTensor<thrust::host_vector<real_type> >(*this);
    }
    virtual std::vector<thrust::host_vector<real_type> > do_compute_map()const{
        std::vector<thrust::host_vector<real_type> > map(2);
        map[0] = dg::evaluate(dg::cooX2d, *this);
        map[1] = dg::evaluate(dg::cooY2d, *this);
        return map;
    }
 
 
};
 
template<class real_type>
struct aRealGeometry3d : public aRealTopology<real_type,3>
{
    SparseTensor<thrust::host_vector<real_type> > jacobian()const{
        return do_compute_jacobian();
    }
    SparseTensor<thrust::host_vector<real_type> > metric()const {
        return do_compute_metric();
    }
    std::vector<thrust::host_vector<real_type> > map()const{
        return do_compute_map();
    }
    virtual aRealGeometry3d* clone()const=0;
    virtual ~aRealGeometry3d() = default;
    protected:
    using aRealTopology<real_type,3>::aRealTopology;
    aRealGeometry3d( const aRealGeometry3d& src) = default;
    aRealGeometry3d& operator=( const aRealGeometry3d& src) = default;
    private:
    virtual SparseTensor<thrust::host_vector<real_type> > do_compute_metric()const {
        return SparseTensor<thrust::host_vector<real_type> >(*this);
    }
    virtual SparseTensor<thrust::host_vector<real_type> > do_compute_jacobian()const {
        return SparseTensor<thrust::host_vector<real_type> >(*this);
    }
    virtual std::vector<thrust::host_vector<real_type> > do_compute_map()const{
        std::vector<thrust::host_vector<real_type> > map(3);
        map[0] = dg::evaluate(dg::cooX3d, *this);
        map[1] = dg::evaluate(dg::cooY3d, *this);
        map[2] = dg::evaluate(dg::cooZ3d, *this);
        return map;
    }
};
 
template<class real_type>
struct aRealProductGeometry3d : public aRealGeometry3d<real_type>
{
    aRealGeometry2d<real_type>* perp_grid()const{
        return do_perp_grid();
    }
    virtual ~aRealProductGeometry3d() = default;
    virtual aRealProductGeometry3d* clone()const=0;
    protected:
    using aRealGeometry3d<real_type>::aRealGeometry3d;
    aRealProductGeometry3d( const aRealProductGeometry3d& src) = default;
    aRealProductGeometry3d& operator=( const aRealProductGeometry3d& src) = default;
    private:
    virtual aRealGeometry2d<real_type>* do_perp_grid()const=0;
};
 
template<class real_type>
struct RealCartesianGrid2d: public dg::aRealGeometry2d<real_type>
{
    RealCartesianGrid2d() = default;
    RealCartesianGrid2d( real_type x0, real_type x1, real_type y0, real_type y1, unsigned n, unsigned Nx, unsigned Ny, bc bcx = PER, bc bcy = PER): dg::aRealGeometry2d<real_type>({x0,y0},{x1,y1},{n,n},{Nx,Ny},{bcx,bcy}){}
 
    RealCartesianGrid2d( RealGrid1d<real_type> gx, RealGrid1d<real_type> gy):
        dg::aRealGeometry2d<real_type>(std::array{gx,gy}){}
    RealCartesianGrid2d( const dg::RealGrid2d<real_type>& g):
        dg::aRealGeometry2d<real_type>(std::array{g.gx(), g.gy()}){}
    virtual RealCartesianGrid2d<real_type>* clone()const override final{
        return new RealCartesianGrid2d<real_type>(*this);
    }
    private:
    virtual void do_set(std::array<unsigned,2> new_n, std::array<unsigned,2> new_N) override final{
        aRealTopology<real_type,2>::do_set(new_n,new_N);
    }
    virtual void do_set_pq( std::array<real_type,2> new_x0, std::array<real_type,2> new_x1) override final{
        aRealTopology<real_type,2>::do_set_pq(new_x0,new_x1);
    }
    virtual void do_set( std::array<dg::bc,2> new_bcs) override final{
        aRealTopology<real_type,2>::do_set(new_bcs);
    }
};
 
template<class real_type>
struct RealCartesianGrid3d: public dg::aRealProductGeometry3d<real_type>
{
    using perpendicular_grid = RealCartesianGrid2d<real_type>;
    RealCartesianGrid3d() = default;
    RealCartesianGrid3d( real_type x0, real_type x1, real_type y0, real_type y1, real_type z0, real_type z1,
    unsigned n, unsigned Nx, unsigned Ny, unsigned Nz,
    bc bcx = PER, bc bcy = PER, bc bcz = PER):
        dg::aRealProductGeometry3d<real_type>({x0,y0,z0},{x1,y1,z1},{n,n,1},{Nx,Ny,Nz},{bcx,bcy,bcz})
        {}
 
    RealCartesianGrid3d( RealGrid1d<real_type> gx, RealGrid1d<real_type> gy, RealGrid1d<real_type> gz):
        dg::aRealProductGeometry3d<real_type>(std::array{gx,gy,gz}){}
    RealCartesianGrid3d( const dg::RealGrid3d<real_type>& g):
        dg::aRealProductGeometry3d<real_type>(std::array{g.gx(), g.gy(), g.gz()}){}
    virtual RealCartesianGrid3d* clone()const override final{
        return new RealCartesianGrid3d(*this);
    }
    private:
    virtual RealCartesianGrid2d<real_type>* do_perp_grid() const override final{
        return new RealCartesianGrid2d<real_type>(this->gx(), this->gy());
    }
    virtual void do_set(std::array<unsigned,3> new_n, std::array<unsigned,3> new_N) override final{
        aRealTopology<real_type,3>::do_set(new_n,new_N);
    }
    virtual void do_set_pq( std::array<real_type,3> new_x0, std::array<real_type,3> new_x1) override final{
        aRealTopology<real_type,3>::do_set_pq(new_x0,new_x1);
    }
    virtual void do_set( std::array<dg::bc,3> new_bcs) override final{
        aRealTopology<real_type,3>::do_set(new_bcs);
    }
};
 
template<class real_type>
struct RealCylindricalGrid3d: public dg::aRealProductGeometry3d<real_type>
{
    using perpendicular_grid = RealCartesianGrid2d<real_type>;
    RealCylindricalGrid3d() = default;
    RealCylindricalGrid3d( real_type x0, real_type x1,
        real_type y0, real_type y1, real_type z0, real_type z1,
        unsigned n, unsigned Nx, unsigned Ny, unsigned Nz,
        bc bcx = PER, bc bcy = PER, bc bcz = PER):
    dg::aRealProductGeometry3d<real_type>({x0,y0,z0},{x1,y1,z1},{n,n,1},{Nx,Ny,Nz},{bcx,bcy,bcz})
    {}
    RealCylindricalGrid3d( RealGrid1d<real_type> gx, RealGrid1d<real_type> gy, RealGrid1d<real_type> gz):
        dg::aRealProductGeometry3d<real_type>(std::array{gx,gy,gz}){}
    virtual RealCylindricalGrid3d* clone()const override final{
        return new RealCylindricalGrid3d(*this);
    }
    private:
    virtual RealCartesianGrid2d<real_type>* do_perp_grid() const override final{
        return new RealCartesianGrid2d<real_type>(this->gx(), this->gy());
    }
    virtual SparseTensor<thrust::host_vector<real_type> > do_compute_metric()const override final{
        SparseTensor<thrust::host_vector<real_type> > metric(*this);
        thrust::host_vector<real_type> R = dg::evaluate(dg::cooX3d, *this);
        for( unsigned i = 0; i<this->size(); i++)
            R[i] = 1./R[i]/R[i];
        metric.idx(2,2)=2;
        metric.values().push_back( R);
        return metric;
    }
    virtual void do_set(std::array<unsigned,3> new_n, std::array<unsigned,3> new_N) override final{
        aRealTopology3d<real_type>::do_set(new_n,new_N);
    }
    virtual void do_set_pq( std::array<real_type,3> new_x0, std::array<real_type,3> new_x1) override final{
        aRealTopology<real_type,3>::do_set_pq(new_x0,new_x1);
    }
    virtual void do_set( std::array<dg::bc,3> new_bcs) override final{
        aRealTopology<real_type,3>::do_set(new_bcs);
    }
};
 
 
//TODO With C++17 class type inference we could rethink these?
using aGeometry2d           = dg::aRealGeometry2d<double>;
using aGeometry3d           = dg::aRealGeometry3d<double>;
using aProductGeometry3d    = dg::aRealProductGeometry3d<double>;
using CartesianGrid2d       = dg::RealCartesianGrid2d<double>;
using CartesianGrid3d       = dg::RealCartesianGrid3d<double>;
using CylindricalGrid3d     = dg::RealCylindricalGrid3d<double>;
#ifndef MPI_VERSION
namespace x{
using aGeometry2d           = aGeometry2d           ;
using aGeometry3d           = aGeometry3d           ;
using aProductGeometry3d    = aProductGeometry3d    ;
using CartesianGrid2d       = CartesianGrid2d       ;
using CartesianGrid3d       = CartesianGrid3d       ;
using CylindricalGrid3d     = CylindricalGrid3d     ;
}//namespace x
#endif //MPI_VERSION
 
 
} //namespace dg