(*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]);
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()
}
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;
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)
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);
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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.
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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.