From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 30 Aug 2022 07:33:09 +0000 (+0200)
Subject: Enable SolverIDR for block vectors
X-Git-Tag: v9.5.0-rc1~996^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a67d62d9ada0fb30b1a82e2c28d807523ee27a69;p=dealii.git

Enable SolverIDR for block vectors
---

diff --git a/include/deal.II/lac/solver_idr.h b/include/deal.II/lac/solver_idr.h
index 64bc3c6793..3c6c66daba 100644
--- a/include/deal.II/lac/solver_idr.h
+++ b/include/deal.II/lac/solver_idr.h
@@ -24,6 +24,7 @@
 #include <deal.II/base/subscriptor.h>
 #include <deal.II/base/utilities.h>
 
+#include <deal.II/lac/block_vector_base.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/solver.h>
 #include <deal.II/lac/solver_control.h>
@@ -226,6 +227,53 @@ namespace internal
         }
       return *data[i];
     }
+
+
+
+    template <typename VectorType,
+              std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
+                * = nullptr>
+    unsigned int
+    n_blocks(const VectorType &)
+    {
+      return 1;
+    }
+
+
+
+    template <typename VectorType,
+              std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
+                nullptr>
+    unsigned int
+    n_blocks(const VectorType &vector)
+    {
+      return vector.n_blocks();
+    }
+
+
+
+    template <typename VectorType,
+              std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
+                * = nullptr>
+    VectorType &
+    block(VectorType &vector, const unsigned int b)
+    {
+      AssertDimension(b, 0);
+      (void)b;
+      return vector;
+    }
+
+
+
+    template <typename VectorType,
+              std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
+                nullptr>
+    typename VectorType::BlockType &
+    block(VectorType &vector, const unsigned int b)
+    {
+      return vector.block(b);
+    }
+
   } // namespace SolverIDRImplementation
 } // namespace internal
 
@@ -328,8 +376,13 @@ SolverIDR<VectorType>::solve(const MatrixType &        A,
       VectorType &tmp_q = Q(i, x);
       if (i != 0)
         {
-          for (auto indx : tmp_q.locally_owned_elements())
-            tmp_q(indx) = normal_distribution(rng);
+          for (unsigned int b = 0;
+               b < internal::SolverIDRImplementation::n_blocks(tmp_q);
+               ++b)
+            for (auto indx : internal::SolverIDRImplementation::block(tmp_q, b)
+                               .locally_owned_elements())
+              internal::SolverIDRImplementation::block(tmp_q, b)(indx) =
+                normal_distribution(rng);
           tmp_q.compress(VectorOperation::insert);
         }
       else