local_results.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2023 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
 
// Forward declaration
#ifndef DOXYGEN
class BlockIndices;
#endif
 
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(const unsigned int i);
 
    number
    value(const unsigned int i) const;
 
    BlockVector<number> &
    vector(const unsigned int i);
 
    const BlockVector<number> &
    vector(const unsigned int i) const;
 
    MatrixBlock<FullMatrix<number>> &
    matrix(const unsigned int i, const bool external = false);
 
    const MatrixBlock<FullMatrix<number>> &
    matrix(const unsigned int i, const bool external = false) const;
 
    Table<2, number> &
    quadrature_values();
 
    number &
    quadrature_value(const unsigned int k, const unsigned int i);
 
    number
    quadrature_value(const unsigned int k, const unsigned int i) const;
 
    void
    initialize_numbers(const unsigned int n);
 
    void
    initialize_vectors(const unsigned int n);
 
    void
    initialize_matrices(const 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(const unsigned int np, const 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(const 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(const unsigned int np,
                                              const 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(const unsigned int i)
  {
    AssertIndexRange(i, J.size());
    return J[i];
  }
 
 
  template <typename number>
  inline BlockVector<number> &
  LocalResults<number>::vector(const unsigned int i)
  {
    AssertIndexRange(i, R.size());
    return R[i];
  }
 
 
  template <typename number>
  inline MatrixBlock<FullMatrix<number>> &
  LocalResults<number>::matrix(const unsigned int i, const 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(const unsigned int k,
                                         const 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(const unsigned int i) const
  {
    AssertIndexRange(i, J.size());
    return J[i];
  }
 
 
  template <typename number>
  inline const BlockVector<number> &
  LocalResults<number>::vector(const unsigned int i) const
  {
    AssertIndexRange(i, R.size());
    return R[i];
  }
 
 
  template <typename number>
  inline const MatrixBlock<FullMatrix<number>> &
  LocalResults<number>::matrix(const unsigned int i, const 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(const unsigned int k,
                                         const 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