UmfpackSolver: Add unit test for MultiTypeBlockMatrix

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