From 1602fce6b9922bb3f598db18dbb6ccf12abe74ec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 26 Sep 2014 15:03:59 +0200
Subject: [PATCH] Add DuneMatrix class.

This is a hack to get a more efficient constructor for dune-istl
matrices from Eigen matrices.
---
 CMakeLists_files.cmake                      |  1 +
 opm/autodiff/DuneMatrix.hpp                 | 73 +++++++++++++++++++++
 opm/autodiff/NewtonIterationBlackoilCPR.cpp |  7 +-
 3 files changed, 80 insertions(+), 1 deletion(-)
 create mode 100644 opm/autodiff/DuneMatrix.hpp

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 4fc90bb46..d3ee985be 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -96,6 +96,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/CPRPreconditioner.hpp
+	opm/autodiff/DuneMatrix.hpp
 	opm/autodiff/GeoProps.hpp
 	opm/autodiff/GridHelpers.hpp
 	opm/autodiff/ImpesTPFAAD.hpp
diff --git a/opm/autodiff/DuneMatrix.hpp b/opm/autodiff/DuneMatrix.hpp
new file mode 100644
index 000000000..821a5b943
--- /dev/null
+++ b/opm/autodiff/DuneMatrix.hpp
@@ -0,0 +1,73 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_DUNEMATRIX_HEADER_INCLUDED
+#define OPM_DUNEMATRIX_HEADER_INCLUDED
+
+#ifdef DUNE_BCRSMATRIX_HH
+#error This header must be included before any bcrsmatrix.hh is included (directly or indirectly)
+#endif
+
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
+
+#include <dune/common/fmatrix.hh>
+// Include matrix header with hackery to make it possible to inherit.
+#define private protected
+#include <dune/istl/bcrsmatrix.hh>
+#undef private
+
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
+
+namespace Opm
+{
+
+    template <class B>
+    class DuneMatrix : public Dune::BCRSMatrix<B>
+    {
+    public:
+        DuneMatrix(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
+            // : build_mode(unknown), ready(built), n(rows), m(cols), nnz(ia[rows]),
+            //   allocationSize(nnz), r(0), a(0),
+            //   avg(0), overflowsize(-1.0)
+        {
+            typedef Dune::BCRSMatrix<B> Super;
+            this->build_mode = Super::unknown;
+            this->ready = Super::built;
+            this->n = rows;
+            this->m = cols;
+            this->nnz = ia[rows];
+            this->allocationSize = this->nnz;
+            this->avg = 0;
+            this->overflowsize = -1.0;
+
+            this->a = new B[this->nnz];
+            static_assert(sizeof(B) == sizeof(double), "This constructor requires a block type that is the same as a double.");
+            std::copy(sa, sa + this->nnz, reinterpret_cast<double*>(this->a));
+            this->j.reset(new typename Super::size_type[this->nnz]);
+            std::copy(ja, ja +this-> nnz, this->j.get());
+            this->r = new typename Super::row_type[rows];
+            for (int row = 0; row < rows; ++row) {
+                this->r[row].set(ia[row+1] - ia[row], this->a + ia[row], this->j.get() + ia[row]);
+            }
+        }
+    };
+
+} // namespace Opm
+
+#endif // OPM_DUNEMATRIX_HEADER_INCLUDED
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index f7d3b2281..f6ec85704 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -19,6 +19,8 @@
 
 #include <config.h>
 
+#include <opm/autodiff/DuneMatrix.hpp>
+
 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -29,7 +31,7 @@
 #include <opm/core/utility/platform_dependent/disable_warnings.h>
 
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
+// #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/operators.hh>
 #include <dune/istl/io.hh>
 #include <dune/istl/owneroverlapcopy.hh>
@@ -465,6 +467,7 @@ namespace Opm
             const int* ia = matrix.outerIndexPtr();
             const int* ja = matrix.innerIndexPtr();
             const double* sa = matrix.valuePtr();
+#if 0
             Mat A(size, size, nonzeros, Mat::row_wise);
             for (Mat::CreateIterator row = A.createbegin(); row != A.createend(); ++row) {
                 const int ri = row.index();
@@ -478,6 +481,8 @@ namespace Opm
                 }
             }
             return A;
+#endif
+            return Opm::DuneMatrix<MatrixBlockType>(size, size, ia, ja, sa);
         }