Commit 87365c05 authored by Patrick Jaap's avatar Patrick Jaap

istlBackend: create MatrixBuilder for MultiTypeBlockMatrix types

parent 47c20c52
Pipeline #28359 passed with stage
in 12 minutes and 26 seconds
......@@ -7,6 +7,8 @@
- assembleGlobalBasisTransferMatrix:
- Support for all bases in `dune-functions` compatible form added
- `istlMatrixBackend` can now be used with `MultiTypeBlockMatrix`
## Deprecations and removals
- ...
......@@ -6,10 +6,13 @@
#include <dune/common/indices.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/functions/common/indexaccess.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
......@@ -134,6 +137,58 @@ private:
};
/**
* \brief Helper class for building matrix pattern
*
* Version for MultiTypeBlockMatrix with BCRSMatrix-like blocks.
* It uses a 2D std::array of MatrixIndexSets
*/
template<class Row0, class... Rows>
class MatrixBuilder<MultiTypeBlockMatrix<Row0,Rows...>>
{
public:
using Matrix = MultiTypeBlockMatrix<Row0,Rows...>;
MatrixBuilder(Matrix& matrix) :
matrix_(matrix)
{}
template<class RowSizeInfo, class ColSizeInfo>
void resize(const RowSizeInfo& rowSizeInfo, const ColSizeInfo& colSizeInfo)
{
for(size_t i=0; i<rowBlocks_; i++)
for(size_t j=0; j<colBlocks_; j++)
indices_[i][j].resize(rowSizeInfo.size({i}), colSizeInfo.size({j}));
}
template<class RowIndex, class ColIndex>
void insertEntry(const RowIndex& rowIndex, const ColIndex& colIndex)
{
indices_[rowIndex[0]][colIndex[0]].add(rowIndex[1], colIndex[1]);
}
void setupMatrix()
{
Hybrid::forEach(Hybrid::integralRange(std::integral_constant<std::size_t,rowBlocks_>()), [&](auto&& i) {
Hybrid::forEach(Hybrid::integralRange(std::integral_constant<std::size_t,colBlocks_>()), [&](auto&& j) {
indices_[i][j].exportIdx(matrix_[i][j]);
});
});
}
private:
//! number of multi type rows
static constexpr size_t rowBlocks_ = Matrix::N();
//! number of multi type cols (we assume the matrix is well-formed)
static constexpr size_t colBlocks_ = Matrix::M();
//! 2D array of IndexSets
std::array<std::array<Dune::MatrixIndexSet,rowBlocks_>,colBlocks_> indices_;
Matrix& matrix_;
};
template<class M, class E=typename M::field_type>
class ISTLMatrixBackend
......
......@@ -6,8 +6,15 @@
#include <dune/istl/matrix.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/fufem/assemblers/istlbackend.hh>
#include <dune/fufem/test/common.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
// Due to a formerly unnoticed bug, activate bounds checking:
#define DUNE_CHECK_BOUNDS 1
......@@ -25,6 +32,64 @@ TestSuite testISTLBackend(Matrix& mat, RowIndex&& row, ColIndex&& col) {
return suite;
}
template<class Matrix>
TestSuite testMultiType(Matrix& mat) {
TestSuite suite;
auto backend = Fufem::istlMatrixBackend(mat);
auto pattern = backend.patternBuilder();
// create small grid
auto* grid = constructCoarseYaspGrid<3>();
auto gridView = grid->leafGridView();
// create small basis
namespace BB = Functions::BasisBuilder;
// create a composite basis matching the matrix type
const auto basis = BB::makeBasis(
gridView,
composite(
BB::power<3>(
BB::lagrange<1>()
),
BB::power<5>(
BB::lagrange<0>()
)
)
);
pattern.resize(basis,basis);
// set the active entries
auto localView = basis.localView();
for(const auto& e : elements(gridView))
{
localView.bind(e);
for(size_t i=0; i<localView.size(); i++)
{
auto idx = localView.index(i);
pattern.insertEntry(idx,idx);
}
}
pattern.setupMatrix();
mat = 1.0;
for(const auto& e : elements(gridView))
{
localView.bind(e);
for(size_t i=0; i<localView.size(); i++)
{
auto idx = localView.index(i);
suite.check(backend(idx, idx) == 1.0,
"Check matrix access") << "Type" << typeid(mat).name() << " failed";
}
}
return suite;
}
int main(int argc, char** argv) {
MPIHelper::instance(argc, argv);
......@@ -78,5 +143,15 @@ int main(int argc, char** argv) {
suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
}
// MultiTypeMatrix
{
using Matrix = MultiTypeBlockMatrix<MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,3,3>>,BCRSMatrix<FieldMatrix<double,3,5>>>,
MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,5,3>>,BCRSMatrix<FieldMatrix<double,5,5>>>>;
Matrix mat;
suite.subTest(testMultiType(mat));
}
return suite.exit();
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment