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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/sparse_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>

#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_block_sparse_matrix.h>
#endif


DEAL_II_NAMESPACE_OPEN

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

  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 () = default;

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

private:

  typedef DataOutBase::Patch<2,2> Patch;

  std::vector<Patch> patches;

  std::string name;

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

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

  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
  {
    namespace
    {
      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);
      }
    }
  }
}



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     ---------------------------*/