]> https://gitweb.dealii.org/ - dealii.git/commitdiff
test DG with constraints and remove text
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 21 Aug 2013 19:07:25 +0000 (19:07 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 21 Aug 2013 19:07:25 +0000 (19:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@30380 0785d39b-7218-0410-832d-ea1e28bc413d

tests/multigrid/step-16.cc
tests/multigrid/step-39-02a.cc [new file with mode: 0644]
tests/multigrid/step-39-02a/cmp/generic [new file with mode: 0644]

index 1d881e62b5f41400bb878eee7fc5ce82f101219b..01bd6949bb11e12c8302ecf571d5baeb3b5ebc16 100644 (file)
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
 
-// As discussed in the introduction, most of
-// this program is copied almost verbatim
-// from step-6, which itself is only a slight
-// modification of step-5. Consequently, a
-// significant part of this program is not
-// new if you've read all the material up to
-// step-6, and we won't comment on that part
-// of the functionality that is
-// unchanged. Rather, we will focus on those
-// aspects of the program that have to do
-// with the multigrid functionality which
-// forms the new aspect of this tutorial
-// program.
-
-// @sect3{Include files}
-
-// Again, the first few include files
-// are already known, so we won't
-// comment on them:
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/error_estimator.h>
 
-// These, now, are the include necessary for
-// the multi-level methods. The first two
-// declare classes that allow us to enumerate
-// degrees of freedom not only on the finest
-// mesh level, but also on intermediate
-// levels (that's what the MGDoFHandler class
-// does) as well as allow to access this
-// information (iterators and accessors over
-// these cells).
-//
-// The rest of the include files deals with
-// the mechanics of multigrid as a linear
-// operator (solver or preconditioner).
 #include <deal.II/multigrid/mg_dof_handler.h>
 #include <deal.II/multigrid/multigrid.h>
 #include <deal.II/multigrid/mg_transfer.h>
 #include <deal.II/multigrid/mg_smoother.h>
 #include <deal.II/multigrid/mg_matrix.h>
 
-// This is C++:
 #include <fstream>
 #include <sstream>
 
-// The last step is as in all
-// previous programs:
 using namespace dealii;
 
-
-// @sect3{The <code>LaplaceProblem</code> class template}
-
-// This main class is basically the same
-// class as in step-6. As far as member
-// functions is concerned, the only addition
-// is the <code>assemble_multigrid</code>
-// function that assembles the matrices that
-// correspond to the discrete operators on
-// intermediate levels:
 template <int dim>
 class LaplaceProblem
 {
@@ -140,37 +95,6 @@ private:
 
   const unsigned int degree;
 
-  // The following three objects are the
-  // only additional member variables,
-  // compared to step-6. They represent the
-  // operators that act on individual
-  // levels of the multilevel hierarchy,
-  // rather than on the finest mesh as do
-  // the objects above.
-  //
-  // To facilitate having objects on each
-  // level of a multilevel hierarchy,
-  // deal.II has the MGLevelObject class
-  // template that provides storage for
-  // objects on each level. What we need
-  // here are matrices on each level, which
-  // implies that we also need sparsity
-  // patterns on each level. As outlined in
-  // the @ref mg_paper, the operators
-  // (matrices) that we need are actually
-  // twofold: one on the interior of each
-  // level, and one at the interface
-  // between each level and that part of
-  // the domain where the mesh is
-  // coarser. In fact, we will need the
-  // latter in two versions: for the
-  // direction from coarse to fine mesh and
-  // from fine to coarse. Fortunately,
-  // however, we here have a self-adjoint
-  // problem for which one of these is the
-  // transpose of the other, and so we only
-  // have to build one; we choose the one
-  // from coarse to fine.
   MGLevelObject<SparsityPattern>       mg_sparsity_patterns;
   MGLevelObject<SparseMatrix<double> > mg_matrices;
   MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
@@ -178,13 +102,6 @@ private:
 };
 
 
-
-// @sect3{Nonconstant coefficients}
-
-// The implementation of nonconstant
-// coefficients is copied verbatim
-// from step-5 and step-6:
-
 template <int dim>
 class Coefficient : public Function<dim>
 {
@@ -231,33 +148,6 @@ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
 }
 
 
-// @sect3{The <code>LaplaceProblem</code> class implementation}
-
-// @sect4{LaplaceProblem::LaplaceProblem}
-
-// The constructor is left mostly
-// unchanged. We take the polynomial degree
-// of the finite elements to be used as a
-// constructor argument and store it in a
-// member variable.
-//
-// By convention, all adaptively refined
-// triangulations in deal.II never change by
-// more than one level across a face between
-// cells. For our multigrid algorithms,
-// however, we need a slightly stricter
-// guarantee, namely that the mesh also does
-// not change by more than refinement level
-// across vertices that might connect two
-// cells. In other words, we must prevent the
-// following situation:
-//
-// @image html limit_level_difference_at_vertices.png ""
-//
-// This is achieved by passing the
-// Triangulation::limit_level_difference_at_vertices
-// flag to the constructor of the
-// triangulation class.
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
   :
@@ -269,22 +159,10 @@ LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
 {}
 
 
-
-// @sect4{LaplaceProblem::setup_system}
-
-// The following function extends what the
-// corresponding one in step-6 did. The top
-// part, apart from the additional output,
-// does the same:
 template <int dim>
 void LaplaceProblem<dim>::setup_system ()
 {
   mg_dof_handler.distribute_dofs (fe);
-
-  // Here we output not only the
-  // degrees of freedom on the finest
-  // level, but also in the
-  // multilevel structure
   deallog << "Number of degrees of freedom: "
           << mg_dof_handler.n_dofs();
 
@@ -303,28 +181,6 @@ void LaplaceProblem<dim>::setup_system ()
   solution.reinit (mg_dof_handler.n_dofs());
   system_rhs.reinit (mg_dof_handler.n_dofs());
 
-  // But it starts to be a wee bit different
-  // here, although this still doesn't have
-  // anything to do with multigrid
-  // methods. step-6 took care of boundary
-  // values and hanging nodes in a separate
-  // step after assembling the global matrix
-  // from local contributions. This works,
-  // but the same can be done in a slightly
-  // simpler way if we already take care of
-  // these constraints at the time of copying
-  // local contributions into the global
-  // matrix. To this end, we here do not just
-  // compute the constraints do to hanging
-  // nodes, but also due to zero boundary
-  // conditions. Both kinds of constraints
-  // can be put into the same object
-  // (<code>constraints</code>), and we will
-  // use this set of constraints later on to
-  // help us copy local contributions
-  // correctly into the global linear system
-  // right away, without the need for a later
-  // clean-up stage:
   constraints.clear ();
   hanging_node_constraints.clear ();
   DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
@@ -345,21 +201,6 @@ void LaplaceProblem<dim>::setup_system ()
 
   mg_constrained_dofs.clear();
   mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary);
-  // Now for the things that concern the
-  // multigrid data structures. First, we
-  // resize the multi-level objects to hold
-  // matrices and sparsity patterns for every
-  // level. The coarse level is zero (this is
-  // mandatory right now but may change in a
-  // future revision). Note that these
-  // functions take a complete, inclusive
-  // range here (not a starting index and
-  // size), so the finest level is
-  // <code>n_levels-1</code>.  We first have
-  // to resize the container holding the
-  // SparseMatrix classes, since they have to
-  // release their SparsityPattern before the
-  // can be destroyed upon resizing.
   const unsigned int n_levels = triangulation.n_levels();
 
   mg_interface_matrices.resize(0, n_levels-1);
@@ -368,32 +209,6 @@ void LaplaceProblem<dim>::setup_system ()
   mg_matrices.clear ();
   mg_sparsity_patterns.resize(0, n_levels-1);
 
-  // Now, we have to provide a matrix on each
-  // level. To this end, we first use the
-  // MGTools::make_sparsity_pattern function
-  // to first generate a preliminary
-  // compressed sparsity pattern on each
-  // level (see the @ref Sparsity module for
-  // more information on this topic) and then
-  // copy it over to the one we really
-  // want. The next step is to initialize
-  // both kinds of level matrices with these
-  // sparsity patterns.
-  //
-  // It may be worth pointing out that the
-  // interface matrices only have entries for
-  // degrees of freedom that sit at or next
-  // to the interface between coarser and
-  // finer levels of the mesh. They are
-  // therefore even sparser than the matrices
-  // on the individual levels of our
-  // multigrid hierarchy. If we were more
-  // concerned about memory usage (and
-  // possibly the speed with which we can
-  // multiply with these matrices), we should
-  // use separate and different sparsity
-  // patterns for these two kinds of
-  // matrices.
   for (unsigned int level=0; level<n_levels; ++level)
     {
       CompressedSparsityPattern csp;
@@ -409,17 +224,6 @@ void LaplaceProblem<dim>::setup_system ()
 }
 
 
-// @sect4{LaplaceProblem::assemble_system}
-
-// The following function assembles the
-// linear system on the finesh level of the
-// mesh. It is almost exactly the same as in
-// step-6, with the exception that we don't
-// eliminate hanging nodes and boundary
-// values after assembling, but while copying
-// local contributions into the global
-// matrix. This is not only simpler but also
-// more efficient for large problems.
 template <int dim>
 void LaplaceProblem<dim>::assemble_system ()
 {
@@ -475,23 +279,6 @@ void LaplaceProblem<dim>::assemble_system ()
 }
 
 
-// @sect4{LaplaceProblem::assemble_multigrid}
-
-// The next function is the one that builds
-// the linear operators (matrices) that
-// define the multigrid method on each level
-// of the mesh. The integration core is the
-// same as above, but the loop below will go
-// over all existing cells instead of just
-// the active ones, and the results must be
-// entered into the correct matrix. Note also
-// that since we only do multi-level
-// preconditioning, no right-hand side needs
-// to be assembled here.
-//
-// Before we go there, however, we have to
-// take care of a significant amount of book
-// keeping:
 template <int dim>
 void LaplaceProblem<dim>::assemble_multigrid ()
 {
@@ -511,70 +298,11 @@ void LaplaceProblem<dim>::assemble_multigrid ()
   const Coefficient<dim> coefficient;
   std::vector<double>    coefficient_values (n_q_points);
 
-  // Next a few things that are specific to
-  // building the multigrid data structures
-  // (since we only need them in the current
-  // function, rather than also elsewhere, we
-  // build them here instead of the
-  // <code>setup_system</code>
-  // function). Some of the following may be
-  // a bit obscure if you're not familiar
-  // with the algorithm actually implemented
-  // in deal.II to support multilevel
-  // algorithms on adaptive meshes; if some
-  // of the things below seem strange, take a
-  // look at the @ref mg_paper.
-  //
-  // Our first job is to identify those
-  // degrees of freedom on each level that
-  // are located on interfaces between
-  // adaptively refined levels, and those
-  // that lie on the interface but also on
-  // the exterior boundary of the domain. As
-  // in many other parts of the library, we
-  // do this by using boolean masks,
-  // i.e. vectors of booleans each element of
-  // which indicates whether the
-  // corresponding degree of freedom index is
-  // an interface DoF or not:
   std::vector<std::vector<bool> > interface_dofs
     = mg_constrained_dofs.get_refinement_edge_indices ();
   std::vector<std::vector<bool> > boundary_interface_dofs
     = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
 
-
-  // The indices just identified will later
-  // be used to impose zero boundary
-  // conditions for the operator that we will
-  // apply on each level. On the other hand,
-  // we also have to impose zero boundary
-  // conditions on the external boundary of
-  // each level. So let's identify these
-  // nodes as well (this time as a set of
-  // degrees of freedom, rather than a
-  // boolean mask; the reason for this being
-  // that we will not need fast tests whether
-  // a certain degree of freedom is in the
-  // boundary list, though we will need such
-  // access for the interface degrees of
-  // freedom further down below):
-
-  // The third step is to construct
-  // constraints on all those degrees of
-  // freedom: their value should be zero
-  // after each application of the level
-  // operators. To this end, we construct
-  // ConstraintMatrix objects for each level,
-  // and add to each of these constraints for
-  // each degree of freedom. Due to the way
-  // the ConstraintMatrix stores its data,
-  // the function to add a constraint on a
-  // single degree of freedom and force it to
-  // be zero is called
-  // Constraintmatrix::add_line(); doing so
-  // for several degrees of freedom at once
-  // can be done using
-  // Constraintmatrix::add_lines():
   std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
   std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
   for (unsigned int level=0; level<triangulation.n_levels(); ++level)
@@ -588,20 +316,6 @@ void LaplaceProblem<dim>::assemble_multigrid ()
       boundary_interface_constraints[level].close ();
     }
 
-  // Now that we're done with most of our
-  // preliminaries, let's start the
-  // integration loop. It looks mostly like
-  // the loop in
-  // <code>assemble_system</code>, with two
-  // exceptions: (i) we don't need a right
-  // han side, and more significantly (ii) we
-  // don't just loop over all active cells,
-  // but in fact all cells, active or
-  // not. Consequently, the correct iterator
-  // to use is MGDoFHandler::cell_iterator
-  // rather than
-  // MGDoFHandler::active_cell_iterator. Let's
-  // go about it:
   typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
                                             endc = mg_dof_handler.end();
 
@@ -621,35 +335,8 @@ void LaplaceProblem<dim>::assemble_multigrid ()
                                  fe_values.shape_grad(j,q_point) *
                                  fe_values.JxW(q_point));
 
-      // The rest of the assembly is again
-      // slightly different. This starts with
-      // a gotcha that is easily forgotten:
-      // The indices of global degrees of
-      // freedom we want here are the ones
-      // for current level, not for the
-      // global matrix. We therefore need the
-      // function
-      // MGDoFAccessorLLget_mg_dof_indices,
-      // not MGDoFAccessor::get_dof_indices
-      // as used in the assembly of the
-      // global system:
       cell->get_mg_dof_indices (local_dof_indices);
 
-      // Next, we need to copy local
-      // contributions into the level
-      // objects. We can do this in the same
-      // way as in the global assembly, using
-      // a constraint object that takes care
-      // of constrained degrees (which here
-      // are only boundary nodes, as the
-      // individual levels have no hanging
-      // node constraints). Note that the
-      // <code>boundary_constraints</code>
-      // object makes sure that the level
-      // matrices contains no contributions
-      // from degrees of freedom at the
-      // interface between cells of different
-      // refinement level.
       boundary_constraints[cell->level()]
       .distribute_local_to_global (cell_matrix,
                                    local_dof_indices,
diff --git a/tests/multigrid/step-39-02a.cc b/tests/multigrid/step-39-02a.cc
new file mode 100644 (file)
index 0000000..6eb5daf
--- /dev/null
@@ -0,0 +1,767 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Add edge matrices and MGConstraints to make sure they are empty and do not mess things up
+
+#include "../tests.h"
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/block_vector.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/multigrid/mg_dof_handler.h>
+
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/assembler.h>
+#include <deal.II/meshworker/loop.h>
+
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_matrix.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <iostream>
+#include <fstream>
+
+namespace Step39
+{
+  using namespace dealii;
+
+  Functions::SlitSingularityFunction<2> exact_solution;
+
+
+
+
+  template <int dim>
+  class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
+  {
+  public:
+    void cell(MeshWorker::DoFInfo<dim> &dinfo,
+              typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void boundary(MeshWorker::DoFInfo<dim> &dinfo,
+                  typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void face(MeshWorker::DoFInfo<dim> &dinfo1,
+              MeshWorker::DoFInfo<dim> &dinfo2,
+              typename MeshWorker::IntegrationInfo<dim> &info1,
+              typename MeshWorker::IntegrationInfo<dim> &info2) const;
+  };
+
+
+  template <int dim>
+  void MatrixIntegrator<dim>::cell(
+    MeshWorker::DoFInfo<dim> &dinfo,
+    typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
+  }
+
+
+  template <int dim>
+  void MatrixIntegrator<dim>::boundary(
+    MeshWorker::DoFInfo<dim> &dinfo,
+    typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
+    LocalIntegrators::Laplace::nitsche_matrix(
+      dinfo.matrix(0,false).matrix, info.fe_values(0),
+      LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg));
+  }
+
+  template <int dim>
+  void MatrixIntegrator<dim>::face(
+    MeshWorker::DoFInfo<dim> &dinfo1,
+    MeshWorker::DoFInfo<dim> &dinfo2,
+    typename MeshWorker::IntegrationInfo<dim> &info1,
+    typename MeshWorker::IntegrationInfo<dim> &info2) const
+  {
+    const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
+    LocalIntegrators::Laplace::ip_matrix(
+      dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix,
+      dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
+      info1.fe_values(0), info2.fe_values(0),
+      LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+  }
+
+  template <int dim>
+  class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
+  {
+  public:
+    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void face(MeshWorker::DoFInfo<dim> &dinfo1,
+              MeshWorker::DoFInfo<dim> &dinfo2,
+              typename MeshWorker::IntegrationInfo<dim> &info1,
+              typename MeshWorker::IntegrationInfo<dim> &info2) const;
+  };
+
+
+  template <int dim>
+  void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &, typename MeshWorker::IntegrationInfo<dim> &) const
+  {}
+
+
+  template <int dim>
+  void RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const FEValuesBase<dim> &fe = info.fe_values();
+    Vector<double> &local_vector = dinfo.vector(0).block(0);
+
+    std::vector<double> boundary_values(fe.n_quadrature_points);
+    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
+
+    const unsigned int deg = fe.get_fe().tensor_degree();
+    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+        local_vector(i) += (- fe.shape_value(i,k) * penalty * boundary_values[k]
+                            + (fe.normal_vector(k) * fe.shape_grad(i,k)) * boundary_values[k])
+                           * fe.JxW(k);
+  }
+
+
+  template <int dim>
+  void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
+                                MeshWorker::DoFInfo<dim> &,
+                                typename MeshWorker::IntegrationInfo<dim> &,
+                                typename MeshWorker::IntegrationInfo<dim> &) const
+  {}
+
+
+  template <int dim>
+  class Estimator : public MeshWorker::LocalIntegrator<dim>
+  {
+  public:
+    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void face(MeshWorker::DoFInfo<dim> &dinfo1,
+              MeshWorker::DoFInfo<dim> &dinfo2,
+              typename MeshWorker::IntegrationInfo<dim> &info1,
+              typename MeshWorker::IntegrationInfo<dim> &info2) const;
+  };
+
+
+  template <int dim>
+  void Estimator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const FEValuesBase<dim> &fe = info.fe_values();
+
+    const std::vector<Tensor<2,dim> > &DDuh = info.hessians[0][0];
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      {
+        const double t = dinfo.cell->diameter() * trace(DDuh[k]);
+        dinfo.value(0) +=  t*t * fe.JxW(k);
+      }
+    dinfo.value(0) = std::sqrt(dinfo.value(0));
+  }
+
+  template <int dim>
+  void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const FEValuesBase<dim> &fe = info.fe_values();
+
+    std::vector<double> boundary_values(fe.n_quadrature_points);
+    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
+
+    const std::vector<double> &uh = info.values[0][0];
+
+    const unsigned int deg = fe.get_fe().tensor_degree();
+    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      dinfo.value(0) += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
+                        * fe.JxW(k);
+    dinfo.value(0) = std::sqrt(dinfo.value(0));
+  }
+
+
+  template <int dim>
+  void Estimator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
+                            MeshWorker::DoFInfo<dim> &dinfo2,
+                            typename MeshWorker::IntegrationInfo<dim> &info1,
+                            typename MeshWorker::IntegrationInfo<dim> &info2) const
+  {
+    const FEValuesBase<dim> &fe = info1.fe_values();
+    const std::vector<double> &uh1 = info1.values[0][0];
+    const std::vector<double> &uh2 = info2.values[0][0];
+    const std::vector<Tensor<1,dim> > &Duh1 = info1.gradients[0][0];
+    const std::vector<Tensor<1,dim> > &Duh2 = info2.gradients[0][0];
+
+    const unsigned int deg = fe.get_fe().tensor_degree();
+    const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+    const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
+    const double penalty = penalty1 + penalty2;
+    const double h = dinfo1.face->measure();
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      {
+        double diff1 = uh1[k] - uh2[k];
+        double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
+        dinfo1.value(0) += (penalty * diff1*diff1 + h * diff2*diff2)
+                           * fe.JxW(k);
+      }
+    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
+    dinfo2.value(0) = dinfo1.value(0);
+  }
+
+
+
+  template <int dim>
+  class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
+  {
+  public:
+    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+    void face(MeshWorker::DoFInfo<dim> &dinfo1,
+              MeshWorker::DoFInfo<dim> &dinfo2,
+              typename MeshWorker::IntegrationInfo<dim> &info1,
+              typename MeshWorker::IntegrationInfo<dim> &info2) const;
+  };
+
+
+  template <int dim>
+  void ErrorIntegrator<dim>::cell(
+    MeshWorker::DoFInfo<dim> &dinfo,
+    typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const FEValuesBase<dim> &fe = info.fe_values();
+    std::vector<Tensor<1,dim> > exact_gradients(fe.n_quadrature_points);
+    std::vector<double> exact_values(fe.n_quadrature_points);
+
+    exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients);
+    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
+
+    const std::vector<Tensor<1,dim> > &Duh = info.gradients[0][0];
+    const std::vector<double> &uh = info.values[0][0];
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      {
+        double sum = 0;
+        for (unsigned int d=0; d<dim; ++d)
+          {
+            const double diff = exact_gradients[k][d] - Duh[k][d];
+            sum += diff*diff;
+          }
+        const double diff = exact_values[k] - uh[k];
+        dinfo.value(0) +=  sum * fe.JxW(k);
+        dinfo.value(1) +=  diff*diff * fe.JxW(k);
+      }
+    dinfo.value(0) = std::sqrt(dinfo.value(0));
+    dinfo.value(1) = std::sqrt(dinfo.value(1));
+  }
+
+
+  template <int dim>
+  void ErrorIntegrator<dim>::boundary(
+    MeshWorker::DoFInfo<dim> &dinfo,
+    typename MeshWorker::IntegrationInfo<dim> &info) const
+  {
+    const FEValuesBase<dim> &fe = info.fe_values();
+
+    std::vector<double> exact_values(fe.n_quadrature_points);
+    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
+
+    const std::vector<double> &uh = info.values[0][0];
+
+    const unsigned int deg = fe.get_fe().tensor_degree();
+    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      {
+        const double diff = exact_values[k] - uh[k];
+        dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
+      }
+    dinfo.value(0) = std::sqrt(dinfo.value(0));
+  }
+
+
+  template <int dim>
+  void ErrorIntegrator<dim>::face(
+    MeshWorker::DoFInfo<dim> &dinfo1,
+    MeshWorker::DoFInfo<dim> &dinfo2,
+    typename MeshWorker::IntegrationInfo<dim> &info1,
+    typename MeshWorker::IntegrationInfo<dim> &info2) const
+  {
+    const FEValuesBase<dim> &fe = info1.fe_values();
+    const std::vector<double> &uh1 = info1.values[0][0];
+    const std::vector<double> &uh2 = info2.values[0][0];
+
+    const unsigned int deg = fe.get_fe().tensor_degree();
+    const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+    const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
+    const double penalty = penalty1 + penalty2;
+
+    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
+      {
+        double diff = uh1[k] - uh2[k];
+        dinfo1.value(0) += (penalty * diff*diff)
+                           * fe.JxW(k);
+      }
+    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
+    dinfo2.value(0) = dinfo1.value(0);
+  }
+
+
+
+
+  template <int dim>
+  class InteriorPenaltyProblem
+  {
+  public:
+    typedef MeshWorker::IntegrationInfo<dim> CellInfo;
+
+    InteriorPenaltyProblem(const FiniteElement<dim> &fe);
+
+    void run(unsigned int n_steps);
+
+  private:
+    void setup_system ();
+    void assemble_matrix ();
+    void assemble_mg_matrix ();
+    void assemble_right_hand_side ();
+    void error ();
+    double estimate ();
+    void solve ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>        triangulation;
+    const MappingQ1<dim>      mapping;
+    const FiniteElement<dim> &fe;
+    MGDoFHandler<dim>         mg_dof_handler;
+    DoFHandler<dim>          &dof_handler;
+      MGConstrainedDoFs mg_constraints;
+      
+    SparsityPattern      sparsity;
+    SparseMatrix<double> matrix;
+    Vector<double>       solution;
+    Vector<double>       right_hand_side;
+    BlockVector<double>  estimates;
+
+    MGLevelObject<SparsityPattern> mg_sparsity;
+    MGLevelObject<SparseMatrix<double> > mg_matrix;
+
+    MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
+    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_down;
+    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_up;
+    MGLevelObject<SparseMatrix<double> > mg_matrix_in_out;
+  };
+
+
+  template <int dim>
+  InteriorPenaltyProblem<dim>::InteriorPenaltyProblem(const FiniteElement<dim> &fe)
+    :
+    mapping(),
+    fe(fe),
+    mg_dof_handler(triangulation),
+    dof_handler(mg_dof_handler),
+    estimates(1)
+  {
+    GridGenerator::hyper_cube_slit(triangulation, -1, 1);
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+    unsigned int n_dofs = dof_handler.n_dofs();
+    solution.reinit(n_dofs);
+    right_hand_side.reinit(n_dofs);
+    
+    mg_constraints.clear();
+    mg_constraints.initialize(dof_handler);
+    
+    CompressedSparsityPattern c_sparsity(n_dofs);
+    DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
+    sparsity.copy_from(c_sparsity);
+    matrix.reinit(sparsity);
+
+    const unsigned int n_levels = triangulation.n_levels();
+    mg_matrix.resize(0, n_levels-1);
+    mg_matrix.clear();
+    mg_matrix_dg_up.resize(0, n_levels-1);
+    mg_matrix_dg_up.clear();
+    mg_matrix_dg_down.resize(0, n_levels-1);
+    mg_matrix_dg_down.clear();
+    mg_matrix_in_out.resize(0, n_levels-1);
+    mg_matrix_in_out.clear();
+    mg_sparsity.resize(0, n_levels-1);
+    mg_sparsity_dg_interface.resize(0, n_levels-1);
+
+    for (unsigned int level=mg_sparsity.min_level();
+         level<=mg_sparsity.max_level(); ++level)
+      {
+        CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level));
+        MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level);
+        mg_sparsity[level].copy_from(c_sparsity);
+        mg_matrix[level].reinit(mg_sparsity[level]);
+        mg_matrix_in_out[level].reinit(mg_sparsity[level]);
+
+        if (level>0)
+          {
+            CompressedSparsityPattern ci_sparsity;
+            ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level));
+            MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level);
+            mg_sparsity_dg_interface[level].copy_from(ci_sparsity);
+            mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
+            mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
+          }
+      }
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::assemble_matrix()
+  {
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags update_flags = update_values | update_gradients;
+    info_box.add_update_flags_all(update_flags);
+    info_box.initialize(fe, mapping);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::MatrixSimple<SparseMatrix<double> > assembler;
+    assembler.initialize(matrix);
+
+    MatrixIntegrator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim>(
+      dof_handler.begin_active(), dof_handler.end(),
+      dof_info, info_box,
+      integrator, assembler);
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::assemble_mg_matrix()
+  {
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags update_flags = update_values | update_gradients;
+    info_box.add_update_flags_all(update_flags);
+    info_box.initialize(fe, mapping);
+
+    MeshWorker::DoFInfo<dim> dof_info(mg_dof_handler);
+
+    MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
+    assembler.initialize(mg_matrix);
+    assembler.initialize(mg_constraints);
+    assembler.initialize_interfaces(mg_matrix_in_out, mg_matrix_in_out);
+    assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
+
+    MatrixIntegrator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim> (
+      mg_dof_handler.begin(), mg_dof_handler.end(),
+      dof_info, info_box,
+      integrator, assembler);
+
+    for (unsigned int level=mg_matrix_in_out.min_level();
+         level<=mg_matrix_in_out.min_level(); ++level)
+      if (mg_matrix_in_out[level].frobenius_norm() != 0.)
+        deallog << "Oops!" << std::endl;
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::assemble_right_hand_side()
+  {
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
+    info_box.add_update_flags_all(update_flags);
+    info_box.initialize(fe, mapping);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::ResidualSimple<Vector<double> > assembler;
+    NamedData<Vector<double>* > data;
+    Vector<double> *rhs = &right_hand_side;
+    data.add(rhs, "RHS");
+    assembler.initialize(data);
+
+    RHSIntegrator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim>(
+      dof_handler.begin_active(), dof_handler.end(),
+      dof_info, info_box,
+      integrator, assembler);
+
+    right_hand_side *= -1.;
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::solve()
+  {
+    SolverControl control(1000, 1.e-12);
+    SolverCG<Vector<double> > solver(control);
+
+    MGTransferPrebuilt<Vector<double> > mg_transfer;
+    mg_transfer.build_matrices(mg_dof_handler);
+
+    FullMatrix<double> coarse_matrix;
+    coarse_matrix.copy_from (mg_matrix[0]);
+    MGCoarseGridHouseholder<double, Vector<double> > mg_coarse;
+    mg_coarse.initialize(coarse_matrix);
+
+    GrowingVectorMemory<Vector<double> > mem;
+    typedef PreconditionSOR<SparseMatrix<double> > RELAXATION;
+    mg::SmootherRelaxation<RELAXATION, Vector<double> >
+    mg_smoother;
+    RELAXATION::AdditionalData smoother_data(1.);
+    mg_smoother.initialize(mg_matrix, smoother_data);
+
+    mg_smoother.set_steps(2);
+    mg_smoother.set_symmetric(true);
+    mg_smoother.set_variable(false);
+
+    MGMatrix<SparseMatrix<double>, Vector<double> > mgmatrix(&mg_matrix);
+    MGMatrix<SparseMatrix<double>, Vector<double> > mgdown(&mg_matrix_dg_down);
+    MGMatrix<SparseMatrix<double>, Vector<double> > mgup(&mg_matrix_dg_up);
+    MGMatrix<SparseMatrix<double>, Vector<double> > mgedge(&mg_matrix_in_out);
+
+    Multigrid<Vector<double> > mg(mg_dof_handler, mgmatrix,
+                                  mg_coarse, mg_transfer,
+                                  mg_smoother, mg_smoother);
+    mg.set_edge_flux_matrices(mgdown, mgup);
+    mg.set_edge_matrices(mgedge, mgedge);
+
+    PreconditionMG<dim, Vector<double>,
+                   MGTransferPrebuilt<Vector<double> > >
+                   preconditioner(mg_dof_handler, mg, mg_transfer);
+    solver.solve(matrix, solution, right_hand_side, preconditioner);
+  }
+
+
+  template <int dim>
+  double
+  InteriorPenaltyProblem<dim>::estimate()
+  {
+    std::vector<unsigned int> old_user_indices;
+    triangulation.save_user_indices(old_user_indices);
+
+    estimates.block(0).reinit(triangulation.n_active_cells());
+    unsigned int i=0;
+    for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell,++i)
+      cell->set_user_index(i);
+
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+    info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
+
+    NamedData<Vector<double>* > solution_data;
+    solution_data.add(&solution, "solution");
+
+    info_box.cell_selector.add("solution", false, false, true);
+    info_box.boundary_selector.add("solution", true, true, false);
+    info_box.face_selector.add("solution", true, true, false);
+
+    info_box.add_update_flags_boundary(update_quadrature_points);
+    info_box.initialize(fe, mapping, solution_data);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::CellsAndFaces<double> assembler;
+    NamedData<BlockVector<double>* > out_data;
+    BlockVector<double> *est = &estimates;
+    out_data.add(est, "cells");
+    assembler.initialize(out_data, false);
+
+    Estimator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim> (
+      dof_handler.begin_active(), dof_handler.end(),
+      dof_info, info_box,
+      integrator, assembler);
+
+    triangulation.load_user_indices(old_user_indices);
+    return estimates.block(0).l2_norm();
+  }
+
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::error()
+  {
+    BlockVector<double> errors(2);
+    errors.block(0).reinit(triangulation.n_active_cells());
+    errors.block(1).reinit(triangulation.n_active_cells());
+    unsigned int i=0;
+    for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell,++i)
+      cell->set_user_index(i);
+
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+    info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
+
+    NamedData<Vector<double>* > solution_data;
+    solution_data.add(&solution, "solution");
+
+    info_box.cell_selector.add("solution", true, true, false);
+    info_box.boundary_selector.add("solution", true, false, false);
+    info_box.face_selector.add("solution", true, false, false);
+
+    info_box.add_update_flags_cell(update_quadrature_points);
+    info_box.add_update_flags_boundary(update_quadrature_points);
+    info_box.initialize(fe, mapping, solution_data);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::CellsAndFaces<double> assembler;
+    NamedData<BlockVector<double>* > out_data;
+    BlockVector<double> *est = &errors;
+    out_data.add(est, "cells");
+    assembler.initialize(out_data, false);
+
+    ErrorIntegrator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim> (
+      dof_handler.begin_active(), dof_handler.end(),
+      dof_info, info_box,
+      integrator, assembler);
+
+    deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl;
+    deallog << "L2-error:     " << errors.block(1).l2_norm() << std::endl;
+  }
+
+
+  template <int dim>
+  void InteriorPenaltyProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    char *fn = new char[100];
+    sprintf(fn, "step-39-02/sol-%02d", cycle);
+
+    std::string filename(fn);
+    filename += ".gnuplot";
+    deallog << "Writing solution to <" << filename << ">..."
+            << std::endl << std::endl;
+    std::ofstream gnuplot_output (filename.c_str());
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "u");
+    data_out.add_data_vector (estimates.block(0), "est");
+
+    data_out.build_patches ();
+
+    data_out.write_gnuplot(gnuplot_output);
+  }
+
+  template <int dim>
+  void
+  InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
+  {
+    deallog << "Element: " << fe.get_name() << std::endl;
+    for (unsigned int s=0; s<n_steps; ++s)
+      {
+        deallog << "Step " << s << std::endl;
+        if (estimates.block(0).size() == 0)
+          triangulation.refine_global(1);
+        else
+          {
+            GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
+                                                               estimates.block(0),
+                                                               0.5, 0.0);
+            triangulation.execute_coarsening_and_refinement ();
+          }
+
+        deallog << "Triangulation "
+                << triangulation.n_active_cells() << " cells, "
+                << triangulation.n_levels() << " levels" << std::endl;
+
+        setup_system();
+        deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, level dofs";
+        for (unsigned int l=0; l<triangulation.n_levels(); ++l)
+          deallog << ' ' << mg_dof_handler.n_dofs(l);
+        deallog << std::endl;
+
+        deallog << "Assemble matrix" << std::endl;
+        assemble_matrix();
+        deallog << "Assemble multilevel matrix" << std::endl;
+        assemble_mg_matrix();
+        deallog << "Assemble right hand side" << std::endl;
+        assemble_right_hand_side();
+        deallog << "Solve" << std::endl;
+        solve();
+        error();
+        deallog << "Estimate " << estimate() << std::endl;
+        //output_results(s);
+      }
+  }
+}
+
+
+
+int main()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step39;
+      initlog(__FILE__);
+
+      FE_DGQ<2> fe1(2);
+      InteriorPenaltyProblem<2> test1(fe1);
+      test1.run(6);
+    }
+  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-39-02a/cmp/generic b/tests/multigrid/step-39-02a/cmp/generic
new file mode 100644 (file)
index 0000000..0af7215
--- /dev/null
@@ -0,0 +1,74 @@
+
+DEAL::Element: FE_DGQ<2>(2)
+DEAL::Step 0
+DEAL::Triangulation 16 cells, 2 levels
+DEAL::DoFHandler 144 dofs, level dofs 36 144
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 23.0730
+DEAL:cg::Convergence step 11 value 5.78162e-14
+DEAL::energy-error: 0.439211
+DEAL::L2-error:     0.0109342
+DEAL::Estimate 0.979555
+DEAL::Step 1
+DEAL::Triangulation 25 cells, 3 levels
+DEAL::DoFHandler 225 dofs, level dofs 36 144 108
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 23.0730
+DEAL:cg::Convergence step 11 value 9.37499e-13
+DEAL::energy-error: 0.332582
+DEAL::L2-error:     0.00548996
+DEAL::Estimate 0.838060
+DEAL::Step 2
+DEAL::Triangulation 37 cells, 4 levels
+DEAL::DoFHandler 333 dofs, level dofs 36 144 144 108
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 23.0730
+DEAL:cg::Convergence step 12 value 3.41649e-13
+DEAL::energy-error: 0.237705
+DEAL::L2-error:     0.00250869
+DEAL::Estimate 0.611449
+DEAL::Step 3
+DEAL::Triangulation 58 cells, 5 levels
+DEAL::DoFHandler 522 dofs, level dofs 36 144 180 180 144
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 23.4898
+DEAL:cg::Convergence step 12 value 5.58516e-13
+DEAL::energy-error: 0.170860
+DEAL::L2-error:     0.00120805
+DEAL::Estimate 0.452418
+DEAL::Step 4
+DEAL::Triangulation 85 cells, 6 levels
+DEAL::DoFHandler 765 dofs, level dofs 36 144 360 180 180 108
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 25.9490
+DEAL:cg::Convergence step 13 value 2.71720e-13
+DEAL::energy-error: 0.122755
+DEAL::L2-error:     0.000602816
+DEAL::Estimate 0.332525
+DEAL::Step 5
+DEAL::Triangulation 130 cells, 7 levels
+DEAL::DoFHandler 1170 dofs, level dofs 36 144 432 432 180 180 144
+DEAL::Assemble matrix
+DEAL::Assemble multilevel matrix
+DEAL::Assemble right hand side
+DEAL::Solve
+DEAL:cg::Starting value 27.1421
+DEAL:cg::Convergence step 13 value 2.50045e-13
+DEAL::energy-error: 0.0869445
+DEAL::L2-error:     0.000292783
+DEAL::Estimate 0.236647

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.