Add DuneMatrix class.

This is a hack to get a more efficient constructor for dune-istl
matrices from Eigen matrices.
This commit is contained in:
Atgeirr Flø Rasmussen
2014-09-26 15:03:59 +02:00
parent 7e472d66e6
commit 1602fce6b9
3 changed files with 80 additions and 1 deletions

View File

@@ -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

View File

@@ -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

View File

@@ -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);
}