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 ea6d8981ac713c457bd638c5c55bd253133ba193..ab62a4aca28bca729b636ba1dfdb8df31a8e6079 100644 (file)
--- 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 3d01770e18f68e7c363ea113d9a80e0ad1fcd9b5..8333cf21e83e1b9d41526a752b0abd3a978aec55 100644 (file)
--- 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);
+        };
+      }
   }