Tests: Use proper obstacles for ::ProjectedBlockGSStep

dune/solvers/test/gssteptest.cc

@@ -46,7 +46,7 @@ struct GSTestSuite {
     std::vector<BoxConstraint<double, BlockSize>> obstacles(size);
     for (size_t i = 0; i < size; ++i)
       for (size_t j = 0; j < BlockSize; ++j)
-          obstacles[i][j] = { -1, +1 };
+          obstacles[i][j] = { p.u_ex[i][j] + 2.0, p.u_ex[i][j] + 3.0 };
 
     using LoopSolver = ::LoopSolver<Vector>;
     using Step = LinearIterationStep<Matrix, Vector>;
@@ -182,6 +182,8 @@ struct GSTestSuite {
     // test projected block GS
     if (trivialDirichletOnly) // TODO: missing feature in ProjectedBlockGS
     {
+      // Test 1: Without obstacles, we should get the same solution as with an
+      // unconstrained solver
       size_t size = p.u.size();
       using GSStep = ProjectedBlockGSStep<Matrix, Vector>;
       typename GSStep::HasObstacle hasObstacle(size, false);
@@ -191,11 +193,50 @@ struct GSTestSuite {
         test(&gsStep, "ProjectedBlockGS free");
       }
       hasObstacle.setAll();
+      // Test2: Different initial iterates should lead to the same solution,
+      // which should, furthermore, not intersect the obstacles
       {
-        ProjectedBlockGSStep<Matrix, Vector> gsStep;
-        gsStep.hasObstacle_ = &hasObstacle;
-        gsStep.obstacles_ = &obstacles;
-        test(&gsStep, "ProjectedBlockGS obstacle");
+        ProjectedBlockGSStep<Matrix, Vector> pGSStep1;
+        pGSStep1.hasObstacle_ = &hasObstacle;
+        pGSStep1.obstacles_ = &obstacles;
+        ProjectedBlockGSStep<Matrix, Vector> pGSStep2;
+        pGSStep2.hasObstacle_ = &hasObstacle;
+        pGSStep2.obstacles_ = &obstacles;
+
+        Vector u_copy1 = p.u;
+        for (size_t i = 0; i < p.u.size(); ++i)
+          for (size_t j = 0; j < blocksize; ++j)
+            if (p.ignore[i][j])
+              u_copy1[i][j] += 2.25;
+        Vector u_copy2 = p.u;
+        for (size_t i = 0; i < p.u.size(); ++i)
+          for (size_t j = 0; j < blocksize; ++j)
+            if (p.ignore[i][j])
+              u_copy2[i][j] += 2.75;
+
+        auto result1 = solve(&pGSStep1, u_copy1, 1e-12, 2000);
+        auto result2 = solve(&pGSStep2, u_copy2, 1e-12, 2000);
+
+        for (size_t i = 0; i < p.u.size(); ++i)
+          for (size_t j = 0; j < blocksize; ++j)
+            if (std::abs(result1[i][j] -
+                         obstacles[i][j].projectIn(result1[i][j])) > 1e-12 or
+                std::abs(result2[i][j] -
+                         obstacles[i][j].projectIn(result2[i][j])) > 1e-12) {
+              std::cerr << "### error: obstacles not respected!" << std::endl;
+              passed = false;
+              break;
+            }
+        auto normDiff = p.energyNorm.diff(result1, result2);
+        // We cannot expect too much here (hence the low tolerance):
+        // We are solving a potentially large system using nothing but
+        // a smoother. Convergence will be very slow.
+        if (normDiff > 1e-6) {
+          std::cerr << "### error: Different initial iterates lead to "
+                       "different solutions! (diff = "
+                    << normDiff << ")" << std::endl;
+          passed = false;
+        }
       }
     }