reshape.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2024-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef RESHAPE_HH
#define RESHAPE_HH
 
#include "dune/istl/bvector.hh"
 
#include <cassert>
 
namespace Kaskade
{
 
  template <class Entry>
  class DynamicMatrix;
 
 
  template <int n, class Entry, int m, class A>
  Dune::BlockVector<Dune::FieldVector<Entry,n>,A>
  reshapeBlocks(Dune::BlockVector<Dune::FieldVector<Entry,m>,A> const& x)
  {
    assert( (x.N()*m) % n == 0);
    Dune::BlockVector<Dune::FieldVector<Entry,n>,A> y(x.N()*m/n);
    for (int i=0; i<x.N(); ++i)
      for (int j=0; j<m; ++j)
      {
        int k = (m*i+j) / n;
        int l = (m*i+j) % n;
        y[k][l] = x[i][j];
      }
 
    return y;
  }
 
  // ----------------------------------------------------------------------------------------------
 
  template <class Kb, int nb, int mb, class Ka, int na, int ma>
  DynamicMatrix<Dune::FieldMatrix<Kb,nb,mb>>
  reshapeBlocks(DynamicMatrix<Dune::FieldMatrix<Ka,na,ma>> const& A)
  {
    assert( (A.N()*na) % nb == 0 );
    assert( (A.M()*ma) % mb == 0 );
 
    DynamicMatrix<Dune::FieldMatrix<Kb,nb,mb>> B(A.N()*na/nb,A.M()*ma/mb);
    for (int j=0; j<A.M(); ++j)
      for (int i=0; i<A.N(); ++i)
        for (int k=0; k<na; ++k)
          for (int l=0; l<ma; ++l)
          {
            int r = i*na+k;
            int c = j*ma+l;
            B[r/nb][c/mb][r%nb][c%mb] = A[i][j][k][l];
          }
    return B;
  }
 
}
 
#endif