From: Matthias Maier <tamiko@43-1.org>
Date: Fri, 14 Jul 2017 03:08:36 +0000 (-0500)
Subject: Bugfix: Use operator exemplars for preconditioners
X-Git-Tag: v9.0.0-rc1~1423^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4606%2Fhead;p=dealii.git

Bugfix: Use operator exemplars for preconditioners

Further, work around an issue with MG wrappers not exposing vector space
information.
---

diff --git a/include/deal.II/multigrid/mg_coarse.h b/include/deal.II/multigrid/mg_coarse.h
index ea6d8981ac..ab62a4aca2 100644
--- a/include/deal.II/multigrid/mg_coarse.h
+++ b/include/deal.II/multigrid/mg_coarse.h
@@ -377,7 +377,7 @@ MGCoarseGridLACIteration<SolverType, VectorType>
   solver(&s, typeid(*this).name())
 {
   matrix = linear_operator<VectorType>(m);
-  precondition = linear_operator<VectorType>(p);
+  precondition = linear_operator<VectorType>(m, p);
 }
 
 
@@ -398,12 +398,8 @@ MGCoarseGridLACIteration<SolverType, VectorType>
               const PreconditionerType &p)
 {
   solver = &s;
-  if (matrix)
-    delete matrix;
   matrix = linear_operator<VectorType>(m);
-  if (precondition)
-    delete precondition;
-  precondition = linear_operator<VectorType>(p);
+  precondition = linear_operator<VectorType>(m, p);
 }
 
 
@@ -437,8 +433,6 @@ void
 MGCoarseGridLACIteration<SolverType, VectorType>
 ::set_matrix(const MatrixType &m)
 {
-  if (matrix)
-    delete matrix;
   matrix = linear_operator<VectorType>(m);
 }
 
diff --git a/include/deal.II/multigrid/mg_matrix.h b/include/deal.II/multigrid/mg_matrix.h
index 3d01770e18..8333cf21e8 100644
--- a/include/deal.II/multigrid/mg_matrix.h
+++ b/include/deal.II/multigrid/mg_matrix.h
@@ -183,8 +183,23 @@ namespace mg
   Matrix<VectorType>::initialize (const MGLevelObject<MatrixType> &p)
   {
     matrices.resize(p.min_level(), p.max_level());
-    for (unsigned int level=p.min_level(); level <= p.max_level(); ++level)
-      matrices[level] = linear_operator<VectorType>(p[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);
+        };
+      }
   }