<commit_before>/* Author: <NAME>, University of Heidelberg, 2003  */
/*         <NAME>, University of Heidelberg, 2010 */
/*         <NAME>, Texas A&M University, 2010   */

/*    $Id$       */
/*                                                                */
/*    Copyright (C) 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors                   */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

// parallel geometric multi-grid. work in progress!

// 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/base/utilities.h>

#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/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/solver_gmres.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_out.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/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>

#include <deal.II/base/index_set.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>

#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_precondition.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/mg_constrained_dofs.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>

// This is C++:
#include <fstream>
#include <sstream>

// The last step is as in all
// previous programs:
namespace Step50
{
  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
  {
  public:
    LaplaceProblem (const unsigned int deg);
    void run ();

  private:
    void setup_system ();
    void assemble_system ();
    void assemble_multigrid ();
    void solve ();
    void refine_grid ();
    void output_results (const unsigned int cycle) const;

    parallel::distributed::Triangulation<dim>   triangulation;
    FE_Q<dim>            fe;
    DoFHandler<dim>    mg_dof_handler;

    typedef TrilinosWrappers::SparseMatrix matrix_t;
    typedef TrilinosWrappers::MPI::Vector vector_t;

    SparsityPattern      sparsity_pattern;
    matrix_t system_matrix;

    // We need an additional object for the
    // hanging nodes constraints. They are
    // handed to the transfer object in the
    // multigrid. Since we call a compress
    // inside the multigrid these constraints
    // are not allowed to be inhomogeneous so
    // we store them in different ConstraintMatrix
    // objects.
    ConstraintMatrix     hanging_node_constraints;
    ConstraintMatrix     constraints;

    vector_t       solution;
    vector_t       system_rhs;

    const unsigned int degree;

    // The following four objects are the
    // only additional member variables,
    // compared to step-6. They first three
    // represent the
    // operators that act on individual
    // levels of the multilevel hierarchy,
    // rather than on the finest mesh as do
    // the objects above while the last object
    // stores information about the boundary
    // indices on each level and information
    // about indices lying on a refinement
    // edge between two different refinement
    // levels.
    //
    // 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<matrix_t> mg_matrices;
    MGLevelObject<matrix_t> mg_interface_matrices;
    MGConstrainedDoFs                    mg_constrained_dofs;
  };



  // @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>
  {
  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]);
  }


  // @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)
    :
    triangulation (MPI_COMM_WORLD,Triangulation<dim>::
                   limit_level_difference_at_vertices,
                   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
    fe (degree),
    mg_dof_handler (triangulation),
    degree(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 ()
  {
    std::string filename = "grid-"
      + Utilities::int_to_string(triangulation.n_levels(), 2)
      + "-"
      + Utilities::int_to_string (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD), 3)
      + ".svg";
    
    std::ofstream out (filename.c_str());
    GridOut grid_out;
    GridOutFlags::Svg svg_flags;
    svg_flags.polar_angle = 60;
    svg_flags.label_level_subdomain_id = true;
    svg_flags.convert_level_number_to_height = true;
    grid_out.set_flags(svg_flags);
    grid_out.write_svg (triangulation, out);
    
    mg_dof_handler.distribute_dofs (fe);
    mg_dof_handler.distribute_mg_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();

    for (unsigned int l=0; l<triangulation.n_levels(); ++l)
      deallog << "   " << 'L' << l << ": "
              << mg_dof_handler.n_dofs(l);
    deallog  << std::endl;

    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());
    solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
    system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);

    // 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. 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, hanging_node_constraints);
    DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);

    typename FunctionMap<dim>::type      dirichlet_boundary;
    ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
    dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
    VectorTools::interpolate_boundary_values (static_cast<const DoFHandler<dim>&>(mg_dof_handler),
                                              dirichlet_boundary,
                                              constraints);
    constraints.close ();
    hanging_node_constraints.close ();
    constraints.condense (sparsity_pattern);
    sparsity_pattern.compress();
    system_matrix.reinit (sparsity_pattern);

    // The multigrid constraints have to be
    // initialized. They need to know about
    // the boundary values as well, so we
    // pass the <code>dirichlet_boundary</code>
    // here as well.
    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);
    mg_interface_matrices.clear ();
    mg_matrices.resize(0, n_levels-1);
    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;
        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]);
      }
  }


  // @sect4{LaplaceProblem::assemble_system}

  // The following function assembles the<commit_msg>step-50: better matrix construction<commit_after>/* Author: <NAME>, University of Heidelberg, 2003  */
/*         <NAME>, University of Heidelberg, 2010 */
/*         <NAME>, Texas A&M University, 2010   */

/*    $Id$       */
/*                                                                */
/*    Copyright (C) 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors                   */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

// parallel geometric multi-grid. work in progress!

// 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/base/utilities.h>

#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/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/solver_gmres.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_out.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/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>

#include <deal.II/base/index_set.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>

#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_precondition.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/mg_constrained_dofs.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>

// This is C++:
#include <fstream>
#include <sstream>

// The last step is as in all
// previous programs:
namespace Step50
{
  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
  {
  public:
    LaplaceProblem (const unsigned int deg);
    void run ();

  private:
    void setup_system ();
    void assemble_system ();
    void assemble_multigrid ();
    void solve ();
    void refine_grid ();
    void output_results (const unsigned int cycle) const;

    parallel::distributed::Triangulation<dim>   triangulation;
    FE_Q<dim>            fe;
    DoFHandler<dim>    mg_dof_handler;

    typedef TrilinosWrappers::SparseMatrix matrix_t;
    typedef TrilinosWrappers::MPI::Vector vector_t;

    matrix_t system_matrix;

    // We need an additional object for the
    // hanging nodes constraints. They are
    // handed to the transfer object in the
    // multigrid. Since we call a compress
    // inside the multigrid these constraints
    // are not allowed to be inhomogeneous so
    // we store them in different ConstraintMatrix
    // objects.
    ConstraintMatrix     hanging_node_constraints;
    ConstraintMatrix     constraints;

    vector_t       solution;
    vector_t       system_rhs;

    const unsigned int degree;

    // The following four objects are the
    // only additional member variables,
    // compared to step-6. They first three
    // represent the
    // operators that act on individual
    // levels of the multilevel hierarchy,
    // rather than on the finest mesh as do
    // the objects above while the last object
    // stores information about the boundary
    // indices on each level and information
    // about indices lying on a refinement
    // edge between two different refinement
    // levels.
    //
    // 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<matrix_t> mg_matrices;
    MGLevelObject<matrix_t> mg_interface_matrices;
    MGConstrainedDoFs                    mg_constrained_dofs;
  };



  // @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>
  {
  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]);
  }


  // @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)
    :
    triangulation (MPI_COMM_WORLD,Triangulation<dim>::
                   limit_level_difference_at_vertices,
                   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
    fe (degree),
    mg_dof_handler (triangulation),
    degree(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 ()
  {
    std::string filename = "grid-"
      + Utilities::int_to_string(triangulation.n_levels(), 2)
      + "-"
      + Utilities::int_to_string (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD), 3)
      + ".svg";
    
    std::ofstream out (filename.c_str());
    GridOut grid_out;
    GridOutFlags::Svg svg_flags;
    svg_flags.polar_angle = 60;
    svg_flags.label_level_subdomain_id = true;
    svg_flags.convert_level_number_to_height = true;
    grid_out.set_flags(svg_flags);
    grid_out.write_svg (triangulation, out);
    
    mg_dof_handler.distribute_dofs (fe);
    mg_dof_handler.distribute_mg_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();

    for (unsigned int l=0; l<triangulation.n_levels(); ++l)
      deallog << "   " << 'L' << l << ": "
              << mg_dof_handler.n_dofs(l);
    deallog  << std::endl;


    //solution.reinit (mg_dof_handler.n_dofs());
    //system_rhs.reinit (mg_dof_handler.n_dofs());
    solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
    system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);

    // 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. 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, hanging_node_constraints);
    DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);

    typename FunctionMap<dim>::type      dirichlet_boundary;
    ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
    dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
    VectorTools::interpolate_boundary_values (static_cast<const DoFHandler<dim>&>(mg_dof_handler),
                                              dirichlet_boundary,
                                              constraints);
    constraints.close ();
    hanging_node_constraints.close ();

    CompressedSimpleSparsityPattern csp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern (mg_dof_handler, csp, constraints);
    system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), csp, MPI_COMM_WORLD, true);

    // The multigrid constraints have to be
    // initialized. They need to know about
    // the boundary values as well, so we
    // pass the <code>dirichlet_boundary</code>
    // here as well.
    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);
    mg_interface_matrices.clear ();
    mg_matrices.resize(0, n_levels-1);
    mg_matrices.clear ();

    // 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;
        csp.reinit(mg_dof_handler.n_dofs(level),
                   mg_dof_handler.n_dofs(level));
        MGTools::make_sparsity_pattern(mg_dof_handler, csp, level);

        mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
            mg_dof_handler.locally_owned_mg_dofs(level),
            csp,
            MPI_COMM_WORLD, true);
      }
  }


  // @sect4{LaplaceProblem::assemble_system}

  // The following function assembles the
{'fraction_non_alphanumeric': 0.09648735574194361, 'fraction_numerical': 0.006180235857980702, 'mean_word_length': 2.619194339837955, 'pattern_counts': {'":': 2, '<': 216, '<?xml version=': 0, '>': 176, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1193768', 'n_tokens_mistral': 9826, 'n_tokens_neox': 9225, 'n_words': 3031}