From: Matthias Maier <tamiko@43-1.org>
Date: Fri, 14 Jul 2017 19:10:50 +0000 (-0500)
Subject: Bugfix: Do not assume a fully functional matrix interface for mg matrices
X-Git-Tag: v9.0.0-rc1~1416^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d07820c8c61050e1ac459e213c2a8efc5f83e5d2;p=dealii.git

Bugfix: Do not assume a fully functional matrix interface for mg matrices
---

diff --git a/include/deal.II/multigrid/mg_block_smoother.h b/include/deal.II/multigrid/mg_block_smoother.h
index 15c780af93..3d5ccd7032 100644
--- a/include/deal.II/multigrid/mg_block_smoother.h
+++ b/include/deal.II/multigrid/mg_block_smoother.h
@@ -194,33 +194,12 @@ MGSmootherBlock<MatrixType, RelaxationType, number>::initialize (const MGMatrixT
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      // Workaround: The matrix objects supplied to this class do not
-      // expose information how to reinitialize a vector. Therefore,
-      // populate vmult and Tvmult by hand...
-
-      matrices[i] = LinearOperator<BlockVector<number>>();
-      matrices[i].vmult = [&m, i](BlockVector<number> &v,
-                                  const BlockVector<number> &u)
-      {
-        m[i].vmult(v, u);
-      };
-      matrices[i].Tvmult = [&m, i](BlockVector<number> &v,
-                                   const BlockVector<number> &u)
-      {
-        m[i].Tvmult(v, u);
-      };
-
-      smoothers[i] = LinearOperator<BlockVector<number>>();
-      smoothers[i].vmult = [&s, i](BlockVector<number> &v,
-                                   const BlockVector<number> &u)
-      {
-        s[i].vmult(v, u);
-      };
-      smoothers[i].Tvmult = [&s, i](BlockVector<number> &v,
-                                    const BlockVector<number> &u)
-      {
-        s[i].Tvmult(v, u);
-      };
+      // Workaround: Unfortunately, not every "m[i]" object has a
+      // rich enough interface to populate reinit_(domain|range)_vector.
+      // Thus, apply an empty LinearOperator exemplar.
+      matrices[i] = linear_operator<BlockVector<number> >(
+                      LinearOperator<BlockVector<number> >(), m[i]);
+      smoothers[i] = linear_operator<BlockVector<number> >(matrices[i], s[i]);
     }
 }
 
diff --git a/include/deal.II/multigrid/mg_coarse.h b/include/deal.II/multigrid/mg_coarse.h
index 04b28f05f2..f086cecc74 100644
--- a/include/deal.II/multigrid/mg_coarse.h
+++ b/include/deal.II/multigrid/mg_coarse.h
@@ -376,8 +376,11 @@ MGCoarseGridLACIteration<SolverType, VectorType>
   :
   solver(&s, typeid(*this).name())
 {
-  matrix = linear_operator<VectorType>(m);
-  precondition = linear_operator<VectorType>(m, p);
+  // Workaround: Unfortunately, not every "m" object has a rich enough
+  // interface to populate reinit_(domain|range)_vector. Thus, supply an
+  // empty LinearOperator exemplar.
+  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
+  precondition = linear_operator<VectorType>(matrix, p);
 }
 
 
@@ -398,8 +401,11 @@ MGCoarseGridLACIteration<SolverType, VectorType>
               const PreconditionerType &p)
 {
   solver = &s;
-  matrix = linear_operator<VectorType>(m);
-  precondition = linear_operator<VectorType>(m, p);
+  // Workaround: Unfortunately, not every "m" object has a rich enough
+  // interface to populate reinit_(domain|range)_vector. Thus, supply an
+  // empty LinearOperator exemplar.
+  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
+  precondition = linear_operator<VectorType>(matrix, p);
 }
 
 
@@ -433,7 +439,10 @@ void
 MGCoarseGridLACIteration<SolverType, VectorType>
 ::set_matrix(const MatrixType &m)
 {
-  matrix = linear_operator<VectorType>(m);
+  // Workaround: Unfortunately, not every "m" object has a rich enough
+  // interface to populate reinit_(domain|range)_vector. Thus, supply an
+  // empty LinearOperator exemplar.
+  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
 }
 
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
diff --git a/include/deal.II/multigrid/mg_matrix.h b/include/deal.II/multigrid/mg_matrix.h
index 8333cf21e8..3f354be286 100644
--- a/include/deal.II/multigrid/mg_matrix.h
+++ b/include/deal.II/multigrid/mg_matrix.h
@@ -185,20 +185,11 @@ namespace mg
     matrices.resize(p.min_level(), p.max_level());
     for (unsigned int level = p.min_level(); level <= p.max_level(); ++level)
       {
-        // Workaround: The matrix objects supplied to this class do not
-        // expose information how to reinitialize a vector. Therefore,
-        // populate vmult and Tvmult by hand...
-        matrices[level] = LinearOperator<VectorType>();
-        matrices[level].vmult = [&p, level](
-                                  VectorType &v, const VectorType &u)
-        {
-          p[level].vmult(v, u);
-        };
-        matrices[level].Tvmult = [&p, level](
-                                   VectorType &v, const VectorType &u)
-        {
-          p[level].Tvmult(v, u);
-        };
+        // Workaround: Unfortunately, not every "p[level]" object has a
+        // rich enough interface to populate reinit_(domain|range)_vector.
+        // Thus, apply an empty LinearOperator exemplar.
+        matrices[level] =
+          linear_operator<VectorType>(LinearOperator<VectorType>(), p[level]);
       }
   }
 
diff --git a/include/deal.II/multigrid/mg_smoother.h b/include/deal.II/multigrid/mg_smoother.h
index 29726dc887..3c577b66ad 100644
--- a/include/deal.II/multigrid/mg_smoother.h
+++ b/include/deal.II/multigrid/mg_smoother.h
@@ -805,7 +805,11 @@ MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = linear_operator<VectorType>(m[i]);
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] =
+        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
       smoothers[i].initialize(m[i], data);
     }
 }
@@ -830,7 +834,11 @@ MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = linear_operator<VectorType>(m[i]);
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] =
+        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
       smoothers[i].initialize(m[i], data[i]);
     }
 }
@@ -852,7 +860,11 @@ MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = &(m[i].block(row, col));
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
+                                                m[i].block(row, col));
       smoothers[i].initialize(m[i].block(row, col), data);
     }
 }
@@ -879,7 +891,11 @@ MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = &(m[i].block(row, col));
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
+                                                m[i].block(row, col));
       smoothers[i].initialize(m[i].block(row, col), data[i]);
     }
 }
@@ -1009,7 +1025,11 @@ MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = linear_operator<VectorType>(m[i]);
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] =
+        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
       smoothers[i].initialize(m[i], data);
     }
 }
@@ -1036,7 +1056,11 @@ MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize
 
   for (unsigned int i=min; i<=max; ++i)
     {
-      matrices[i] = linear_operator<VectorType>(m[i]);
+      // Workaround: Unfortunately, not every "m[i]" object has a rich
+      // enough interface to populate reinit_(domain|range)_vector. Thus,
+      // apply an empty LinearOperator exemplar.
+      matrices[i] =
+        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
       smoothers[i].initialize(m[i], data[i]);
     }
 }