matrix_out.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 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_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;

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

  virtual ~MatrixOut() override = default;

  template <class Matrix>
  void
  build_patches(const Matrix &     matrix,
                const std::string &name,
                const Options      options = Options(false, 1, false));

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 gridpoints_x = (matrix.n() / options.block_size +
                            (matrix.n() % options.block_size != 0 ? 1 : 0)),
            gridpoints_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)
    {
      --gridpoints_x;
      --gridpoints_y;
    }

  // first clear old data and set it
  // to virgin state
  patches.clear();
  patches.resize((gridpoints_x) * (gridpoints_y));

  // now build the patches
  size_type index = 0;
  for (size_type i = 0; i < gridpoints_y; ++i)
    for (size_type j = 0; j < gridpoints_x; ++j, ++index)
      {
        // 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 rim, 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 (unsigned int v = 0; v < 4; ++v)
          patches[index].vertices[v] *= 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);
          }
      };

  // finally set the name
  this->name = name;
}



/*----------------------------   matrix_out.h     ---------------------------*/

DEAL_II_NAMESPACE_CLOSE

#endif
/*----------------------------   matrix_out.h     ---------------------------*/