]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Work on comments. Fix and test W and F cycles.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 16 Sep 2016 10:09:50 +0000 (12:09 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 28 Sep 2016 06:39:11 +0000 (08:39 +0200)
include/deal.II/multigrid/multigrid.templates.h
tests/multigrid/cycles.output
tests/multigrid/step-16-05.cc [new file with mode: 0644]
tests/multigrid/step-16-05.output [new file with mode: 0644]
tests/multigrid/step-16-06.cc [new file with mode: 0644]
tests/multigrid/step-16-06.output [new file with mode: 0644]

index 8400c635216e268e4d79e36d4c49411bb29995e3..d5694e1cb215c38349097dc9b1369134de0a1e70 100644 (file)
@@ -121,71 +121,61 @@ Multigrid<VectorType>::level_v_step (const unsigned int level)
       (*coarse)(level, solution[level], defect[level]);
       return;
     }
+
+  // smoothing of the residual
   if (debug>1)
     deallog << "Smoothing on     level " << level << std::endl;
-  // smoothing of the residual by
-  // modifying s
-//  defect[level].print(std::cout, 2,false);
-//  std::cout<<std::endl;
+  // defect[level].print(std::cout, 2,false);
+  // std::cout<<std::endl;
   pre_smooth->smooth(level, solution[level], defect[level]);
-//  solution[level].print(std::cout, 2,false);
+  // solution[level].print(std::cout, 2,false);
 
   if (debug>2)
     deallog << "Solution norm          " << solution[level].l2_norm()
             << std::endl;
 
+  // compute residual on level, which includes the (CG) edge matrix
   if (debug>1)
     deallog << "Residual on      level " << level << std::endl;
-  // t = A*solution[level]
   matrix->vmult(level, t[level], solution[level]);
-
-  if (debug>2)
-    deallog << "Residual norm          " << t[level].l2_norm()
-            << std::endl;
-//  std::cout<<std::endl;
-//  t[level].print(std::cout, 2,false);
-
+  // std::cout<<std::endl;
+  // t[level].print(std::cout, 2,false);
   if (edge_out != 0)
     {
       edge_out->vmult_add(level, t[level], solution[level]);
       if (debug>2)
         deallog << "Norm     t[" << level << "] " << t[level].l2_norm() << std::endl;
     }
-
   t[level].sadd(-1.0, 1.0, defect[level]);
 
-  // transfer the residual to the coarser level
+  if (debug>2)
+    deallog << "Residual norm          " << t[level].l2_norm()
+            << std::endl;
+
+  // Get the defect on the next coarser level as part of the (DG) edge matrix
+  // and then the main part by the restriction of the transfer
   if (edge_down != 0)
     {
       edge_down->vmult(level, t[level-1], solution[level]);
       defect[level-1] -= t[level-1];
     }
-
   transfer->restrict_and_add(level, defect[level-1], t[level]);
 
   // do recursion
-  solution[level-1] = 0.;
   level_v_step(level-1);
 
-  // reset size of the auxiliary
-  // vector, since it has been
-  // resized in the recursive call to
-  // level_v_step directly above
-  t[level] = 0.;
-
   // do coarse grid correction
   transfer->prolongate(level, t[level], solution[level-1]);
   if (debug>2)
     deallog << "Prolongate norm        " << t[level].l2_norm() << std::endl;
-
   solution[level] += t[level];
 
+  // get in contribution from edge matrices to the defect
   if (edge_in != 0)
     {
       edge_in->Tvmult(level, t[level], solution[level]);
       defect[level] -= t[level];
     }
-
   if (edge_up != 0)
     {
       edge_up->Tvmult(level, t[level], solution[level-1]);
@@ -196,15 +186,10 @@ Multigrid<VectorType>::level_v_step (const unsigned int level)
     deallog << "V-cycle  Defect norm   " << defect[level].l2_norm()
             << std::endl;
 
+  // post-smoothing
   if (debug>1)
     deallog << "Smoothing on     level " << level << std::endl;
-  // post-smoothing
-
-//  std::cout<<std::endl;
-//  defect[level].print(std::cout, 2,false);
   post_smooth->smooth(level, solution[level], defect[level]);
-//  solution[level].print(std::cout, 2,false);
-//  std::cout<<std::endl;
 
   if (debug>2)
     deallog << "Solution norm          " << solution[level].l2_norm()
@@ -236,112 +221,106 @@ Multigrid<VectorType>::level_step(const unsigned int level,
     }
 
   if (debug>0)
-    deallog << cychar << "-cycle entering level " << level << std::endl;
+    deallog << cychar << "-cycle entering level  " << level << std::endl;
 
-  // Not actually the defect yet, but
-  // we do not want to spend yet
-  // another vector.
-  if (level>minlevel)
-    {
-      defect2[level-1] = 0.;
-      transfer->restrict_and_add (level, defect2[level-1], defect2[level]);
-    }
-
-  // We get an update of the defect
-  // from the previous level in t and
-  // from two levels above in
-  // defect2. This is subtracted from
-  // the original defect.
-  t[level].equ(-1.,defect2[level]);
-  t[level] += defect[level];
+  // Combine the defect from the initial copy_to_mg with the one that has come
+  // from the finer level by the transfer
+  defect2[level] += defect[level];
+  defect[level] = 0;
 
   if (debug>2)
-    deallog << cychar << "-cycle defect norm    " << t[level].l2_norm()
+    deallog << cychar << "-cycle defect norm     " << defect2[level].l2_norm()
             << std::endl;
 
   if (level == minlevel)
     {
       if (debug>0)
-        deallog << cychar << "-cycle coarse level   " << level << std::endl;
+        deallog << cychar << "-cycle coarse level    " << level << std::endl;
 
-      (*coarse)(level, solution[level], t[level]);
+      (*coarse)(level, solution[level], defect2[level]);
       return;
     }
+
+  // smoothing of the residual
   if (debug>1)
     deallog << cychar << "-cycle smoothing level " << level << std::endl;
-  // smoothing of the residual by
-  // modifying s
-  pre_smooth->smooth(level, solution[level], t[level]);
+  pre_smooth->smooth(level, solution[level], defect2[level]);
 
   if (debug>2)
-    deallog << cychar << "-cycle solution norm    "
+    deallog << cychar << "-cycle solution norm   "
             << solution[level].l2_norm() << std::endl;
 
+  // compute residual on level, which includes the (CG) edge matrix
   if (debug>1)
-    deallog << cychar << "-cycle residual level   " << level << std::endl;
-  // t = A*solution[level]
+    deallog << cychar << "-cycle residual level  " << level << std::endl;
   matrix->vmult(level, t[level], solution[level]);
-  // make t rhs of lower level The
-  // non-refined parts of the
-  // coarse-level defect already
-  // contain the global defect, the
-  // refined parts its restriction.
   if (edge_out != 0)
     edge_out->vmult_add(level, t[level], solution[level]);
+  t[level].sadd(-1.0, 1.0, defect2[level]);
 
+  if (debug>2)
+    deallog << cychar << "-cycle residual norm   " << t[level].l2_norm()
+            << std::endl;
+
+  // Get the defect on the next coarser level as part of the (DG) edge matrix
+  // and then the main part by the restriction of the transfer
   if (edge_down != 0)
-    edge_down->vmult_add(level, defect2[level-1], solution[level]);
+    edge_down->vmult(level, defect2[level-1], solution[level]);
+  else
+    defect2[level-1] = 0;
 
   transfer->restrict_and_add (level, defect2[level-1], t[level]);
-  // do recursion
-  solution[level-1] = 0.;
-  // Every cycle need one recursion,
-  // the V-cycle, which is included
-  // here for the sake of the
-  // F-cycle, needs only one,
+
+  // Every cycle needs one recursion, the V-cycle, which is included here for
+  // the sake of the F-cycle, needs only one. As opposed to the V-cycle above
+  // where we do not need to set the solution vector to zero (as that has been
+  // done in the init functions), we also need to set the content in the
+  // solution from previous cycles to zero.
+  solution[level-1] = 0;
   level_step(level-1, cycle);
-  // If we solve exactly, then we do
-  // not need a second coarse grid
-  // step.
+
+  // For W and F-cycle, repeat the process on the next coarser level except
+  // for the coarse solver which we invoke just once
   if (level>minlevel+1)
     {
       // while the W-cycle repeats itself
       if (cycle == w_cycle)
         level_step(level-1, cycle);
-      // and the F-cycle does a V-cycle
-      // after an F-cycle.
+
+      // the F-cycle does a V-cycle after an F-cycle.
       else if (cycle == f_cycle)
         level_step(level-1, v_cycle);
     }
 
-  // reset size of the auxiliary
-  // vector, since it has been
-  // resized in the recursive call to
-  // level_v_step directly above
-  t[level] = 0.;
   // do coarse grid correction
   transfer->prolongate(level, t[level], solution[level-1]);
-
   solution[level] += t[level];
 
-
+  // get in contribution from edge matrices to the defect
   if (edge_in != 0)
-    edge_in->Tvmult(level, t[level], solution[level]);
+    {
+      edge_in->Tvmult(level, t[level], solution[level]);
+      defect2[level] -= t[level];
+    }
 
   if (edge_up != 0)
-    edge_up->Tvmult(level, t[level], solution[level-1]);
-
-  t[level].sadd(-1.,-1.,defect2[level]);
-  t[level] += defect[level];
+    {
+      edge_up->Tvmult(level, t[level], solution[level-1]);
+      defect2[level] -= t[level];
+    }
 
   if (debug>2)
-    deallog << cychar << "-cycle  Defect norm    " << t[level].l2_norm()
+    deallog << cychar << "-cycle  Defect norm    " << defect2[level].l2_norm()
             << std::endl;
 
+  // post-smoothing
   if (debug>1)
     deallog << cychar << "-cycle smoothing level " << level << std::endl;
-  // post-smoothing
-  post_smooth->smooth(level, solution[level], t[level]);
+  post_smooth->smooth(level, solution[level], defect2[level]);
+
+  if (debug>2)
+    deallog << cychar << "-cycle solution norm   " << solution[level].l2_norm()
+            << std::endl;
 
   if (debug>1)
     deallog << cychar << "-cycle leaving level   " << level << std::endl;
@@ -353,9 +332,8 @@ template <typename VectorType>
 void
 Multigrid<VectorType>::cycle()
 {
-  // The defect vector has been
-  // initialized by copy_to_mg. Now
-  // adjust the other vectors.
+  // The defect vector has been initialized by copy_to_mg. Now adjust the
+  // other vectors.
   solution.resize(minlevel, maxlevel);
   t.resize(minlevel, maxlevel);
   if (cycle_type != v_cycle)
@@ -381,9 +359,8 @@ template <typename VectorType>
 void
 Multigrid<VectorType>::vcycle()
 {
-  // The defect vector has been
-  // initialized by copy_to_mg. Now
-  // adjust the other vectors.
+  // The defect vector has been initialized by copy_to_mg. Now adjust the
+  // other vectors.
   solution.resize(minlevel, maxlevel);
   t.resize(minlevel, maxlevel);
 
index 4a3ca5cc1dca89825b978119641ab74be687f79b..fa9dab6bb5e5cc0df303e5930d7dddc7776ef7a0 100644 (file)
@@ -47,242 +47,292 @@ DEAL::Smoothing on     level 4
 DEAL::Solution norm          0
 DEAL::V-cycle leaving  level 4
 DEAL::
-DEAL::W-cycle entering level 4
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  4
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 4
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   4
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  4
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 2
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   2
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  2
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 2
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   2
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 2
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   2
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  2
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 2
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 2
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   2
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  2
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 2
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   2
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 2
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   2
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  2
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
-DEAL::W-cycle entering level 1
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  1
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 1
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   1
-DEAL::W-cycle entering level 0
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  1
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  0
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    0
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 1
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   1
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 2
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 4
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   4
 DEAL::
-DEAL::F-cycle entering level 4
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle entering level  4
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 4
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   4
-DEAL::F-cycle entering level 3
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  4
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  3
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 3
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   3
-DEAL::F-cycle entering level 2
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  3
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  2
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 2
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   2
-DEAL::F-cycle entering level 1
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  2
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  1
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 1
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   1
-DEAL::F-cycle entering level 0
-DEAL::F-cycle defect norm    0
-DEAL::F-cycle coarse level   0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  1
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  0
+DEAL::F-cycle defect norm     0
+DEAL::F-cycle coarse level    0
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 1
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   1
-DEAL::V-cycle entering level 1
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle entering level  1
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 1
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   1
-DEAL::V-cycle entering level 0
-DEAL::V-cycle defect norm    0
-DEAL::V-cycle coarse level   0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  1
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  0
+DEAL::V-cycle defect norm     0
+DEAL::V-cycle coarse level    0
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 1
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   1
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 2
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   2
-DEAL::V-cycle entering level 2
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle entering level  2
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 2
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   2
-DEAL::V-cycle entering level 1
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  2
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  1
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 1
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   1
-DEAL::V-cycle entering level 0
-DEAL::V-cycle defect norm    0
-DEAL::V-cycle coarse level   0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  1
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  0
+DEAL::V-cycle defect norm     0
+DEAL::V-cycle coarse level    0
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 1
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   1
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 2
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   2
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 3
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   3
-DEAL::V-cycle entering level 3
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle entering level  3
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 3
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   3
-DEAL::V-cycle entering level 2
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  3
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  2
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 2
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   2
-DEAL::V-cycle entering level 1
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  2
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  1
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 1
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   1
-DEAL::V-cycle entering level 0
-DEAL::V-cycle defect norm    0
-DEAL::V-cycle coarse level   0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  1
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  0
+DEAL::V-cycle defect norm     0
+DEAL::V-cycle coarse level    0
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 1
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   1
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 2
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   2
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 3
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   3
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 4
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   4
 DEAL::V-cycle entering level 5
 DEAL::V-cycle  Defect norm   0
@@ -321,129 +371,155 @@ DEAL::Smoothing on     level 5
 DEAL::Solution norm          0
 DEAL::V-cycle leaving  level 5
 DEAL::
-DEAL::W-cycle entering level 5
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  5
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 5
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   5
-DEAL::W-cycle entering level 4
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  5
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  4
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 4
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   4
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  4
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   2
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   2
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 4
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   4
-DEAL::W-cycle entering level 4
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  4
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 4
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   4
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  4
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   2
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
-DEAL::W-cycle entering level 3
-DEAL::W-cycle defect norm    0
+DEAL::W-cycle entering level  3
+DEAL::W-cycle defect norm     0
 DEAL::W-cycle smoothing level 3
-DEAL::W-cycle solution norm    0
-DEAL::W-cycle residual level   3
-DEAL::W-cycle entering level 2
-DEAL::W-cycle defect norm    0
-DEAL::W-cycle coarse level   2
+DEAL::W-cycle solution norm   0
+DEAL::W-cycle residual level  3
+DEAL::W-cycle residual norm   0
+DEAL::W-cycle entering level  2
+DEAL::W-cycle defect norm     0
+DEAL::W-cycle coarse level    2
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 3
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   3
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 4
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   4
 DEAL::W-cycle  Defect norm    0
 DEAL::W-cycle smoothing level 5
+DEAL::W-cycle solution norm   0
 DEAL::W-cycle leaving level   5
 DEAL::
-DEAL::F-cycle entering level 5
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle entering level  5
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 5
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   5
-DEAL::F-cycle entering level 4
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  5
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  4
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 4
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   4
-DEAL::F-cycle entering level 3
-DEAL::F-cycle defect norm    0
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  4
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  3
+DEAL::F-cycle defect norm     0
 DEAL::F-cycle smoothing level 3
-DEAL::F-cycle solution norm    0
-DEAL::F-cycle residual level   3
-DEAL::F-cycle entering level 2
-DEAL::F-cycle defect norm    0
-DEAL::F-cycle coarse level   2
+DEAL::F-cycle solution norm   0
+DEAL::F-cycle residual level  3
+DEAL::F-cycle residual norm   0
+DEAL::F-cycle entering level  2
+DEAL::F-cycle defect norm     0
+DEAL::F-cycle coarse level    2
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 3
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   3
-DEAL::V-cycle entering level 3
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle entering level  3
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 3
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   3
-DEAL::V-cycle entering level 2
-DEAL::V-cycle defect norm    0
-DEAL::V-cycle coarse level   2
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  3
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  2
+DEAL::V-cycle defect norm     0
+DEAL::V-cycle coarse level    2
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 3
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   3
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 4
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   4
-DEAL::V-cycle entering level 4
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle entering level  4
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 4
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   4
-DEAL::V-cycle entering level 3
-DEAL::V-cycle defect norm    0
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  4
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  3
+DEAL::V-cycle defect norm     0
 DEAL::V-cycle smoothing level 3
-DEAL::V-cycle solution norm    0
-DEAL::V-cycle residual level   3
-DEAL::V-cycle entering level 2
-DEAL::V-cycle defect norm    0
-DEAL::V-cycle coarse level   2
+DEAL::V-cycle solution norm   0
+DEAL::V-cycle residual level  3
+DEAL::V-cycle residual norm   0
+DEAL::V-cycle entering level  2
+DEAL::V-cycle defect norm     0
+DEAL::V-cycle coarse level    2
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 3
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   3
 DEAL::V-cycle  Defect norm    0
 DEAL::V-cycle smoothing level 4
+DEAL::V-cycle solution norm   0
 DEAL::V-cycle leaving level   4
 DEAL::F-cycle  Defect norm    0
 DEAL::F-cycle smoothing level 5
+DEAL::F-cycle solution norm   0
 DEAL::F-cycle leaving level   5
diff --git a/tests/multigrid/step-16-05.cc b/tests/multigrid/step-16-05.cc
new file mode 100644 (file)
index 0000000..d0b5bdd
--- /dev/null
@@ -0,0 +1,509 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Similar to step-16-04 but using a W cycle instead of a V-cycle and no
+// adaptive meshes
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <fstream>
+#include <sstream>
+
+using namespace dealii;
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem (const unsigned int deg);
+  void run ();
+
+private:
+  void setup_system ();
+  void assemble_system ();
+  void assemble_multigrid ();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  DoFHandler<dim>      mg_dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  ConstraintMatrix     hanging_node_constraints;
+  ConstraintMatrix     constraints;
+
+  LinearAlgebra::distributed::Vector<double> solution;
+  LinearAlgebra::distributed::Vector<double> system_rhs;
+
+  const unsigned int degree;
+  const unsigned int min_level;
+
+  MGLevelObject<SparsityPattern>       mg_sparsity_patterns;
+  MGLevelObject<SparseMatrix<double> > mg_matrices;
+  MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
+  MGConstrainedDoFs                    mg_constrained_dofs;
+};
+
+
+template <int dim>
+class Coefficient : public Function<dim>
+{
+public:
+  Coefficient () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+
+  virtual void value_list (const std::vector<Point<dim> > &points,
+                           std::vector<double>            &values,
+                           const unsigned int              component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+                                const unsigned int) const
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+}
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+                                   std::vector<double>            &values,
+                                   const unsigned int              component) const
+{
+  const unsigned int n_points = points.size();
+
+  Assert (values.size() == n_points,
+          ExcDimensionMismatch (values.size(), n_points));
+
+  Assert (component == 0,
+          ExcIndexRange (component, 0, 1));
+
+  for (unsigned int i=0; i<n_points; ++i)
+    values[i] = Coefficient<dim>::value (points[i]);
+}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
+  :
+  triangulation (Triangulation<dim>::
+                 limit_level_difference_at_vertices),
+  fe (degree),
+  mg_dof_handler (triangulation),
+  degree(degree),
+  min_level(2)
+{}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  mg_dof_handler.distribute_dofs(fe);
+  mg_dof_handler.distribute_mg_dofs (fe);
+
+  sparsity_pattern.reinit (mg_dof_handler.n_dofs(),
+                           mg_dof_handler.n_dofs(),
+                           mg_dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern);
+
+  solution.reinit (mg_dof_handler.n_dofs());
+  system_rhs.reinit (mg_dof_handler.n_dofs());
+
+  constraints.clear ();
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints);
+  typename FunctionMap<dim>::type      dirichlet_boundary;
+  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
+  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+  MappingQGeneric<dim> mapping(1);
+  VectorTools::interpolate_boundary_values (mapping,
+                                            mg_dof_handler,
+                                            dirichlet_boundary,
+                                            constraints);
+  constraints.close ();
+  hanging_node_constraints.close ();
+  constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+  system_matrix.reinit (sparsity_pattern);
+
+  mg_constrained_dofs.clear();
+  mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary);
+  const unsigned int n_levels = triangulation.n_levels();
+
+  mg_interface_matrices.resize(min_level, n_levels-1);
+  mg_interface_matrices.clear_elements ();
+  mg_matrices.resize(min_level, n_levels-1);
+  mg_matrices.clear_elements ();
+  mg_sparsity_patterns.resize(min_level, n_levels-1);
+
+  for (unsigned int level=min_level; level<n_levels; ++level)
+    {
+      DynamicSparsityPattern csp;
+      csp.reinit(mg_dof_handler.n_dofs(level),
+                 mg_dof_handler.n_dofs(level));
+      MGTools::make_sparsity_pattern(mg_dof_handler, csp, level);
+
+      mg_sparsity_patterns[level].copy_from (csp);
+
+      mg_matrices[level].reinit(mg_sparsity_patterns[level]);
+      mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]);
+    }
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system ()
+{
+  Vector<double> tmp(system_rhs.size());
+  const QGauss<dim>  quadrature_formula(degree+1);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values    |  update_gradients |
+                           update_quadrature_points  |  update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = mg_dof_handler.begin_active(),
+  endc = mg_dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      fe_values.reinit (cell);
+
+      coefficient.value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              cell_matrix(i,j) += (coefficient_values[q_point] *
+                                   fe_values.shape_grad(i,q_point) *
+                                   fe_values.shape_grad(j,q_point) *
+                                   fe_values.JxW(q_point));
+
+            cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+                            1.0 *
+                            fe_values.JxW(q_point));
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                              local_dof_indices,
+                                              system_matrix, tmp);
+    }
+  for (unsigned int i=0; i<tmp.size(); ++i)
+    system_rhs(i) = tmp(i);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_multigrid ()
+{
+  QGauss<dim>  quadrature_formula(1+degree);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
+  ConstraintMatrix empty_constraints;
+  for (unsigned int level=min_level; level<triangulation.n_levels(); ++level)
+    {
+      boundary_constraints[level].add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
+      boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices(level));
+      boundary_constraints[level].close ();
+    }
+
+  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(min_level),
+                                          endc = mg_dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      fe_values.reinit (cell);
+
+      coefficient.value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            cell_matrix(i,j) += (coefficient_values[q_point] *
+                                 fe_values.shape_grad(i,q_point) *
+                                 fe_values.shape_grad(j,q_point) *
+                                 fe_values.JxW(q_point));
+
+      cell->get_mg_dof_indices (local_dof_indices);
+
+      boundary_constraints[cell->level()]
+      .distribute_local_to_global (cell_matrix,
+                                   local_dof_indices,
+                                   mg_matrices[cell->level()]);
+
+      const unsigned int lvl = cell->level();
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          if (mg_constrained_dofs.at_refinement_edge(lvl, local_dof_indices[i])
+              &&
+              ! mg_constrained_dofs.at_refinement_edge(lvl, local_dof_indices[j])
+              &&
+              (
+                (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+                 &&
+                 !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+                ) // ( !boundary(i) && !boundary(j) )
+                ||
+                (
+                  mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+                  &&
+                  local_dof_indices[i]==local_dof_indices[j]
+                ) // ( boundary(i) && boundary(j) && i==j )
+              )
+             )
+            {
+              // do nothing, so add entries to interface matrix
+            }
+          else
+            {
+              cell_matrix(i,j) = 0;
+            }
+
+
+      empty_constraints
+      .distribute_local_to_global (cell_matrix,
+                                   local_dof_indices,
+                                   mg_interface_matrices[cell->level()]);
+    }
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve ()
+{
+  MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >
+  mg_transfer(hanging_node_constraints, mg_constrained_dofs);
+  mg_transfer.build_matrices(mg_dof_handler);
+
+  SolverControl coarse_solver_control (1000, 1e-10, false, false);
+  SolverCG<LinearAlgebra::distributed::Vector<double> > coarse_solver(coarse_solver_control);
+  PreconditionIdentity id;
+  MGCoarseGridLACIteration<SolverCG<LinearAlgebra::distributed::Vector<double> >,LinearAlgebra::distributed::Vector<double> >
+  coarse_grid_solver(coarse_solver, mg_matrices[min_level], id);
+  deallog << "   Size of coarse grid matrix: " << mg_matrices[min_level].m() << std::endl;
+
+  typedef PreconditionChebyshev<SparseMatrix<double>,LinearAlgebra::distributed::Vector<double> > Smoother;
+  GrowingVectorMemory<LinearAlgebra::distributed::Vector<double> >   vector_memory;
+  MGSmootherPrecondition<SparseMatrix<double>, Smoother, LinearAlgebra::distributed::Vector<double> >
+  mg_smoother;
+  typename Smoother::AdditionalData smoother_data;
+  smoother_data.smoothing_range = 20.;
+  smoother_data.degree = 2;
+  smoother_data.eig_cg_n_iterations = 20;
+  mg_smoother.initialize(mg_matrices, smoother_data);
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_matrix(mg_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_interface_up(mg_interface_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_interface_down(mg_interface_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<double> > mg(mg_matrix,
+                                                            coarse_grid_solver,
+                                                            mg_transfer,
+                                                            mg_smoother,
+                                                            mg_smoother,
+                                                            min_level,
+                                                            triangulation.n_global_levels()-1,
+                                                            Multigrid<LinearAlgebra::distributed::Vector<double> >::w_cycle);
+  mg.set_edge_matrices(mg_interface_down, mg_interface_up);
+
+  PreconditionMG<dim, LinearAlgebra::distributed::Vector<double>, MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> > >
+  preconditioner(mg_dof_handler, mg, mg_transfer);
+
+  SolverControl solver_control (1000, 1e-12);
+  SolverCG<LinearAlgebra::distributed::Vector<double> > cg (solver_control);
+
+  solution = 0;
+
+  cg.solve (system_matrix, solution, system_rhs,
+            preconditioner);
+  constraints.distribute (solution);
+
+  deallog << "   " << solver_control.last_step()
+          << " CG iterations needed to obtain convergence."
+          << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<4; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube (triangulation);
+          triangulation.refine_global (min_level+1);
+        }
+      else
+        triangulation.refine_global(1);
+
+
+      deallog << "   Number of active cells:       "
+              << triangulation.n_active_cells()
+              << std::endl;
+
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+              << mg_dof_handler.n_dofs()
+              << " (by level: ";
+      for (unsigned int level=min_level; level<triangulation.n_levels(); ++level)
+        deallog << mg_dof_handler.n_dofs(level)
+                << (level == triangulation.n_levels()-1
+                    ? ")" : ", ");
+      deallog << std::endl;
+
+      assemble_system ();
+      assemble_multigrid ();
+
+      solve ();
+    }
+}
+
+
+// @sect3{The main() function}
+//
+// This is again the same function as
+// in step-6:
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(4);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
+
+  try
+    {
+
+      LaplaceProblem<2> laplace_problem(1);
+      laplace_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/multigrid/step-16-05.output b/tests/multigrid/step-16-05.output
new file mode 100644 (file)
index 0000000..b7f18bf
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 81 (by level: 25, 81)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.1094
+DEAL:cg::Convergence step 10 value 0
+DEAL::   10 CG iterations needed to obtain convergence.
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       256
+DEAL::   Number of degrees of freedom: 289 (by level: 25, 81, 289)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.05859
+DEAL:cg::Convergence step 10 value 0
+DEAL::   10 CG iterations needed to obtain convergence.
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 1089 (by level: 25, 81, 289, 1089)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.9995
+DEAL:cg::Failure step 20 value 0.1128
+DEAL:cg::Starting value 0.03027
+DEAL:cg::Convergence step 11 value 0
+DEAL::   11 CG iterations needed to obtain convergence.
+DEAL::Cycle 3:
+DEAL::   Number of active cells:       4096
+DEAL::   Number of degrees of freedom: 4225 (by level: 25, 81, 289, 1089, 4225)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.9995
+DEAL:cg::Failure step 20 value 0.1128
+DEAL:cg::Starting value 0.9999
+DEAL:cg::Failure step 20 value 1.273
+DEAL:cg::Starting value 0.01538
+DEAL:cg::Convergence step 11 value 0
+DEAL::   11 CG iterations needed to obtain convergence.
diff --git a/tests/multigrid/step-16-06.cc b/tests/multigrid/step-16-06.cc
new file mode 100644 (file)
index 0000000..8780581
--- /dev/null
@@ -0,0 +1,509 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Similar to step-16-04 but using a multigrid F cycle and without adaptive
+// grid
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <fstream>
+#include <sstream>
+
+using namespace dealii;
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem (const unsigned int deg);
+  void run ();
+
+private:
+  void setup_system ();
+  void assemble_system ();
+  void assemble_multigrid ();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  DoFHandler<dim>      mg_dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  ConstraintMatrix     hanging_node_constraints;
+  ConstraintMatrix     constraints;
+
+  LinearAlgebra::distributed::Vector<double> solution;
+  LinearAlgebra::distributed::Vector<double> system_rhs;
+
+  const unsigned int degree;
+  const unsigned int min_level;
+
+  MGLevelObject<SparsityPattern>       mg_sparsity_patterns;
+  MGLevelObject<SparseMatrix<double> > mg_matrices;
+  MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
+  MGConstrainedDoFs                    mg_constrained_dofs;
+};
+
+
+template <int dim>
+class Coefficient : public Function<dim>
+{
+public:
+  Coefficient () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+
+  virtual void value_list (const std::vector<Point<dim> > &points,
+                           std::vector<double>            &values,
+                           const unsigned int              component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+                                const unsigned int) const
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+}
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+                                   std::vector<double>            &values,
+                                   const unsigned int              component) const
+{
+  const unsigned int n_points = points.size();
+
+  Assert (values.size() == n_points,
+          ExcDimensionMismatch (values.size(), n_points));
+
+  Assert (component == 0,
+          ExcIndexRange (component, 0, 1));
+
+  for (unsigned int i=0; i<n_points; ++i)
+    values[i] = Coefficient<dim>::value (points[i]);
+}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
+  :
+  triangulation (Triangulation<dim>::
+                 limit_level_difference_at_vertices),
+  fe (degree),
+  mg_dof_handler (triangulation),
+  degree(degree),
+  min_level(2)
+{}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  mg_dof_handler.distribute_dofs(fe);
+  mg_dof_handler.distribute_mg_dofs (fe);
+
+  sparsity_pattern.reinit (mg_dof_handler.n_dofs(),
+                           mg_dof_handler.n_dofs(),
+                           mg_dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern);
+
+  solution.reinit (mg_dof_handler.n_dofs());
+  system_rhs.reinit (mg_dof_handler.n_dofs());
+
+  constraints.clear ();
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints);
+  typename FunctionMap<dim>::type      dirichlet_boundary;
+  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
+  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+  MappingQGeneric<dim> mapping(1);
+  VectorTools::interpolate_boundary_values (mapping,
+                                            mg_dof_handler,
+                                            dirichlet_boundary,
+                                            constraints);
+  constraints.close ();
+  hanging_node_constraints.close ();
+  constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+  system_matrix.reinit (sparsity_pattern);
+
+  mg_constrained_dofs.clear();
+  mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary);
+  const unsigned int n_levels = triangulation.n_levels();
+
+  mg_interface_matrices.resize(min_level, n_levels-1);
+  mg_interface_matrices.clear_elements ();
+  mg_matrices.resize(min_level, n_levels-1);
+  mg_matrices.clear_elements ();
+  mg_sparsity_patterns.resize(min_level, n_levels-1);
+
+  for (unsigned int level=min_level; level<n_levels; ++level)
+    {
+      DynamicSparsityPattern csp;
+      csp.reinit(mg_dof_handler.n_dofs(level),
+                 mg_dof_handler.n_dofs(level));
+      MGTools::make_sparsity_pattern(mg_dof_handler, csp, level);
+
+      mg_sparsity_patterns[level].copy_from (csp);
+
+      mg_matrices[level].reinit(mg_sparsity_patterns[level]);
+      mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]);
+    }
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system ()
+{
+  Vector<double> tmp(system_rhs.size());
+  const QGauss<dim>  quadrature_formula(degree+1);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values    |  update_gradients |
+                           update_quadrature_points  |  update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = mg_dof_handler.begin_active(),
+  endc = mg_dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      fe_values.reinit (cell);
+
+      coefficient.value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              cell_matrix(i,j) += (coefficient_values[q_point] *
+                                   fe_values.shape_grad(i,q_point) *
+                                   fe_values.shape_grad(j,q_point) *
+                                   fe_values.JxW(q_point));
+
+            cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+                            1.0 *
+                            fe_values.JxW(q_point));
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                              local_dof_indices,
+                                              system_matrix, tmp);
+    }
+  for (unsigned int i=0; i<tmp.size(); ++i)
+    system_rhs(i) = tmp(i);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_multigrid ()
+{
+  QGauss<dim>  quadrature_formula(1+degree);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
+  ConstraintMatrix empty_constraints;
+  for (unsigned int level=min_level; level<triangulation.n_levels(); ++level)
+    {
+      boundary_constraints[level].add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
+      boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices(level));
+      boundary_constraints[level].close ();
+    }
+
+  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(min_level),
+                                          endc = mg_dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      fe_values.reinit (cell);
+
+      coefficient.value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            cell_matrix(i,j) += (coefficient_values[q_point] *
+                                 fe_values.shape_grad(i,q_point) *
+                                 fe_values.shape_grad(j,q_point) *
+                                 fe_values.JxW(q_point));
+
+      cell->get_mg_dof_indices (local_dof_indices);
+
+      boundary_constraints[cell->level()]
+      .distribute_local_to_global (cell_matrix,
+                                   local_dof_indices,
+                                   mg_matrices[cell->level()]);
+
+      const unsigned int lvl = cell->level();
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          if (mg_constrained_dofs.at_refinement_edge(lvl, local_dof_indices[i])
+              &&
+              ! mg_constrained_dofs.at_refinement_edge(lvl, local_dof_indices[j])
+              &&
+              (
+                (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+                 &&
+                 !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+                ) // ( !boundary(i) && !boundary(j) )
+                ||
+                (
+                  mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+                  &&
+                  local_dof_indices[i]==local_dof_indices[j]
+                ) // ( boundary(i) && boundary(j) && i==j )
+              )
+             )
+            {
+              // do nothing, so add entries to interface matrix
+            }
+          else
+            {
+              cell_matrix(i,j) = 0;
+            }
+
+
+      empty_constraints
+      .distribute_local_to_global (cell_matrix,
+                                   local_dof_indices,
+                                   mg_interface_matrices[cell->level()]);
+    }
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve ()
+{
+  MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >
+  mg_transfer(hanging_node_constraints, mg_constrained_dofs);
+  mg_transfer.build_matrices(mg_dof_handler);
+
+  SolverControl coarse_solver_control (1000, 1e-10, false, false);
+  SolverCG<LinearAlgebra::distributed::Vector<double> > coarse_solver(coarse_solver_control);
+  PreconditionIdentity id;
+  MGCoarseGridLACIteration<SolverCG<LinearAlgebra::distributed::Vector<double> >,LinearAlgebra::distributed::Vector<double> >
+  coarse_grid_solver(coarse_solver, mg_matrices[min_level], id);
+  deallog << "   Size of coarse grid matrix: " << mg_matrices[min_level].m() << std::endl;
+
+  typedef PreconditionChebyshev<SparseMatrix<double>,LinearAlgebra::distributed::Vector<double> > Smoother;
+  GrowingVectorMemory<LinearAlgebra::distributed::Vector<double> >   vector_memory;
+  MGSmootherPrecondition<SparseMatrix<double>, Smoother, LinearAlgebra::distributed::Vector<double> >
+  mg_smoother;
+  typename Smoother::AdditionalData smoother_data;
+  smoother_data.smoothing_range = 20.;
+  smoother_data.degree = 2;
+  smoother_data.eig_cg_n_iterations = 20;
+  mg_smoother.initialize(mg_matrices, smoother_data);
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_matrix(mg_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_interface_up(mg_interface_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> > mg_interface_down(mg_interface_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<double> > mg(mg_matrix,
+                                                            coarse_grid_solver,
+                                                            mg_transfer,
+                                                            mg_smoother,
+                                                            mg_smoother,
+                                                            min_level,
+                                                            triangulation.n_global_levels()-1,
+                                                            Multigrid<LinearAlgebra::distributed::Vector<double> >::f_cycle);
+  mg.set_edge_matrices(mg_interface_down, mg_interface_up);
+
+  PreconditionMG<dim, LinearAlgebra::distributed::Vector<double>, MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> > >
+  preconditioner(mg_dof_handler, mg, mg_transfer);
+
+  SolverControl solver_control (1000, 1e-12);
+  SolverCG<LinearAlgebra::distributed::Vector<double> > cg (solver_control);
+
+  solution = 0;
+
+  cg.solve (system_matrix, solution, system_rhs,
+            preconditioner);
+  constraints.distribute (solution);
+
+  deallog << "   " << solver_control.last_step()
+          << " CG iterations needed to obtain convergence."
+          << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<4; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube (triangulation);
+          triangulation.refine_global (min_level+1);
+        }
+      else
+        triangulation.refine_global (1);
+
+
+      deallog << "   Number of active cells:       "
+              << triangulation.n_active_cells()
+              << std::endl;
+
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+              << mg_dof_handler.n_dofs()
+              << " (by level: ";
+      for (unsigned int level=min_level; level<triangulation.n_levels(); ++level)
+        deallog << mg_dof_handler.n_dofs(level)
+                << (level == triangulation.n_levels()-1
+                    ? ")" : ", ");
+      deallog << std::endl;
+
+      assemble_system ();
+      assemble_multigrid ();
+
+      solve ();
+    }
+}
+
+
+// @sect3{The main() function}
+//
+// This is again the same function as
+// in step-6:
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(4);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
+
+  try
+    {
+
+      LaplaceProblem<2> laplace_problem(1);
+      laplace_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/multigrid/step-16-06.output b/tests/multigrid/step-16-06.output
new file mode 100644 (file)
index 0000000..b7f18bf
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 81 (by level: 25, 81)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.1094
+DEAL:cg::Convergence step 10 value 0
+DEAL::   10 CG iterations needed to obtain convergence.
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       256
+DEAL::   Number of degrees of freedom: 289 (by level: 25, 81, 289)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.05859
+DEAL:cg::Convergence step 10 value 0
+DEAL::   10 CG iterations needed to obtain convergence.
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 1089 (by level: 25, 81, 289, 1089)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.9995
+DEAL:cg::Failure step 20 value 0.1128
+DEAL:cg::Starting value 0.03027
+DEAL:cg::Convergence step 11 value 0
+DEAL::   11 CG iterations needed to obtain convergence.
+DEAL::Cycle 3:
+DEAL::   Number of active cells:       4096
+DEAL::   Number of degrees of freedom: 4225 (by level: 25, 81, 289, 1089, 4225)
+DEAL::   Size of coarse grid matrix: 25
+DEAL:cg::Starting value 0.9798
+DEAL:cg::Convergence step 7 value 0
+DEAL:cg::Starting value 0.9938
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 0.9983
+DEAL:cg::Failure step 20 value 0.001167
+DEAL:cg::Starting value 0.9995
+DEAL:cg::Failure step 20 value 0.1128
+DEAL:cg::Starting value 0.9999
+DEAL:cg::Failure step 20 value 1.273
+DEAL:cg::Starting value 0.01538
+DEAL:cg::Convergence step 11 value 0
+DEAL::   11 CG iterations needed to obtain convergence.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.