From c02d4ffeb8f0bd2ab796d41671ffd30103f4ea51 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Tue, 20 Nov 2012 09:57:17 +0100
Subject: [PATCH] added: MatrixOps::saveAsc

save a matrix as a dense asc file for easy import into e.g. octave
note this is only useful for debugging small systems, as the asc
files grow very big and unmanageable
---
 dune/elasticity/matrixops.hpp      |  6 ++++++
 dune/elasticity/matrixops_impl.hpp | 30 ++++++++++++++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp
index 03ed906..e2aefc1 100644
--- a/dune/elasticity/matrixops.hpp
+++ b/dune/elasticity/matrixops.hpp
@@ -60,6 +60,12 @@ class MatrixOps {
     //! \param[in] The matrix to extract the diagonal from
     //! \returns M = diag(A)
     static Matrix extractDiagonal(const Matrix& A);
+
+    //! \brief Save a matrix as a dense asc file
+    //! \param[in] A The matrix to save
+    //! \param[in] file File name
+    //! \details This is only useful for debugging as the files grow very big
+    static void saveAsc(const Matrix& A, const std::string& file);
 };
 
 #include "matrixops_impl.hpp"
diff --git a/dune/elasticity/matrixops_impl.hpp b/dune/elasticity/matrixops_impl.hpp
index 7da51d7..a4c6de5 100644
--- a/dune/elasticity/matrixops_impl.hpp
+++ b/dune/elasticity/matrixops_impl.hpp
@@ -158,3 +158,33 @@ Matrix MatrixOps::extractDiagonal(const Matrix& A)
 
   return result;
 }
+
+void MatrixOps::saveAsc(const Matrix& A, const std::string& file)
+{
+  std::ofstream f;
+  f.open(file.c_str());
+  f << "% " << A.N() << " " << A.M() << std::endl;
+  int prevrow=-1;
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end(); ++it) {
+    for (int i=0;i<int(it.index())-prevrow-1;++i) {
+      for (size_t j=0;j<A.M();++j)
+        f << "0 ";
+      f << std::endl;
+    }
+    int prevcol=-1;
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2) {
+      for (int j=0;j<int(it2.index())-prevcol-1;++j)
+        f << "0 ";
+      double val = *it2;
+      f << val << " ";
+      prevcol = it2.index();
+    }
+    for (int j=0;j<int(A.M())-prevcol-1;++j)
+      f << "0 ";
+    prevrow = it.index();
+    f << std::endl;
+  }
+  f.close();
+}