Skip to content
Snippets Groups Projects
Commit 63427735 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

[test]Extend test of istlMatrixBackend 

This removed some trivial tests and instead introduces
some more systematic tests with more complicated matrix
types and bases. Since the `MatrixBuilder` has not been
specialized for some of the tested matrix types, the
corresponding matrices must be created manually.
parent a435e04f
No related branches found
No related tags found
1 merge request!149[cleanup]Simplify ISTLMatrixBackend entry access by multi-indices
Hide whitespace changes
Inline Side-by-side
dune/fufem/test/istlbackendtest.cc
+ 246
− 95
View file @ 63427735
#include <config.h>
#include <dune/common/test/testsuite.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/indices.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/classname.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/fufem/backends/istlmatrixbackend.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>
#include <dune/fufem/backends/istlmatrixbackend.hh>
#include <dune/fufem/test/common.hh>
// Due to a formerly unnoticed bug, activate bounds checking:
#define DUNE_CHECK_BOUNDS 1
using namespace Dune;
template<class Entry, class Matrix, class RowIndex, class ColIndex>
TestSuite testISTLBackend(Matrix& mat, RowIndex&& row, ColIndex&& col) {
TestSuite suite;
auto backend = Fufem::istlMatrixBackend<Entry>(mat);
// Build matrix that contains entries for all pairs (i,j)
// of basis functions that share an element.
template<class Matrix, class Basis>
Matrix buildMatrixForBasis(const Basis& basis) {
Matrix matrix;
suite.check(backend(row, col) == 1.0, "Check matrix access") << "Type" << typeid(mat).name() << " failed";
auto backend = Dune::Fufem::istlMatrixBackend(matrix);
auto pattern = backend.patternBuilder();
return suite;
pattern.resize(basis,basis);
auto localView = basis.localView();
for(const auto& e : elements(basis.gridView()))
{
localView.bind(e);
for(auto i : Dune::range(localView.size()))
for(auto j : Dune::range(localView.size()))
pattern.insertEntry(localView.index(i), localView.index(j));
}
pattern.setupMatrix();
return matrix;
}
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();
// Test istlMatrixBackend() for existing matrix.
template<class Entry, class Matrix, class Basis>
Dune::TestSuite testISTLMatrixBackendForBasis(Matrix& matrix, const Basis& basis) {
Dune::TestSuite suite(Dune::className(matrix));
// create small basis
namespace BB = Functions::BasisBuilder;
// create a composite basis matching the matrix type
const auto basis = BB::makeBasis(
gridView,
BB::composite(
BB::power<3>(
BB::lagrange<1>()
),
BB::power<5>(
BB::lagrange<0>()
)
)
);
// FieldTraits<MultiTypeBlockMatrix> are not implemented properly
// in dune-2.9. Hence determining the field type automatically
// in istlMatrixBackend fails and we have to provide the type
// manually.
#if DUNE_VERSION_GT(DUNE_COMMON, 2, 9)
auto backend = Dune::Fufem::istlMatrixBackend(matrix);
#else
auto backend = Dune::Fufem::istlMatrixBackend<Entry>(matrix);
#endif
pattern.resize(basis,basis);
backend.assign(1.0);
// set the active entries
auto localView = basis.localView();
for(const auto& e : elements(gridView))
for(const auto& e : elements(basis.gridView()))
{
localView.bind(e);
for(size_t i=0; i<localView.size(); i++)
for(auto i : Dune::range(localView.size()))
{
auto idx = localView.index(i);
pattern.insertEntry(idx,idx);
auto rowIndex = localView.index(i);
for(auto j : Dune::range(localView.size()))
{
auto colIndex = localView.index(j);
suite.check(backend(rowIndex, colIndex) == 1.0, "Multi-index access")
<< "Unexpected entry (" << rowIndex << "," << colIndex << ")";
}
}
}
pattern.setupMatrix();
mat = 1.0;
return suite;
}
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;
// Build matrix for a basis and test istlMatrixBackend() for it.
template<class Entry, class Matrix, class Basis>
Dune::TestSuite testISTLMatrixBackendForBasis(const Basis& basis) {
Matrix matrix = buildMatrixForBasis<Matrix>(basis);
return testISTLMatrixBackendForBasis<Entry>(matrix, basis);
}
int main(int argc, char** argv) {
MPIHelper::instance(argc, argv);
TestSuite suite;
using K = double;
template<class... V>
using MTBV = Dune::MultiTypeBlockVector<V...>;
// Matrix<FieldMatrix<K, 1,1>
{
using Matrix = Matrix<FieldMatrix<K,1,1>>;
Matrix mat(2,2);
mat =1.0;
template<class... V>
using MTBM = Dune::MultiTypeBlockMatrix<V...>;
// check with multiindex with length 1
{
auto idx = std::array<size_t, 1>{{1}};
template<class K, int n, int m>
using FM = Dune::FieldMatrix<K,n,m>;
suite.subTest(testISTLBackend<K>(mat, idx, idx));
template<class M>
using BCRSM = Dune::BCRSMatrix<M>;
// check const:
const auto& const_mat = mat;
suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
}
// check with multiindex with length 2
{
// someone might set a trailing zero to account for the matrix
// character of FieldMatrix<K, 1,1>:
auto idx = std::array<size_t, 2>{{1, 0}};
template<class M>
using SRM = Dune::Fufem::SingleRowMatrix<M>;
suite.subTest(testISTLBackend<K>(mat, idx, idx));
template<class M>
using SCM = Dune::Fufem::SingleColumnMatrix<M>;
// check const:
const auto& const_mat = mat;
suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
}
int main(int argc, char** argv) {
Dune::MPIHelper::instance(argc, argv);
Dune::TestSuite suite;
using K = double;
using namespace Dune::Functions::BasisBuilder;
// Construct grid for checking backend against bases
auto grid = constructCoarseYaspGrid<3>();
auto gridView = grid->leafGridView();
//****************************************************************************
// Test with dense matrix and scalar entries
//****************************************************************************
{
using Matrix = Dune::Matrix<K>;
const auto basis = makeBasis(
gridView,
lagrange<1>()
);
// MatrixBuilder is not specialized for this type yet.
Matrix matrix(basis.size(), basis.size());
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
}
// Matrix<FieldMatrix<K, n,n> n>1
//****************************************************************************
// Test with dense matrix and scalar entries (as FM<K,1,1>
//****************************************************************************
{
constexpr int n = 2;
using Matrix = Matrix<FieldMatrix<K, n, n>>;
Matrix mat(2,2);
mat = 1.0;
using Matrix = Dune::Matrix<FM<K,1,1>>;
const auto basis = makeBasis(
gridView,
lagrange<1>()
);
Matrix matrix(basis.size(), basis.size());
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
}
//****************************************************************************
// Test with sparse matrix and scalar entries
//****************************************************************************
{
using Matrix = BCRSM<K>;
const auto basis = makeBasis(
gridView,
lagrange<1>()
);
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
}
auto idx = std::array<size_t, 2>{{1, 1}};
//****************************************************************************
// Test with sparse matrix and scalar entries (as FM<K,1,1>)
//****************************************************************************
{
using Matrix = BCRSM<FM<K,1,1>>;
const auto basis = makeBasis(
gridView,
lagrange<1>()
);
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
}
suite.subTest(testISTLBackend<K>(mat, idx, idx));
//****************************************************************************
// Test with dense matrix and block entries
//****************************************************************************
{
using Matrix = Dune::Matrix<FM<K,2,2>>;
const auto basis = makeBasis(
gridView,
power<2>(
lagrange<1>()
)
);
// MatrixBuilder is not specialized for this type yet.
Matrix matrix(basis.size(), basis.size());
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
}
// check const:
const auto& const_mat = mat;
suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
//****************************************************************************
// Test with dense matrix and dynamic block entries
//****************************************************************************
{
using Matrix = Dune::Matrix<Dune::Matrix<K>>;
const auto basis = makeBasis(
gridView,
power<2>(
lagrange<1>()
)
);
// MatrixBuilder is not specialized for this type yet.
Matrix matrix(basis.size(), basis.size());
for (auto i: Dune::range(basis.size()))
for (auto j: Dune::range(basis.size()))
matrix[i][j].setSize(basis.size({i}), basis.size({j}));
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
}
//****************************************************************************
// Test with sparse matrix and block entries
//****************************************************************************
{
using Matrix = BCRSM<FM<K,2,2>>;
const auto basis = makeBasis(
gridView,
power<2>(
lagrange<1>()
)
);
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
}
// MultiTypeMatrix
//****************************************************************************
// Test with MTBM and uniform indices
//****************************************************************************
{
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>>>>;
using Matrix =
MTBM<
MTBV< BCRSM<FM<K,3,3>>, BCRSM<FM<K,3,5>> >,
MTBV< BCRSM<FM<K,5,3>>, BCRSM<FM<K,5,5>> >
>;
const auto basis = makeBasis(
gridView,
composite(
power<3>(
lagrange<1>()
),
power<5>(
lagrange<0>()
)
)
);
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
}
Matrix mat;
//****************************************************************************
// Test with MTBM and non-uniform indices
//****************************************************************************
{
using Matrix =
MTBM<
MTBV< BCRSM<FM<K,3,3>>, BCRSM<SCM<FM<K,3,1>>> >,
MTBV< BCRSM<SRM<FM<K,3,1>>>, BCRSM<FM<K,1,1>> >
>;
const auto basis = makeBasis(
gridView,
composite(
power<3>(
lagrange<2>()
),
lagrange<1>()
)
);
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
}
suite.subTest(testMultiType(mat));
//****************************************************************************
// Test with nested MTBM and non-uniform indices (scalar entries as FM<K,1,1>)
//****************************************************************************
{
using Matrix =
MTBM<
MTBV<
MTBM<
MTBV< Dune::Matrix<FM<K,3,3>>, Dune::Matrix<SCM<FM<K,3,1>>> >,
MTBV< Dune::Matrix<SRM<FM<K,3,1>>>, Dune::Matrix<FM<K,1,1>> >
>
>
>;
const auto basis = makeBasis(
gridView,
composite(
composite(
power<3>(
lagrange<2>()
),
lagrange<1>()
)
)
);
// MatrixBuilder is not specialized for this type yet.
using namespace Dune::Indices;
Matrix matrix;
matrix[_0][_0][_0][_0].setSize(basis.size({0,0}), basis.size({0,0}));
matrix[_0][_0][_0][_1].setSize(basis.size({0,0}), basis.size({0,1}));
matrix[_0][_0][_1][_0].setSize(basis.size({0,1}), basis.size({0,0}));
matrix[_0][_0][_1][_1].setSize(basis.size({0,1}), basis.size({0,1}));
suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
}
return suite.exit();
......
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