local_results.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2018 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_local_results_h
#define dealii_mesh_worker_local_results_h

#include <deal.II/base/config.h>

#include <deal.II/base/geometry_info.h>

#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/matrix_block.h>

#include <deal.II/meshworker/vector_selector.h>

#include <functional>

DEAL_II_NAMESPACE_OPEN

class BlockIndices;

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:
    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;
      }
  }

} // namespace MeshWorker


DEAL_II_NAMESPACE_CLOSE

#endif