00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_DUNEMATRIX_HEADER_INCLUDED
00022 #define OPM_DUNEMATRIX_HEADER_INCLUDED
00023
00024 #ifdef DUNE_BCRSMATRIX_HH
00025 #error This header must be included before any bcrsmatrix.hh is included (directly or indirectly)
00026 #endif
00027
00028 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00029
00030 #include <Eigen/Eigen>
00031 #include <Eigen/Sparse>
00032
00033 #include <dune/common/fmatrix.hh>
00034 #include <dune/common/version.hh>
00035
00036 #include <dune/istl/bcrsmatrix.hh>
00037
00038 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00039
00040 namespace Opm
00041 {
00042
00043 class DuneMatrix : public Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> >
00044 {
00045 public:
00046 DuneMatrix(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
00047 {
00048
00049 init( rows, cols, ia, ja, sa );
00050 }
00051
00053 DuneMatrix( const Eigen::SparseMatrix<double, Eigen::RowMajor>& matrix )
00054 {
00055
00056 const int rows = matrix.rows();
00057 const int cols = matrix.cols();
00058 const int* ia = matrix.outerIndexPtr();
00059 const int* ja = matrix.innerIndexPtr();
00060 const double* sa = matrix.valuePtr();
00061
00062 init( rows, cols, ia, ja, sa );
00063 }
00064
00065 protected:
00066 void init(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
00067 {
00068 typedef Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> > Super;
00069 typedef Super::block_type block_type;
00070 this->build_mode = Super::unknown;
00071 this->ready = Super::built;
00072 this->n = rows;
00073 this->m = cols;
00074
00075 typedef Super::size_type size_type ;
00076 #if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2, 4, 1)
00077 size_type& nnz = this->nnz_;
00078 std::shared_ptr<size_type>& j = this->j_;
00079 #else
00080 size_type& nnz = this->nnz;
00081 std::shared_ptr<size_type>& j = this->j;
00082 #endif
00083
00084 nnz = ia[rows];
00085
00086 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 5)
00087 this->allocationSize_ = nnz;
00088 #else
00089 this->allocationSize = nnz;
00090 #endif
00091 this->avg = 0;
00092 this->overflowsize = -1.0;
00093
00094
00095
00096 this->a = this->allocator_.allocate(nnz);
00097 static_assert(sizeof(block_type) == sizeof(double), "This constructor requires a block type that is the same as a double.");
00098 std::copy(sa, sa + nnz, reinterpret_cast<double*>(this->a));
00099 j.reset(this->sizeAllocator_.allocate(nnz));
00100 std::copy(ja, ja +nnz, j.get());
00101 this->r = rowAllocator_.allocate(rows);
00102 for (int row = 0; row < rows; ++row) {
00103 this->r[row].set(ia[row+1] - ia[row], this->a + ia[row], j.get() + ia[row]);
00104 }
00105 }
00106 };
00107
00108 }
00109
00110 #endif // OPM_DUNEMATRIX_HEADER_INCLUDED