local_results.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#ifndef dealii__mesh_worker_local_results_h
#define dealii__mesh_worker_local_results_h

#include <deal.II/base/config.h>
#include <deal.II/base/std_cxx11/function.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/lac/matrix_block.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/meshworker/vector_selector.h>

DEAL_II_NAMESPACE_OPEN

class BlockIndices;
template<int,int> class DoFHandler;


namespace MeshWorker
{
  template <typename number>
  class LocalResults
  {
  public:
    unsigned int n_values () const;

    unsigned int n_vectors () const;

    unsigned int n_matrices () const;

    unsigned int n_quadrature_points() const;

    unsigned int n_quadrature_values() const;

    number &value(unsigned int i);

    number value(unsigned int i) const;

    BlockVector<number> &vector(unsigned int i);

    const BlockVector<number> &vector(unsigned int i) const;

    MatrixBlock<FullMatrix<number> > &matrix(unsigned int i, bool external = false);

    const MatrixBlock<FullMatrix<number> > &matrix(unsigned int i, bool external = false) const;

    Table<2, number> &quadrature_values();

    number &quadrature_value(unsigned int k, unsigned int i);

    number quadrature_value(unsigned int k, unsigned int i) const;

    void initialize_numbers(const unsigned int n);

    void initialize_vectors(const unsigned int n);

    void initialize_matrices(unsigned int n, bool both);

    template <typename MatrixType>
    void initialize_matrices(const MatrixBlockVector<MatrixType> &matrices,
                             bool both);

    template <typename MatrixType>
    void initialize_matrices(const MGMatrixBlockVector<MatrixType> &matrices,
                             bool both);

    void initialize_quadrature(unsigned int np, unsigned int nv);

    void reinit(const BlockIndices &local_sizes);

    template <class StreamType>
    void print_debug(StreamType &os) const;

    std::size_t memory_consumption () const;

  private:
    void initialize_local(MatrixBlock<FullMatrix<number> > &M,
                          const unsigned int row,
                          const unsigned int col);

    std::vector<number> J;

    std::vector<BlockVector<number> > R;

    std::vector<MatrixBlock<FullMatrix<number> > > M1;

    std::vector<MatrixBlock<FullMatrix<number> > > M2;

    Table<2, number> quadrature_data;
  };

//----------------------------------------------------------------------//

  template <typename number>
  inline void
  LocalResults<number>::initialize_numbers(unsigned int n)
  {
    J.resize(n);
  }


  template <typename number>
  inline void
  LocalResults<number>::initialize_vectors(const unsigned int n)
  {
    R.resize(n);
  }


  template <typename number>
  template <typename MatrixType>
  inline void
  LocalResults<number>::initialize_matrices
  (const MatrixBlockVector<MatrixType> &matrices,
   bool                                 both)
  {
    M1.resize(matrices.size());
    if (both)
      M2.resize(matrices.size());
    for (unsigned int i=0; i<matrices.size(); ++i)
      {
        const unsigned int row = matrices.block(i).row;
        const unsigned int col = matrices.block(i).column;

        M1[i].row = row;
        M1[i].column = col;
        if (both)
          {
            M2[i].row = row;
            M2[i].column = col;
          }
      }
  }


  template <typename number>
  template <typename MatrixType>
  inline void
  LocalResults<number>::initialize_matrices
  (const MGMatrixBlockVector<MatrixType> &matrices,
   bool                                   both)
  {
    M1.resize(matrices.size());
    if (both)
      M2.resize(matrices.size());
    for (unsigned int i=0; i<matrices.size(); ++i)
      {
        const MGLevelObject<MatrixBlock<MatrixType> > &o = matrices.block(i);
        const unsigned int row = o[o.min_level()].row;
        const unsigned int col = o[o.min_level()].column;

        M1[i].row = row;
        M1[i].column = col;
        if (both)
          {
            M2[i].row = row;
            M2[i].column = col;
          }
      }
  }


  template <typename number>
  inline void
  LocalResults<number>::initialize_matrices(const unsigned int n,
                                            const bool both)
  {
    M1.resize(n);
    if (both)
      M2.resize(n);
    for (unsigned int i=0; i<n; ++i)
      {
        M1[i].row = 0;
        M1[i].column = 0;
        if (both)
          {
            M2[i].row = 0;
            M2[i].column = 0;
          }
      }
  }


  template <typename number>
  inline void
  LocalResults<number>::initialize_quadrature(unsigned int np, unsigned int nv)
  {
    quadrature_data.reinit(np, nv);
  }


  template <typename number>
  inline
  unsigned int
  LocalResults<number>::n_values() const
  {
    return J.size();
  }


  template <typename number>
  inline
  unsigned int
  LocalResults<number>::n_vectors() const
  {
    return R.size();
  }


  template <typename number>
  inline
  unsigned int
  LocalResults<number>::n_matrices() const
  {
    return M1.size();
  }


  template <typename number>
  inline
  unsigned int
  LocalResults<number>::n_quadrature_points() const
  {
    return quadrature_data.n_rows();
  }


  template <typename number>
  inline
  unsigned int
  LocalResults<number>::n_quadrature_values() const
  {
    return quadrature_data.n_cols();
  }


  template <typename number>
  inline
  number &
  LocalResults<number>::value(unsigned int i)
  {
    AssertIndexRange(i,J.size());
    return J[i];
  }


  template <typename number>
  inline
  BlockVector<number> &
  LocalResults<number>::vector(unsigned int i)
  {
    AssertIndexRange(i,R.size());
    return R[i];
  }


  template <typename number>
  inline
  MatrixBlock<FullMatrix<number> > &
  LocalResults<number>::matrix(unsigned int i, bool external)
  {
    if (external)
      {
        AssertIndexRange(i,M2.size());
        return M2[i];
      }
    AssertIndexRange(i,M1.size());
    return M1[i];
  }


  template <typename number>
  inline
  number &
  LocalResults<number>::quadrature_value(unsigned int k, unsigned int i)
  {
    return quadrature_data(k,i);
  }


  template <typename number>
  inline
  Table<2, number> &
  LocalResults<number>::quadrature_values()
  {
    return quadrature_data;
  }


  template <typename number>
  inline
  number
  LocalResults<number>::value(unsigned int i) const
  {
    AssertIndexRange(i,J.size());
    return J[i];
  }


  template <typename number>
  inline
  const BlockVector<number> &
  LocalResults<number>::vector(unsigned int i) const
  {
    AssertIndexRange(i,R.size());
    return R[i];
  }


  template <typename number>
  inline
  const MatrixBlock<FullMatrix<number> > &
  LocalResults<number>::matrix(unsigned int i, bool external) const
  {
    if (external)
      {
        AssertIndexRange(i,M2.size());
        return M2[i];
      }
    AssertIndexRange(i,M1.size());
    return M1[i];
  }


  template <typename number>
  inline
  number
  LocalResults<number>::quadrature_value(unsigned int k, unsigned int i) const
  {
    return quadrature_data(k,i);
  }


  template <typename number>
  template <class StreamType>
  void
  LocalResults<number>::print_debug(StreamType &os) const
  {
    os << "J: " << J.size() << std::endl;
    os << "R: " << R.size() << std::endl;
    for (unsigned int i=0; i<R.size(); ++i)
      {
        os << "  " << R[i].n_blocks() << " -";
        for (unsigned int j=0; j<R[i].n_blocks(); ++j)
          os << ' ' << R[i].block(j).size();
        os << std::endl;
      }
    os << "M: " << M1.size() << " face " << M2.size() << std::endl;
    for (unsigned int i=0; i<M1.size(); ++i)
      {
        os << "  " << M1[i].row << "," << M1[i].column
           << " " << M1[i].matrix.m() << 'x' << M1[i].matrix.n();
        if (i < M2.size())
          os << " face " << M2[i].row << "," << M2[i].column
             << " " << M2[i].matrix.m() << 'x' << M2[i].matrix.n();
        os << std::endl;
      }
  }

}


DEAL_II_NAMESPACE_CLOSE

#endif