Make MatrixOut also work for Trilinos matrices.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index aee4e94d2635c8a216e88cffab994d000ec7f809..ec481a65f89c9c661b4fe26bf0dd279116cec8e4 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -493,6 +493,11 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: MatrixOut now also works with Trilinos matrices.
+  <br>
+  (Wolfgang Bangerth, 2015/05/11)
+  </li>
+
   <li> Changed: TrilinosWrappers::Vector, TrilinosWrappers::BlockVector, 
   PETScWrappers::Vector, and PETScWrappers::BlockVector are deprecated. Either 
   use the MPI or the deal.II version of the Vector/BlockVector.

diff --git a/include/deal.II/lac/matrix_out.h b/include/deal.II/lac/matrix_out.h
index 00f21e256b8ae027350d0ac209c8768d9d0c84a1..7be4cbbbab5607b10d575acc28472c88a6e9f755 100644 (file)
--- a/include/deal.II/lac/matrix_out.h
+++ b/include/deal.II/lac/matrix_out.h
@@ -21,6 +21,16 @@
 #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
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/lac/petsc_sparse_matrix.h>
+#  include <deal.II/lac/petsc_block_sparse_matrix.h>
+#endif
+
 DEAL_II_NAMESPACE_OPEN
 
 /**
@@ -120,12 +130,12 @@ public:
    * the fields of this structure for more information.
    *
    * Note that this function requires that we can extract elements of the
-   * matrix, which is done using the get_element() function declared below. By
-   * adding specializations, you can extend this class to other matrix classes
-   * which are not presently supported. Furthermore, we need to be able to
-   * extract the size of the matrix, for which we assume that the matrix type
-   * offers member functions <tt>m()</tt> and <tt>n()</tt>, which return the
-   * number of rows and columns, respectively.
+   * matrix, which is done using the get_element() function declared in an
+   * internal namespace. By adding specializations, you can extend this class
+   * to other matrix classes which are not presently supported. Furthermore,
+   * we need to be able to extract the size of the matrix, for which we assume
+   * that the matrix type offers member functions <tt>m()</tt> and
+   * <tt>n()</tt>, which return the number of rows and columns, respectively.
    */
   template <class Matrix>
   void build_patches (const Matrix      &matrix,
@@ -165,32 +175,6 @@ private:
    */
   virtual std::vector<std::string> get_dataset_names () const;
 
-  /**
-   * Return the element with given indices of a sparse matrix.
-   */
-  template <typename number>
-  static double get_element (const SparseMatrix<number> &matrix,
-                             const size_type             i,
-                             const size_type             j);
-
-  /**
-   * Return the element with given indices of a block sparse matrix.
-   */
-  template <typename number>
-  static double get_element (const BlockSparseMatrix<number> &matrix,
-                             const size_type                  i,
-                             const size_type                  j);
-
-  /**
-   * Return the element with given indices from any matrix type for which no
-   * specialization of this function was declared above. This will call
-   * <tt>operator()</tt> on the matrix.
-   */
-  template <class Matrix>
-  static double get_element (const Matrix       &matrix,
-                             const size_type     i,
-                             const size_type     j);
-
   /**
    * Get the value of the matrix at gridpoint <tt>(i,j)</tt>. Depending on the
    * given flags, this can mean different things, for example if only absolute
@@ -209,40 +193,76 @@ private:
 /* ---------------------- Template and inline functions ------------- */
 
 
-template <typename number>
-inline
-double
-MatrixOut::get_element (const SparseMatrix<number> &matrix,
-                        const size_type             i,
-                        const size_type             j)
+namespace internal
 {
-  return matrix.el(i,j);
-}
+  namespace MatrixOut
+  {
+    namespace
+    {
+      /**
+       * Return the element with given indices of a sparse matrix.
+       */
+      template <typename number>
+      double get_element (const dealii::SparseMatrix<number> &matrix,
+                          const types::global_dof_index             i,
+                          const types::global_dof_index             j)
+      {
+        return matrix.el(i,j);
+      }
 
 
 
+      /**
+       * Return the element with given indices of a block sparse matrix.
+       */
+      template <typename number>
+      double get_element (const dealii::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
+      /**
+       * Return the element with given indices of a sparse matrix.
+       */
+      double get_element (const TrilinosWrappers::SparseMatrix &matrix,
+                          const types::global_dof_index             i,
+                          const types::global_dof_index             j)
+      {
+        return matrix.el(i,j);
+      }
 
-template <typename number>
-inline
-double
-MatrixOut::get_element (const BlockSparseMatrix<number> &matrix,
-                        const size_type                  i,
-                        const size_type                  j)
-{
-  return matrix.el(i,j);
-}
 
 
+      /**
+       * Return the element with given indices of a Trilinos block sparse
+       * matrix.
+       */
+      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
 
 
-template <class Matrix>
-inline
-double
-MatrixOut::get_element (const Matrix   &matrix,
-                        const size_type i,
-                        const size_type j)
-{
-  return matrix(i,j);
+      /**
+       * Return the element with given indices from any matrix type for which no
+       * specialization of this function was declared above. This will call
+       * <tt>operator()</tt> on the matrix.
+       */
+      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);
+      }
+    }
+  }
 }
 
 
@@ -261,9 +281,9 @@ MatrixOut::get_gridpoint_value (const Matrix   &matrix,
   if (options.block_size == 1)
     {
       if (options.show_absolute_values == true)
-        return std::fabs(get_element (matrix, i, j));
+        return std::fabs(internal::MatrixOut::get_element (matrix, i, j));
       else
-        return get_element (matrix, i, j);
+        return internal::MatrixOut::get_element (matrix, i, j);
     }
 
   // if blocksize greater than one,
@@ -277,9 +297,9 @@ MatrixOut::get_gridpoint_value (const Matrix   &matrix,
          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(get_element (matrix, row, col));
+        average += std::fabs(internal::MatrixOut::get_element (matrix, row, col));
       else
-        average += get_element (matrix, row, col);
+        average += internal::MatrixOut::get_element (matrix, row, col);
   average /= n_elements;
   return average;
 }

diff --git a/tests/lac/matrix_out_02.cc b/tests/lac/matrix_out_02.cc
new file mode 100644 (file)
index 0000000..c3a8f98
--- /dev/null
+++ b/tests/lac/matrix_out_02.cc
@@ -0,0 +1,62 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+// like matrix_out.cc, but test for Trilinos matrices
+//
+// also test some of the other options of the MatrixOut::Options class
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/matrix_out.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, numbers::invalid_unsigned_int);
+
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  // test for a rectangular sparse
+  // matrix
+  if (true)
+    {
+      TrilinosWrappers::SparsityPattern sparsity (4,8,7);
+      for (unsigned int i=0; i<4; ++i)
+        for (unsigned int j=0; j<8; ++j)
+          if (i==j+1)
+            sparsity.add (i,j);
+      sparsity.compress ();
+
+      TrilinosWrappers::SparseMatrix sparse_matrix(sparsity);
+      for (unsigned int i=0; i<4; ++i)
+        for (unsigned int j=0; j<8; ++j)
+          if (i==j+1)
+           sparse_matrix.set(i,j, i+3*j);
+
+      MatrixOut matrix_out;
+      matrix_out.build_patches (sparse_matrix, "sparse_matrix",
+                                MatrixOut::Options (true, 1, true));
+      matrix_out.write_gnuplot (logfile);
+    }
+}

diff --git a/tests/lac/matrix_out_02.with_trilinos=true.output b/tests/lac/matrix_out_02.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..c0fa8d8
--- /dev/null
+++ b/tests/lac/matrix_out_02.with_trilinos=true.output
@@ -0,0 +1,232 @@
+
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <sparse_matrix> 
+0.00000 0.00000 0.00000 
+0.00000 -1.00000 0.00000 
+
+1.00000 0.00000 0.00000 
+1.00000 -1.00000 0.00000 
+
+
+1.00000 0.00000 0.00000 
+1.00000 -1.00000 0.00000 
+
+2.00000 0.00000 0.00000 
+2.00000 -1.00000 0.00000 
+
+
+2.00000 0.00000 0.00000 
+2.00000 -1.00000 0.00000 
+
+3.00000 0.00000 0.00000 
+3.00000 -1.00000 0.00000 
+
+
+3.00000 0.00000 0.00000 
+3.00000 -1.00000 0.00000 
+
+4.00000 0.00000 0.00000 
+4.00000 -1.00000 0.00000 
+
+
+4.00000 0.00000 0.00000 
+4.00000 -1.00000 0.00000 
+
+5.00000 0.00000 0.00000 
+5.00000 -1.00000 0.00000 
+
+
+5.00000 0.00000 0.00000 
+5.00000 -1.00000 0.00000 
+
+6.00000 0.00000 0.00000 
+6.00000 -1.00000 0.00000 
+
+
+6.00000 0.00000 0.00000 
+6.00000 -1.00000 0.00000 
+
+7.00000 0.00000 0.00000 
+7.00000 -1.00000 0.00000 
+
+
+7.00000 0.00000 0.00000 
+7.00000 -1.00000 0.00000 
+
+8.00000 0.00000 0.00000 
+8.00000 -1.00000 0.00000 
+
+
+0.00000 -1.00000 1.00000 
+0.00000 -2.00000 1.00000 
+
+1.00000 -1.00000 1.00000 
+1.00000 -2.00000 1.00000 
+
+
+1.00000 -1.00000 0.00000 
+1.00000 -2.00000 0.00000 
+
+2.00000 -1.00000 0.00000 
+2.00000 -2.00000 0.00000 
+
+
+2.00000 -1.00000 0.00000 
+2.00000 -2.00000 0.00000 
+
+3.00000 -1.00000 0.00000 
+3.00000 -2.00000 0.00000 
+
+
+3.00000 -1.00000 0.00000 
+3.00000 -2.00000 0.00000 
+
+4.00000 -1.00000 0.00000 
+4.00000 -2.00000 0.00000 
+
+
+4.00000 -1.00000 0.00000 
+4.00000 -2.00000 0.00000 
+
+5.00000 -1.00000 0.00000 
+5.00000 -2.00000 0.00000 
+
+
+5.00000 -1.00000 0.00000 
+5.00000 -2.00000 0.00000 
+
+6.00000 -1.00000 0.00000 
+6.00000 -2.00000 0.00000 
+
+
+6.00000 -1.00000 0.00000 
+6.00000 -2.00000 0.00000 
+
+7.00000 -1.00000 0.00000 
+7.00000 -2.00000 0.00000 
+
+
+7.00000 -1.00000 0.00000 
+7.00000 -2.00000 0.00000 
+
+8.00000 -1.00000 0.00000 
+8.00000 -2.00000 0.00000 
+
+
+0.00000 -2.00000 0.00000 
+0.00000 -3.00000 0.00000 
+
+1.00000 -2.00000 0.00000 
+1.00000 -3.00000 0.00000 
+
+
+1.00000 -2.00000 5.00000 
+1.00000 -3.00000 5.00000 
+
+2.00000 -2.00000 5.00000 
+2.00000 -3.00000 5.00000 
+
+
+2.00000 -2.00000 0.00000 
+2.00000 -3.00000 0.00000 
+
+3.00000 -2.00000 0.00000 
+3.00000 -3.00000 0.00000 
+
+
+3.00000 -2.00000 0.00000 
+3.00000 -3.00000 0.00000 
+
+4.00000 -2.00000 0.00000 
+4.00000 -3.00000 0.00000 
+
+
+4.00000 -2.00000 0.00000 
+4.00000 -3.00000 0.00000 
+
+5.00000 -2.00000 0.00000 
+5.00000 -3.00000 0.00000 
+
+
+5.00000 -2.00000 0.00000 
+5.00000 -3.00000 0.00000 
+
+6.00000 -2.00000 0.00000 
+6.00000 -3.00000 0.00000 
+
+
+6.00000 -2.00000 0.00000 
+6.00000 -3.00000 0.00000 
+
+7.00000 -2.00000 0.00000 
+7.00000 -3.00000 0.00000 
+
+
+7.00000 -2.00000 0.00000 
+7.00000 -3.00000 0.00000 
+
+8.00000 -2.00000 0.00000 
+8.00000 -3.00000 0.00000 
+
+
+0.00000 -3.00000 0.00000 
+0.00000 -4.00000 0.00000 
+
+1.00000 -3.00000 0.00000 
+1.00000 -4.00000 0.00000 
+
+
+1.00000 -3.00000 0.00000 
+1.00000 -4.00000 0.00000 
+
+2.00000 -3.00000 0.00000 
+2.00000 -4.00000 0.00000 
+
+
+2.00000 -3.00000 9.00000 
+2.00000 -4.00000 9.00000 
+
+3.00000 -3.00000 9.00000 
+3.00000 -4.00000 9.00000 
+
+
+3.00000 -3.00000 0.00000 
+3.00000 -4.00000 0.00000 
+
+4.00000 -3.00000 0.00000 
+4.00000 -4.00000 0.00000 
+
+
+4.00000 -3.00000 0.00000 
+4.00000 -4.00000 0.00000 
+
+5.00000 -3.00000 0.00000 
+5.00000 -4.00000 0.00000 
+
+
+5.00000 -3.00000 0.00000 
+5.00000 -4.00000 0.00000 
+
+6.00000 -3.00000 0.00000 
+6.00000 -4.00000 0.00000 
+
+
+6.00000 -3.00000 0.00000 
+6.00000 -4.00000 0.00000 
+
+7.00000 -3.00000 0.00000 
+7.00000 -4.00000 0.00000 
+
+
+7.00000 -3.00000 0.00000 
+7.00000 -4.00000 0.00000 
+
+8.00000 -3.00000 0.00000 
+8.00000 -4.00000 0.00000 
+
+