prolongation.h

Go to the documentation of this file.

#pragma once
 
#include "grid.h"
#include "weights.h"
#include "dg/blas1.h"
#include "interpolation.h"
#include "projection.h"
namespace dg{
namespace create{
namespace detail{
template< class Topology, size_t M, size_t ...I>
auto do_full_abscissas( const Topology& g, std::array<unsigned, M> map,
        std::index_sequence<I...>)
{
    std::array< decltype(g.abscissas(0)),Topology::ndim()> abs;
    for( unsigned u=0; u<Topology::ndim(); u++)
        abs[u] = g.abscissas(u);
    auto abs_red = abs;
    for( unsigned u=0; u<Topology::ndim(); u++)
        dg::blas1::copy( 1, abs_red[u]);
    std::array<decltype(g.abscissas(0)),M> full_abs;
    for( unsigned u=0; u<M; u++)
    {
        dg::blas1::copy( abs[map[u]], abs_red[map[u]]);
        full_abs[u] = dg::kronecker( dg::Product(), abs_red[I]...);
        dg::blas1::copy( 1, abs_red[map[u]]);
    }
    return full_abs;
}
 
template<class Topology, size_t Md>
std::array<bool,Topology::ndim()> axes2bools( const Topology& g, std::array<unsigned,Md> axes)
{
    // The axes that are averaged
    std::array<bool,Topology::ndim()> bools;
    for( unsigned u=0; u<Topology::ndim(); u++)
        bools[u] = false;
    for( unsigned u=0; u<Md; u++)
    {
        assert( axes[u] < bools.size());
        bools[axes[u]] = true;
    }
    return bools;
}
template<class Topology, size_t Md>
std::array<unsigned,Topology::ndim()-Md> complement( const Topology& g, std::array<unsigned,Md> axes)
{
    // The axes that remain
    std::array<unsigned,Topology::ndim()-Md> comp = {0u};
    auto bools = axes2bools( g, axes);
    unsigned counter = 0;
    for ( unsigned u=0; u<Topology::ndim(); u++)
    {
        if( !bools[u])
        {
            comp[counter] = u;
            counter++;
        }
    }
    return comp;
}
 
} // namespace detail
 
template<class real_type, size_t Nd, size_t Md>
dg::SparseMatrix<int, real_type, thrust::host_vector> prolongation(
    const aRealTopology<real_type,Nd>& g_new,
    std::array<unsigned,Md> axes)
{
    static_assert( Md < Nd && Md > 0 && Nd > 0);
    std::array<unsigned,Nd-Md> remains = detail::complement( g_new, axes);
    auto full_abs = detail::do_full_abscissas( g_new, remains, std::make_index_sequence<Nd>());
    std::vector<dg::SparseMatrix<int,real_type,thrust::host_vector>> matrix(Nd-Md);
    for( unsigned u=0; u<Nd-Md; u++)
    {
        matrix[u] = detail::interpolation1d( dg::xspace, full_abs[u],
            g_new.grid(remains[u]), g_new.bc(remains[u]));
    }
    for( unsigned u=1; u<Nd-Md; u++)
        matrix[0] = dg::tensorproduct_cols( matrix[u], matrix[0]);
    return matrix[0];
}
 
template<class real_type, size_t Nd, size_t Md>
dg::SparseMatrix<int, real_type, thrust::host_vector> reduction(
    std::array<unsigned,Md> axes,
    const aRealTopology<real_type,Nd>& g_old)
{
    return prolongation( g_old, axes).transpose();
}
 
template<class real_type, size_t Nd, size_t Md>
dg::SparseMatrix<int, real_type, thrust::host_vector> projection(
    std::array<unsigned,Md> axes,
    const aRealTopology<real_type,Nd>& g_old)
{
    std::array<unsigned,Nd-Md> remains = detail::complement( g_old, axes);
    std::array<RealGrid<real_type,1>,Nd-Md> gs;
    for( unsigned u=0; u<Nd-Md; u++)
        gs[u] = g_old.grid( remains[u]);
    RealGrid<real_type, Nd-Md> g_new(gs);
    auto w_old = dg::create::weights( g_old);
    auto v_new = dg::create::inv_weights( g_new);
    auto project = prolongation( g_old, axes).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;
}
 
} // namespace create
 
}//namespace dg