matrix_out.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2001 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_matrix_out_h
#  define dealii_matrix_out_h
 
#  include <deal.II/base/config.h>
 
#  include <deal.II/base/data_out_base.h>
 
#  include <deal.II/lac/block_sparse_matrix.h>
#  include <deal.II/lac/sparse_matrix.h>
 
#  ifdef DEAL_II_WITH_TRILINOS
#    include <deal.II/lac/trilinos_block_sparse_matrix.h>
#    include <deal.II/lac/trilinos_sparse_matrix.h>
#  endif
 
 
DEAL_II_NAMESPACE_OPEN

class MatrixOut : public DataOutInterface<2, 2>
{
public:
  using size_type = types::global_dof_index;

  struct Options
  {
    bool show_absolute_values;

    unsigned int block_size;

    bool discontinuous;

    bool create_sparse_plot;

    Options(const bool         show_absolute_values = false,
            const unsigned int block_size           = 1,
            const bool         discontinuous        = false,
            const bool         create_sparse_plot   = true);
  };

  virtual ~MatrixOut() override = default;

  template <class Matrix>
  void
  build_patches(const Matrix      &matrix,
                const std::string &name,
                const Options      options = Options());
 
private:
  using Patch = DataOutBase::Patch<2, 2>;

  std::vector<Patch> patches;

  std::string name;

  virtual const std::vector<Patch> &
  get_patches() const override;

  virtual std::vector<std::string>
  get_dataset_names() const override;

  template <class Matrix>
  static double
  get_gridpoint_value(const Matrix   &matrix,
                      const size_type i,
                      const size_type j,
                      const Options  &options);
};
 
 
/* ---------------------- Template and inline functions ------------- */
 
 
namespace internal
{
  namespace MatrixOutImplementation
  {
    template <typename number>
    double
    get_element(const ::SparseMatrix<number> &matrix,
                const types::global_dof_index       i,
                const types::global_dof_index       j)
    {
      return matrix.el(i, j);
    }
 
 

    template <typename number>
    double
    get_element(const ::BlockSparseMatrix<number> &matrix,
                const types::global_dof_index            i,
                const types::global_dof_index            j)
    {
      return matrix.el(i, j);
    }
 
 
#  ifdef DEAL_II_WITH_TRILINOS
    inline double
    get_element(const TrilinosWrappers::SparseMatrix &matrix,
                const types::global_dof_index         i,
                const types::global_dof_index         j)
    {
      return matrix.el(i, j);
    }
 
 

    inline double
    get_element(const TrilinosWrappers::BlockSparseMatrix &matrix,
                const types::global_dof_index              i,
                const types::global_dof_index              j)
    {
      return matrix.el(i, j);
    }
#  endif
 
 
#  ifdef DEAL_II_WITH_PETSC
    // no need to do anything: PETSc matrix objects do not distinguish
    // between operator() and el(i,j), so we can safely access elements
    // through the generic function below
#  endif
 

    template <class Matrix>
    double
    get_element(const Matrix                 &matrix,
                const types::global_dof_index i,
                const types::global_dof_index j)
    {
      return matrix(i, j);
    }
  } // namespace MatrixOutImplementation
} // namespace internal
 
 
 
template <class Matrix>
inline double
MatrixOut::get_gridpoint_value(const Matrix   &matrix,
                               const size_type i,
                               const size_type j,
                               const Options  &options)
{
  // special case if block size is
  // one since we then don't need all
  // that loop overhead
  if (options.block_size == 1)
    {
      if (options.show_absolute_values == true)
        return std::fabs(
          internal::MatrixOutImplementation::get_element(matrix, i, j));
      else
        return internal::MatrixOutImplementation::get_element(matrix, i, j);
    }
 
  // if blocksize greater than one,
  // then compute average of elements
  double    average    = 0;
  size_type n_elements = 0;
  for (size_type row = i * options.block_size;
       row <
       std::min(size_type(matrix.m()), size_type((i + 1) * options.block_size));
       ++row)
    for (size_type col = j * options.block_size;
         col < std::min(size_type(matrix.m()),
                        size_type((j + 1) * options.block_size));
         ++col, ++n_elements)
      if (options.show_absolute_values == true)
        average += std::fabs(
          internal::MatrixOutImplementation::get_element(matrix, row, col));
      else
        average +=
          internal::MatrixOutImplementation::get_element(matrix, row, col);
  average /= n_elements;
  return average;
}
 
 
 
template <class Matrix>
void
MatrixOut::build_patches(const Matrix      &matrix,
                         const std::string &name,
                         const Options      options)
{
  size_type n_patches_x = (matrix.n() / options.block_size +
                           (matrix.n() % options.block_size != 0 ? 1 : 0)),
            n_patches_y = (matrix.m() / options.block_size +
                           (matrix.m() % options.block_size != 0 ? 1 : 0));
 
  // If continuous, the number of
  // plotted patches is matrix size-1
  if (!options.discontinuous)
    {
      --n_patches_x;
      --n_patches_y;
    }
 
  const size_type n_patches =
    (options.create_sparse_plot ?
       [&]() {
         size_type count = 0;
         for (size_type i = 0; i < n_patches_y; ++i)
           {
             for (size_type j = 0; j < n_patches_x; ++j)
               // Use the same logic as below to determine whether we
               // need to output a patch, and count if we do:
               if ((((options.discontinuous == true) &&
                     (get_gridpoint_value(matrix, i, j, options) != 0)) ||
                    ((options.discontinuous == false) &&
                     ((get_gridpoint_value(matrix, i, j, options) != 0) ||
                      (get_gridpoint_value(matrix, i + 1, j, options) != 0) ||
                      (get_gridpoint_value(matrix, i, j + 1, options) != 0) ||
                      (get_gridpoint_value(matrix, i + 1, j + 1, options) !=
                       0)))))
                 ++count;
           }
         return count;
       }() :
       n_patches_x * n_patches_y);
 
  // first clear old data and re-set the object to a correctly sized state:
  patches.clear();
  try
    {
      patches.resize(n_patches);
    }
  catch (const std::bad_alloc &)
    {
      AssertThrow(false,
                  ExcMessage("You are trying to create a graphical "
                             "representation of a matrix that would "
                             "requiring outputting " +
                             (options.create_sparse_plot ?
                                std::to_string(n_patches) :
                                std::to_string(n_patches_x) + "x" +
                                  std::to_string(n_patches_y)) +
                             " patches. There is not enough memory to output " +
                             "this many patches."));
    }
 
  // now build the patches
  size_type index = 0;
  for (size_type i = 0; i < n_patches_y; ++i)
    for (size_type j = 0; j < n_patches_x; ++j)
      {
        // If we are creating a sparse plot, check whether this patch
        // would have any nonzero values. If not, we can skip the
        // patch:
        if (options.create_sparse_plot &&
            (((options.discontinuous == true) &&
              (get_gridpoint_value(matrix, i, j, options) == 0)) ||
             ((options.discontinuous == false) &&
              (get_gridpoint_value(matrix, i, j, options) == 0) &&
              (get_gridpoint_value(matrix, i + 1, j, options) == 0) &&
              (get_gridpoint_value(matrix, i, j + 1, options) == 0) &&
              (get_gridpoint_value(matrix, i + 1, j + 1, options) == 0))))
          continue;
 
        patches[index].n_subdivisions = 1;
        patches[index].reference_cell = ReferenceCells::Quadrilateral;
 
        // within each patch, order the points in such a way that if some
        // graphical output program (such as gnuplot) plots the quadrilaterals
        // as two triangles, then the diagonal of the quadrilateral which cuts
        // it into the two printed triangles is parallel to the diagonal of the
        // matrix, rather than perpendicular to it. this has the advantage that,
        // for example, the unit matrix is plotted as a straight ridge, rather
        // than as a series of bumps and valleys along the diagonal
        patches[index].vertices[0][0] = j;
        patches[index].vertices[0][1] = -static_cast<signed int>(i);
        patches[index].vertices[1][0] = j;
        patches[index].vertices[1][1] = -static_cast<signed int>(i + 1);
        patches[index].vertices[2][0] = j + 1;
        patches[index].vertices[2][1] = -static_cast<signed int>(i);
        patches[index].vertices[3][0] = j + 1;
        patches[index].vertices[3][1] = -static_cast<signed int>(i + 1);
        // next scale all the patch
        // coordinates by the block
        // size, to get original
        // coordinates
        for (auto &vertex : patches[index].vertices)
          vertex *= options.block_size;
 
        patches[index].n_subdivisions = 1;
 
        patches[index].data.reinit(1, 4);
        if (options.discontinuous)
          {
            patches[index].data(0, 0) =
              get_gridpoint_value(matrix, i, j, options);
            patches[index].data(0, 1) =
              get_gridpoint_value(matrix, i, j, options);
            patches[index].data(0, 2) =
              get_gridpoint_value(matrix, i, j, options);
            patches[index].data(0, 3) =
              get_gridpoint_value(matrix, i, j, options);
          }
        else
          {
            patches[index].data(0, 0) =
              get_gridpoint_value(matrix, i, j, options);
            patches[index].data(0, 1) =
              get_gridpoint_value(matrix, i + 1, j, options);
            patches[index].data(0, 2) =
              get_gridpoint_value(matrix, i, j + 1, options);
            patches[index].data(0, 3) =
              get_gridpoint_value(matrix, i + 1, j + 1, options);
          }
 
        ++index;
      }
 
  // finally set the name
  this->name = name;
}
 
 
 
/*----------------------------   matrix_out.h     ---------------------------*/
 
DEAL_II_NAMESPACE_CLOSE
 
#endif
/*----------------------------   matrix_out.h     ---------------------------*/