diff --git a/dune/solvers/test/umfpacksolvertest.cc b/dune/solvers/test/umfpacksolvertest.cc index c396b26aa60dd68fc6ffeb903b364e387dd767e1..ee2c36dcf57c9e1fb2d7bed0aa3621d96d78a342 100644 --- a/dune/solvers/test/umfpacksolvertest.cc +++ b/dune/solvers/test/umfpacksolvertest.cc @@ -62,6 +62,99 @@ struct UMFPackSolverTestSuite +/** \brief test for the UMFPackSolver class with multitype matrix */ +template <size_t blocksize> +struct UMFPackMultiTypeSolverTestSuite +{ + template <class GridType> + bool check(const GridType& grid) + { + double tol = 1e-11; + + bool passed = true; + + using Problem = + SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>; + Problem p(grid.leafGridView()); + + // create a multitype problem + using BlockRow = MultiTypeBlockVector<typename Problem::Matrix,typename Problem::Matrix>; + using Matrix = MultiTypeBlockMatrix<BlockRow,BlockRow>; + + using Vector = MultiTypeBlockVector<typename Problem::Vector,typename Problem::Vector>; + + Matrix A; + Vector u, rhs; + + using namespace Dune::Indices; + // initialize each block + A[_0][_0] = p.A; + A[_0][_1] = p.A; + A[_1][_0] = p.A; + A[_1][_1] = p.A; + + u[_0] = p.u; + u[_1] = p.u; + rhs[_0] = p.rhs; + rhs[_1] = p.rhs; + + using IgnoreBlock = decltype(p.ignore); + using Ignore = Solvers::TupleVector<IgnoreBlock,IgnoreBlock>; + + Ignore ignore; + ignore[_0] = p.ignore; + ignore[_1] = p.ignore; + + typedef Solvers::UMFPackSolver<Matrix,Vector> Solver; + Solver solver(A,u,rhs); + + solver.setIgnore(ignore); + + // solve blockdiagonal + A[_0][_1] *= 0.0; + A[_1][_0] *= 0.0; + solver.solve(); + + // the solution is ( u_ex, u_ex ) + if (p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) > tol) + { + std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) << std::endl; + passed = false; + } + + // solve a problem more generic problem without ignore field + + A[_0][_1] = p.A; + A[_1][_0] = p.A; + + // make A regular + A[_0][_0] *= 2.0; + + u[_0] = p.u; + u[_1] = p.u; + + // solve problem + Solver solver2(A,u,rhs); + solver2.preprocess(); + solver2.solve(); + + // check residual + auto resisual = rhs; + A.mmv(u,resisual); + + if (p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) > tol) + { + std::cout << "The UMFPackSolver did not produce a satisfactory result. ||residual||=" << p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) << std::endl; + passed = false; + } + + + return passed; + } +}; + + + int main(int argc, char** argv) { Dune::MPIHelper::instance(argc, argv); @@ -71,9 +164,13 @@ int main(int argc, char** argv) UMFPackSolverTestSuite<0> testsuite0; UMFPackSolverTestSuite<1> testsuite1; UMFPackSolverTestSuite<2> testsuite2; + + UMFPackMultiTypeSolverTestSuite<3> testsuite3; + passed = checkWithStandardGrids(testsuite0); passed = checkWithStandardGrids(testsuite1); passed = checkWithStandardGrids(testsuite2); + passed = checkWithStandardGrids(testsuite3); return passed ? 0 : 1; }