projection.h

Go to the documentation of this file.

#pragma once
#include <vector>
#include "grid.h"
#include "interpolation.h"
#include "weights.h"
#include "fem.h"
 
namespace dg{
 
template<class T>
T gcd( T a, T b)
{
    T r2 = std::max(a,b);
    T r1 = std::min(a,b);
    while( r1!=0)
    {
        r2 = r2%r1;
        std::swap( r1, r2);
    }
    return r2;
}
 
template<class T>
T lcm( T a, T b)
{
    T g = gcd( a,b);
    return a/g*b;
}
 
namespace create{
 
template<class real_type>
dg::SparseMatrix< int, real_type, thrust::host_vector> diagonal( const thrust::host_vector<real_type>& diagonal)
{
    unsigned size = diagonal.size();
    thrust::host_vector<int> A_row_offsets(size+1), A_column_indices( size);
    thrust::host_vector<real_type> A_values( size);
 
    A_row_offsets[0] = 0;
    for( unsigned i=0; i<size; i++)
    {
        A_row_offsets[i+1] = i+1;
        A_column_indices[i] = i;
        A_values[i] = diagonal[i];
    }
    return { size, size, A_row_offsets, A_column_indices, A_values};
}
 
 
template<class real_type, size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector> projection(
    const aRealTopology<real_type,Nd>& g_new,
    const aRealTopology<real_type,Nd>& g_old, std::string method = "dg")
{
    // These tests cannot pass for MPI
    //for( unsigned u=0; u<Nd; u++)
    //{
    //    if( fabs(g_new.h(u) / g_old.h(u) -round(g_new.h(u) / g_old.h(u))) > 1e-13)
    //        throw dg::Error( dg::Message(_ping_)
    //                  << "# ATTENTION: you project between incompatible grids!! old h: "
    //                  <<g_old.h(u)<<" new h: "<<g_new.h(u)<<"\n");
    //    // Check that boundaries of old grid conincide with cell boundaries of new
    //    real_type fp = (g_old.p(u) - g_new.p(u) ) / g_new.h(u); // must be integer
    //    if( fabs(fp - round( fp)) > 1e-13 and method == "dg")
    //        throw dg::Error( dg::Message(_ping_)
    //                  << "# ATTENTION: you project between incompatible grids!! old p: "
    //                  <<g_old.p(u)<<" new p: "<<g_new.p(u)<<" fp "<<fp<<"\n");
    //    real_type fq = (g_old.q(u) - g_new.p(u) ) / g_new.h(u); // must be integer
    //    if( fabs(fq - round( fq)) > 1e-13 and method == "dg")
    //        throw dg::Error( dg::Message(_ping_)
    //                  << "# ATTENTION: you project between incompatible grids!! old q: "
    //                  <<g_old.q(u)<<" new p: "<<g_new.p(u)<<" fq "<<fq<<"\n");
    //    if( g_old.n(u) < g_new.n(u))
    //        throw dg::Error( dg::Message(_ping_)
    //                   << "# ATTENTION: you project between incompatible grids!! old n: "
    //                   <<g_old.n(u)<<" new n: "<<g_new.n(u)<<"\n");
    //}
    //form the adjoint
    auto w_old = dg::create::weights( g_old);
    auto v_new = dg::create::inv_weights( g_new);
    auto project = interpolation( g_old, g_new, method).transpose();
    for( unsigned row=0; row<project.num_rows(); row++)
        for ( int jj = project.row_offsets()[row]; jj < project.row_offsets()[row+1]; jj++)
        {
            int col = project.column_indices()[jj];
            // Note that w_old is multiplied before v_new (keeps results backwards reproducible)
            project.values()[jj] = v_new[row] * ( project.values()[jj]* w_old[col]);
        }
    return project;
}
 
template<class real_type, size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector> transformation(
    const aRealTopology<real_type,Nd>& g_new,
    const aRealTopology<real_type,Nd>& g_old)
{
    std::array<unsigned, Nd> n_lcm, N_lcm;
    for( unsigned u=0; u<Nd; u++)
    {
        n_lcm [u] = lcm( g_new.n(u), g_old.n(u));
        N_lcm [u] = lcm( g_new.N(u), g_old.N(u));
    }
    RealGrid<real_type, Nd> g_lcm ( g_new.get_p(), g_new.get_q(), n_lcm, N_lcm, g_new.get_bc());
    return create::projection( g_new, g_lcm)*create::interpolation( g_lcm, g_old);
}
 
 
template<class real_type, size_t Nd>
dg::IHMatrix_t<real_type> backproject( const aRealTopology<real_type,Nd>& g)
{
    std::array<dg::IHMatrix_t<real_type>,Nd> matrix;
    for( unsigned u=0; u<Nd; u++)
    {
        unsigned n=g.n(u);
        dg::RealGrid1d<real_type> g_old( -1., 1., n, 1);
        dg::RealGrid1d<real_type> g_new( -1., 1., 1, n);
        auto block = dg::create::transformation( g_new, g_old);
        dg::SquareMatrix<real_type> op(n, 0.);
        for( unsigned i=0; i<block.num_rows(); i++)
            for( unsigned j=block.row_offsets()[i]; j<(unsigned)block.row_offsets()[i+1]; j++)
                op( i, block.column_indices()[j]) = block.values()[j];
        matrix[u] = (dg::IHMatrix_t<real_type>)dg::tensorproduct( g.N(u), op);
 
    }
    for( unsigned u=1; u<Nd; u++)
        matrix[0] = dg::tensorproduct( matrix[u], matrix[0]);
    return matrix[0];
}
 
template<class real_type, size_t Nd>
dg::IHMatrix_t<real_type> inv_backproject( const aRealTopology<real_type,Nd>& g)
{
    std::array<dg::IHMatrix_t<real_type>,Nd> matrix;
    for( unsigned u=0; u<Nd; u++)
    {
        unsigned n=g.n(u);
        dg::RealGrid1d<real_type> g_old( -1., 1., n, 1);
        dg::RealGrid1d<real_type> g_new( -1., 1., 1, n);
        auto block = dg::create::transformation( g_new, g_old);
        dg::SquareMatrix<real_type> op(n, 0.);
        for( unsigned i=0; i<block.num_rows(); i++)
            for( unsigned j=block.row_offsets()[i]; j<(unsigned)block.row_offsets()[i+1]; j++)
                op( i, block.column_indices()[j]) = block.values()[j];
        matrix[u] = (dg::IHMatrix_t<real_type>)dg::tensorproduct( g.N(u), dg::invert(op));
 
    }
    for( unsigned u=1; u<Nd; u++)
        matrix[0] = dg::tensorproduct( matrix[u], matrix[0]);
    return matrix[0];
}
 
 
}//namespace create
}//namespace dg