DuneMatrix.hpp
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3  Copyright 2014 IRIS AS
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_DUNEMATRIX_HEADER_INCLUDED
22 #define OPM_DUNEMATRIX_HEADER_INCLUDED
23 
24 #ifdef DUNE_BCRSMATRIX_HH
25 #error This header must be included before any bcrsmatrix.hh is included (directly or indirectly)
26 #endif
27 
28 #include <opm/common/utility/platform_dependent/disable_warnings.h>
29 
30 #include <Eigen/Eigen>
31 #include <Eigen/Sparse>
32 
33 #include <dune/common/fmatrix.hh>
34 #include <dune/common/version.hh>
35 
36 #include <dune/istl/bcrsmatrix.hh>
37 
38 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
39 
40 namespace Opm
41 {
42 
43  class DuneMatrix : public Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> >
44  {
45  public:
46  DuneMatrix(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
47  {
48  // create BCRSMatrix from given CSR storage
49  init( rows, cols, ia, ja, sa );
50  }
51 
53  DuneMatrix( const Eigen::SparseMatrix<double, Eigen::RowMajor>& matrix )
54  {
55  // Create ISTL matrix.
56  const int rows = matrix.rows();
57  const int cols = matrix.cols();
58  const int* ia = matrix.outerIndexPtr();
59  const int* ja = matrix.innerIndexPtr();
60  const double* sa = matrix.valuePtr();
61  // create BCRSMatrix from Eigen matrix
62  init( rows, cols, ia, ja, sa );
63  }
64 
65  protected:
66  void init(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
67  {
68  typedef Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> > Super;
69  typedef Super::block_type block_type;
70  this->build_mode = Super::unknown;
71  this->ready = Super::built;
72  this->n = rows;
73  this->m = cols;
74 
75  typedef Super::size_type size_type ;
76 #if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2, 4, 1)
77  size_type& nnz = this->nnz_;
78  std::shared_ptr<size_type>& j = this->j_;
79 #else
80  size_type& nnz = this->nnz;
81  std::shared_ptr<size_type>& j = this->j;
82 #endif
83 
84  nnz = ia[rows];
85 
86 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 5)
87  this->allocationSize_ = nnz;
88 #else
89  this->allocationSize = nnz;
90 #endif
91  this->avg = 0;
92  this->overflowsize = -1.0;
93 
94  // make sure to use the allocators of this matrix
95  // because the same allocators are used to deallocate the data
96  this->a = this->allocator_.allocate(nnz);
97  static_assert(sizeof(block_type) == sizeof(double), "This constructor requires a block type that is the same as a double.");
98  std::copy(sa, sa + nnz, reinterpret_cast<double*>(this->a));
99  j.reset(this->sizeAllocator_.allocate(nnz));
100  std::copy(ja, ja +nnz, j.get());
101  this->r = rowAllocator_.allocate(rows);
102  for (int row = 0; row < rows; ++row) {
103  this->r[row].set(ia[row+1] - ia[row], this->a + ia[row], j.get() + ia[row]);
104  }
105  }
106  };
107 
108 } // namespace Opm
109 
110 #endif // OPM_DUNEMATRIX_HEADER_INCLUDED
DuneMatrix(const Eigen::SparseMatrix< double, Eigen::RowMajor > &matrix)
create an ISTL BCRSMatrix from a Eigen::SparseMatrix
Definition: DuneMatrix.hpp:53
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
Definition: DuneMatrix.hpp:43