]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixing problems with MGConstrainedDoFs
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 31 Jul 2014 11:23:16 +0000 (13:23 +0200)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 31 Jul 2014 11:23:16 +0000 (13:23 +0200)
include/deal.II/multigrid/mg_constrained_dofs.h
include/deal.II/multigrid/mg_tools.h
source/multigrid/mg_tools.cc
source/multigrid/mg_tools.inst.in
tests/multigrid/step-16-02.cc
tests/multigrid/step-16-bdry1.cc [new file with mode: 0644]
tests/multigrid/step-16-bdry1.output [new file with mode: 0644]
tests/multigrid/step-16.cc

index 0547380d321612e16e891a9254819a70c3d81eb5..d7cac5aa110ee187284500755ead48c581673dbd 100644 (file)
@@ -42,18 +42,6 @@ class MGConstrainedDoFs : public Subscriptor
 public:
 
   typedef std::vector<std::set<types::global_dof_index> >::size_type size_dof;
-  /**
-   * Fill the internal data structures with values extracted from the dof
-   * handler object.
-   *
-   * This function ensures that on every level, degrees of freedom at interior
-   * edges of a refinement level are treated corrected but leaves degrees of
-   * freedom at the boundary of the domain untouched assuming that no
-   * Dirichlet boundary conditions for them exist.
-   */
-  template <int dim, int spacedim>
-  void initialize (const DoFHandler<dim,spacedim> &dof);
-
   /**
    * Fill the internal data structures with values extracted from the dof
    * handler object and apply the boundary values provided.
@@ -74,18 +62,17 @@ public:
 
   /**
    * Determine whether a dof index is subject to a boundary
-   * constraint. (In other words, whether it is on boundary of domain.)
+   * constraint.
    */
   bool is_boundary_index (const unsigned int level,
                           const types::global_dof_index index) const;
 
   /**
-   * Determine whether a dof index is at an edge that is not a
+   * @deprecated Determine whether a dof index is at an edge that is not a
    * refinement edge.
    */
-// TODO: remove
   bool non_refinement_edge_index (const unsigned int level,
-                                  const types::global_dof_index index) const;
+                                  const types::global_dof_index index) const DEAL_II_DEPRECATED;
 
   /**
    * Determine whether a dof index is at the refinement edge.
@@ -94,16 +81,19 @@ public:
                            const types::global_dof_index index) const;
 
   /**
-   * Determine whether a dof index is at the refinement edge and
-   * subject to a boundary constraint .
-   * = is_boundary_index() && at_refinement_edge()
+   * @deprecated Use is_boundary_index() instead. The logic behind
+   * this function here is unclear and for practical purposes, the
+   * other is needed.
+   *
+   * Determine whether a dof index is subject to a boundary
+   * constraint.
    */
   bool at_refinement_edge_boundary (const unsigned int level,
-                                    const types::global_dof_index index) const;
+                                    const types::global_dof_index index) const DEAL_II_DEPRECATED;
 
   /**
-   * Return the indices of dofs for each level that lie on the
-   * boundary of the domain.
+   * Return the indices of dofs for each level that are subject to
+   * boundary constraints.
    */
   const std::vector<std::set<types::global_dof_index> > &
   get_boundary_indices () const;
@@ -137,13 +127,6 @@ public:
   const std::vector<std::vector<bool> > &
   get_refinement_edge_boundary_indices () const DEAL_II_DEPRECATED;
 
-  /**
-   * Return indices of all dofs that are on boundary faces on the given level 
-   * if the cell has refinement edge indices (i.e. has a coarser neighbor).
-   */
-  const IndexSet &
-  get_refinement_edge_boundary_indices (unsigned int level) const;
-
   /**
    * @deprecated The function is_boundary_index() now returns false if
    * no boundary values are set.
@@ -154,6 +137,21 @@ public:
 
 private:
 
+  /**
+   * @warning This function generates an inconsistent object if not
+   * called from the other initialize() in this class.
+   *
+   * Fill the internal data structures with values extracted from the dof
+   * handler object.
+   *
+   * This function ensures that on every level, degrees of freedom at interior
+   * edges of a refinement level are treated corrected but leaves degrees of
+   * freedom at the boundary of the domain untouched assuming that no
+   * Dirichlet boundary conditions for them exist.
+   */
+  template <int dim, int spacedim>
+  void initialize (const DoFHandler<dim,spacedim> &dof);
+
   /**
    * The indices of boundary dofs for each level.
    */
@@ -165,14 +163,6 @@ private:
    */
   std::vector<IndexSet> refinement_edge_indices;
 
-  /**
-   * The degrees of freedom on the refinement edge between a level and
-   * coarser cells, which are also on the boundary.
-   *
-   * This is a subset of #refinement_edge_indices.
-   */
-  std::vector<IndexSet> refinement_edge_boundary_indices;
-
   /**
    * old data structure only filled on demand
    */
@@ -195,7 +185,6 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim,spacedim> &dof)
   boundary_indices.resize(nlevels);
 
   refinement_edge_indices.resize(nlevels);
-  refinement_edge_boundary_indices.resize(nlevels);
   refinement_edge_indices_old.clear();
   refinement_edge_boundary_indices_old.clear();
   for (unsigned int l=0; l<nlevels; ++l)
@@ -203,11 +192,9 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim,spacedim> &dof)
       boundary_indices[l].clear();
 
       refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
-      refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l));
     }
   
-  MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices,
-                                         refinement_edge_boundary_indices);
+  MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices);
 }
 
 
@@ -237,9 +224,6 @@ MGConstrainedDoFs::clear()
   for (unsigned int l=0; l<refinement_edge_indices.size(); ++l)
     refinement_edge_indices[l].clear();
 
-  for (unsigned int l=0; l<refinement_edge_boundary_indices.size(); ++l)
-    refinement_edge_boundary_indices[l].clear();
-
   refinement_edge_indices_old.clear();
   refinement_edge_boundary_indices_old.clear();
 }
@@ -281,9 +265,7 @@ bool
 MGConstrainedDoFs::at_refinement_edge_boundary (const unsigned int level,
                                                 const types::global_dof_index index) const
 {
-  AssertIndexRange(level, refinement_edge_boundary_indices.size());
-
-  return refinement_edge_boundary_indices[level].is_element(index);
+  return is_boundary_index(level, index);
 }
 
 inline
@@ -319,30 +301,27 @@ MGConstrainedDoFs::get_refinement_edge_indices (unsigned int level) const
   return refinement_edge_indices[level];
 }
 
+
 inline
 const std::vector<std::vector<bool> > &
 MGConstrainedDoFs::get_refinement_edge_boundary_indices () const
 {
-  if (refinement_edge_boundary_indices_old.size()!=refinement_edge_boundary_indices.size())
-    {
-      unsigned int n_levels = refinement_edge_boundary_indices.size();
-      refinement_edge_boundary_indices_old.resize(n_levels);
-      for (unsigned int l=0;l<n_levels;++l)
-        {
-          refinement_edge_boundary_indices_old[l].resize(refinement_edge_boundary_indices[l].size(), false);
-          refinement_edge_boundary_indices[l].fill_binary_vector(refinement_edge_boundary_indices_old[l]);
-        }
-    }
+  Assert(false, ExcMessage("Timo fixes this"));
+  
+  // if (refinement_edge_boundary_indices_old.size()!=refinement_edge_boundary_indices.size())
+  //   {
+  //     unsigned int n_levels = refinement_edge_boundary_indices.size();
+  //     refinement_edge_boundary_indices_old.resize(n_levels);
+  //     for (unsigned int l=0;l<n_levels;++l)
+  //       {
+  //         refinement_edge_boundary_indices_old[l].resize(refinement_edge_boundary_indices[l].size(), false);
+  //         refinement_edge_boundary_indices[l].fill_binary_vector(refinement_edge_boundary_indices_old[l]);
+  //       }
+  //   }
 
   return refinement_edge_boundary_indices_old;
 }
 
-inline
-const IndexSet &
-MGConstrainedDoFs::get_refinement_edge_boundary_indices (unsigned int level) const
-{
-  return refinement_edge_boundary_indices[level];
-}
 
 inline
 bool
index 143bffb74651b74f19fee956dce90d09f3b75f1b..61bc2a3119407d4f4895f09fd0fd64c86db8672e 100644 (file)
@@ -267,31 +267,10 @@ namespace MGTools
                          const bool preserve_symmetry) DEAL_II_DEPRECATED;
 
   /**
-   * For each level in a multigrid
-   * hierarchy, produce a boolean
-   * mask that indicates which of
-   * the degrees of freedom are
-   * along interfaces of this level
-   * to cells that only exist on
-   * coarser levels. The function
-   * returns the subset of these
-   * indices in the last argument
-   * that are not only on interior
-   * interfaces (i.e. between cells
-   * of a given level and adjacent
-   * coarser levels) but also on
-   * the external boundary of the
-   * domain.
-   */
-  template <int dim, int spacedim>
-  void
-  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<IndexSet>  &interface_dofs,
-                                std::vector<IndexSet>  &boundary_interface_dofs);
-
-  /**
-   * Does the same as the function above,
-   * but fills only the interface_dofs.
+   * For each level in a multigrid hierarchy, produce an IndexSet
+   * that indicates which of the degrees of freedom are along
+   * interfaces of this level to cells that only exist on coarser
+   * levels.
    */
   template <int dim, int spacedim>
   void
index 9117d499167b67371a63831e2b8de1068d14bba1..4f942a4663d846a3494c6657e5f6483d36fe2f07 100644 (file)
@@ -1418,72 +1418,7 @@ namespace MGTools
       }
   }
 
-
-
-
-  template <int dim, int spacedim>
-  void
-  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<IndexSet>  &interface_dofs)
-  {
-    Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
-            ExcDimensionMismatch (interface_dofs.size(),
-                                  mg_dof_handler.get_tria().n_global_levels()));
-
-    for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
-      {
-        Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
-                ExcDimensionMismatch (interface_dofs[l].size(),
-                                      mg_dof_handler.n_dofs(l)));
-       interface_dofs[l].clear();
-      }
-
-    const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
-
-    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-    const unsigned int   dofs_per_face   = fe.dofs_per_face;
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-    std::vector<bool> cell_dofs(dofs_per_cell, false);
-
-    typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                            endc = mg_dof_handler.end();
-
-    for (; cell!=endc; ++cell)
-      {
-        if (mg_dof_handler.get_tria().locally_owned_subdomain()!=numbers::invalid_subdomain_id
-            && cell->level_subdomain_id()!=mg_dof_handler.get_tria().locally_owned_subdomain())
-          continue;
-
-        std::fill (cell_dofs.begin(), cell_dofs.end(), false);
-
-        for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
-          {
-            const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
-            if (!face->at_boundary())
-              {
-                //interior face
-                const typename DoFHandler<dim>::cell_iterator
-                neighbor = cell->neighbor(face_nr);
-
-                if (neighbor->level() < cell->level())
-                  {
-                    for (unsigned int j=0; j<dofs_per_face; ++j)
-                      cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
-
-                  }
-              }
-          }
-
-        const unsigned int level = cell->level();
-        cell->get_mg_dof_indices (local_dof_indices);
-
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          if (cell_dofs[i])
-            interface_dofs[level].add_index(local_dof_indices[i]);
-      }
-  }
-
+  
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
@@ -1579,20 +1514,14 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<IndexSet>  &interface_dofs,
-                                std::vector<IndexSet>  &boundary_interface_dofs)
+                                std::vector<IndexSet>  &interface_dofs)
   {
     Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (interface_dofs.size(),
                                   mg_dof_handler.get_tria().n_global_levels()));
-    Assert (boundary_interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
-            ExcDimensionMismatch (boundary_interface_dofs.size(),
-                                  mg_dof_handler.get_tria().n_global_levels()));
 
     std::vector<std::vector<types::global_dof_index> >
       tmp_interface_dofs(interface_dofs.size());
-    std::vector<std::vector<types::global_dof_index> >
-      tmp_boundary_interface_dofs(interface_dofs.size());
 
     const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
 
@@ -1600,10 +1529,8 @@ namespace MGTools
     const unsigned int   dofs_per_face   = fe.dofs_per_face;
 
     std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-    std::vector<types::global_dof_index> face_dof_indices (dofs_per_face);
 
     std::vector<bool> cell_dofs(dofs_per_cell, false);
-    std::vector<bool> boundary_cell_dofs(dofs_per_cell, false);
 
     typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
                                             endc = mg_dof_handler.end();
@@ -1617,7 +1544,6 @@ namespace MGTools
         bool has_coarser_neighbor = false;
 
         std::fill (cell_dofs.begin(), cell_dofs.end(), false);
-        std::fill (boundary_cell_dofs.begin(), boundary_cell_dofs.end(), false);
 
         for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
           {
@@ -1642,14 +1568,7 @@ namespace MGTools
 
         if (has_coarser_neighbor == false)
           continue;
-
-        for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
-          if (cell->at_boundary(face_nr))
-            for (unsigned int j=0; j<dofs_per_face; ++j)
-//            if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) //is this necessary?
-              boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
-
-
+       
         const unsigned int level = cell->level();
         cell->get_mg_dof_indices (local_dof_indices);
 
@@ -1657,13 +1576,9 @@ namespace MGTools
           {
             if (cell_dofs[i])
               tmp_interface_dofs[level].push_back(local_dof_indices[i]);
-
-            if (boundary_cell_dofs[i])
-              tmp_boundary_interface_dofs[level].push_back(local_dof_indices[i]);
           }
       }
 
-
     for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
       {
         interface_dofs[l].clear();
@@ -1672,13 +1587,6 @@ namespace MGTools
                                       std::unique(tmp_interface_dofs[l].begin(),
                                                   tmp_interface_dofs[l].end()));
         interface_dofs[l].compress();
-        boundary_interface_dofs[l].clear();
-        std::sort(tmp_boundary_interface_dofs[l].begin(),
-                  tmp_boundary_interface_dofs[l].end());
-        boundary_interface_dofs[l].add_indices(tmp_boundary_interface_dofs[l].begin(),
-                                               std::unique(tmp_boundary_interface_dofs[l].begin(),
-                                                           tmp_boundary_interface_dofs[l].end()));
-        boundary_interface_dofs[l].compress();
       }
 
   }
index 7f25f462f0714b93e4086c04905e78466077b06e..63f6353149114f47927255d647f54c2f4b3f6f90 100644 (file)
@@ -115,11 +115,6 @@ for (deal_II_dimension : DIMENSIONS)
        std::vector<IndexSet>&,
        const ComponentMask &);
 
-    template
-       void
-       extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
-                                     std::vector<IndexSet>  &interface_dofs,
-                                     std::vector<IndexSet>  &boundary_interface_dofs);
       template
        void
        extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
index 60e181bcf9a3d23413c1897d3ce6129d3adf542f..cd4fc73163572f2e72f003b62723287e9e07d114 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-/*
- * Authors: Guido Kanschat, University of Heidelberg, 2003
- *          Baerbel Janssen, University of Heidelberg, 2010
- *          Wolfgang Bangerth, Texas A&M University, 2010
- */
+// Multigrid for continuous finite elements using MeshWorker
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
diff --git a/tests/multigrid/step-16-bdry1.cc b/tests/multigrid/step-16-bdry1.cc
new file mode 100644 (file)
index 0000000..95359d1
--- /dev/null
@@ -0,0 +1,637 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2003 - 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.
+//
+// ---------------------------------------------------------------------
+
+// Multigrid for continuous finite elements using MeshWorker
+// Investigate a seeming inconsistency in MGConstrainedDoFs at the boundary
+
+#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/numbers.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/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/integrators/laplace.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/multigrid/mg_dof_handler.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 <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/simple.h>
+#include <deal.II/meshworker/output.h>
+#include <deal.II/meshworker/loop.h>
+
+#include <fstream>
+#include <sstream>
+
+using namespace dealii;
+using namespace LocalIntegrators;
+
+template <int dim>
+class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  LaplaceMatrix();
+  virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+  virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+  virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
+                   MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const;
+};
+
+
+template <int dim>
+LaplaceMatrix<dim>::LaplaceMatrix()
+               :
+               MeshWorker::LocalIntegrator<dim>(true, false, false)
+{}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
+{
+  AssertDimension (dinfo.n_matrices(), 1);  
+  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0));
+}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::boundary(MeshWorker::DoFInfo<dim>& /*dinfo*/,
+                                 typename MeshWorker::IntegrationInfo<dim>& /*info*/) const
+{
+//  const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
+//  Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0),
+//                       Laplace::compute_penalty(dinfo, dinfo, deg, deg));
+}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::face(
+  MeshWorker::DoFInfo<dim>& /*dinfo1*/, MeshWorker::DoFInfo<dim>& /*dinfo2*/,
+  MeshWorker::IntegrationInfo<dim>& /*info1*/, MeshWorker::IntegrationInfo<dim>& /*info2*/) const
+{
+//  const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
+//  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),
+//                  Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem (const unsigned int deg);
+  void run ();
+
+private:
+  void setup_system ();
+  void assemble_system ();
+  void assemble_multigrid (const bool &use_mw);
+  void solve ();
+  void refine_grid (const std::string& reftype);
+  void output_results (const unsigned int cycle) const;
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  MGDoFHandler<dim>    mg_dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  ConstraintMatrix     hanging_node_constraints;
+  ConstraintMatrix     constraints;
+
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+
+  const unsigned int degree;
+  LaplaceMatrix<dim> matrix_integrator;
+
+  MGLevelObject<SparsityPattern>       mg_sparsity_patterns;
+  MGLevelObject<SparseMatrix<double> > mg_matrices;
+  MGLevelObject<SparseMatrix<double> > mg_interface_in;
+  MGLevelObject<SparseMatrix<double> > mg_interface_out;
+  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),
+  matrix_integrator()
+{}
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  mg_dof_handler.distribute_dofs (fe);
+  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 (
+    static_cast<const DoFHandler<dim>&>(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;
+  MappingQ1<dim> mapping;
+  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_in.resize(0, n_levels-1);
+  mg_interface_in.clear ();
+  mg_interface_out.resize(0, n_levels-1);
+  mg_interface_out.clear ();
+  mg_matrices.resize(0, n_levels-1);
+  mg_matrices.clear ();
+  mg_sparsity_patterns.resize(0, n_levels-1);
+
+  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_in[level].reinit(mg_sparsity_patterns[level]);
+      mg_interface_out[level].reinit(mg_sparsity_patterns[level]);
+    }
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system ()
+{
+  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 MGDoFHandler<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, system_rhs);
+    }
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_multigrid (const bool& use_mw)
+{
+  if(use_mw == true)
+  {
+    mg_matrices = 0.;
+
+    MappingQ1<dim> mapping;
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags update_flags = update_values | update_gradients | update_hessians;
+    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_constrained_dofs);
+    assembler.initialize(mg_matrices);
+    assembler.initialize_interfaces(mg_interface_in, mg_interface_out);
+
+    MeshWorker::integration_loop<dim, dim> (
+        mg_dof_handler.begin(), mg_dof_handler.end(),
+        dof_info, info_box, matrix_integrator, assembler);
+
+    const unsigned int nlevels = triangulation.n_levels();
+    for (unsigned int level=0;level<nlevels;++level)
+    {
+      for(unsigned int i=0; i<mg_dof_handler.n_dofs(level); ++i)
+        if(mg_matrices[level].diag_element(i)==0)
+          mg_matrices[level].set(i,i,1.);
+    }
+  }
+  else
+  {
+    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<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 ();
+
+    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)
+    {
+      boundary_constraints[level].add_lines (interface_dofs[level]);
+      boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
+      boundary_constraints[level].close ();
+
+      boundary_interface_constraints[level]
+        .add_lines (boundary_interface_dofs[level]);
+      boundary_interface_constraints[level].close ();
+    }
+
+    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+             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()]);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          if ( !(interface_dofs[cell->level()][local_dof_indices[i]]==true &&
+                interface_dofs[cell->level()][local_dof_indices[j]]==false))
+            cell_matrix(i,j) = 0;
+
+      boundary_interface_constraints[cell->level()]
+        .distribute_local_to_global (cell_matrix,
+            local_dof_indices,
+            mg_interface_in[cell->level()]);
+    }
+  }
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve ()
+{
+  MGTransferPrebuilt<Vector<double> > mg_transfer(hanging_node_constraints, mg_constrained_dofs);
+  mg_transfer.build_matrices(mg_dof_handler);
+
+  FullMatrix<double> coarse_matrix;
+  coarse_matrix.copy_from (mg_matrices[0]);
+  MGCoarseGridHouseholder<> coarse_grid_solver;
+  coarse_grid_solver.initialize (coarse_matrix);
+
+  typedef PreconditionSOR<SparseMatrix<double> > Smoother;
+  MGSmootherRelaxation<SparseMatrix<double>, Smoother, Vector<double> >
+  mg_smoother;
+  mg_smoother.initialize(mg_matrices);
+  mg_smoother.set_steps(2);
+  mg_smoother.set_symmetric(true);
+
+  MGMatrix<> mg_matrix(&mg_matrices);
+  MGMatrix<> mg_interface_up(&mg_interface_in);
+  MGMatrix<> mg_interface_down(&mg_interface_in);
+
+  Multigrid<Vector<double> > mg(mg_dof_handler,
+                                mg_matrix,
+                                coarse_grid_solver,
+                                mg_transfer,
+                                mg_smoother,
+                                mg_smoother);
+  mg.set_edge_matrices(mg_interface_down, mg_interface_up);
+
+  PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double> > >
+  preconditioner(mg_dof_handler, mg, mg_transfer);
+
+  SolverControl solver_control (1000, 1e-12);
+  SolverCG<>    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>::refine_grid (const std::string& reftype)
+{
+  bool cell_refined = false;
+  if (reftype == "center" || !cell_refined)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator
+            cell = triangulation.begin_active();
+          cell != triangulation.end(); ++cell)
+       for (unsigned int vertex=0;
+            vertex < GeometryInfo<dim>::vertices_per_cell;
+            ++vertex)
+         {
+           {
+             const Point<dim> p = cell->vertex(vertex);
+             const Point<dim> origin = (dim == 2 ?
+                                        Point<dim>(0,0) :
+                                        Point<dim>(0,0,0));
+             const double dist = p.distance(origin);
+             if(dist<0.25/numbers::PI)
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+         }
+    }
+  if (reftype=="global" || !cell_refined)
+    triangulation.refine_global(1);
+  else
+    triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (mg_dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+
+  std::ostringstream filename;
+  filename << "solution-"
+           << cycle
+           << ".vtk";
+
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<8; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube (triangulation, -1, 1, true);
+
+//          static const HyperBallBoundary<dim> boundary;
+//          triangulation.set_boundary (0, boundary);
+
+          triangulation.refine_global (1);
+        }
+      else
+        refine_grid ("center");
+
+
+      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=0; level<triangulation.n_levels(); ++level)
+        deallog << mg_dof_handler.n_dofs(level)
+                << (level == triangulation.n_levels()-1
+                    ? ")" : ", ");
+      deallog << std::endl;
+
+      assemble_system ();
+      assemble_multigrid (true);
+
+      solve ();
+      output_results (cycle);
+    }
+}
+
+
+int main ()
+{
+  initlog();
+
+  try
+    {
+      deallog.depth_console (0);
+
+      LaplaceProblem<2> laplace_problem(1);
+      laplace_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/multigrid/step-16-bdry1.output b/tests/multigrid/step-16-bdry1.output
new file mode 100644 (file)
index 0000000..a0be674
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       4
+DEAL::Number of degrees of freedom: 9   L0: 4   L1: 9
+DEAL::   Number of degrees of freedom: 9 (by level: 4, 9)
+DEAL:cg::Starting value 1.00000
+DEAL:cg::Convergence step 1 value 0.00000
+DEAL::   1 CG iterations needed to obtain convergence.
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       16
+DEAL::Number of degrees of freedom: 25   L0: 4   L1: 9   L2: 25
+DEAL::   Number of degrees of freedom: 25 (by level: 4, 9, 25)
+DEAL:cg::Starting value 0.750000
+DEAL:cg::Convergence step 5 value 1.17215e-16
+DEAL::   5 CG iterations needed to obtain convergence.
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       28
+DEAL::Number of degrees of freedom: 41   L0: 4   L1: 9   L2: 25   L3: 25
+DEAL::   Number of degrees of freedom: 41 (by level: 4, 9, 25, 25)
+DEAL:cg::Starting value 0.628894
+DEAL:cg::Convergence step 7 value 2.44070e-14
+DEAL::   7 CG iterations needed to obtain convergence.
+DEAL::Cycle 3:
+DEAL::   Number of active cells:       40
+DEAL::Number of degrees of freedom: 57   L0: 4   L1: 9   L2: 25   L3: 25   L4: 25
+DEAL::   Number of degrees of freedom: 57 (by level: 4, 9, 25, 25, 25)
+DEAL:cg::Starting value 0.620541
+DEAL:cg::Convergence step 7 value 5.62398e-13
+DEAL::   7 CG iterations needed to obtain convergence.
+DEAL::Cycle 4:
+DEAL::   Number of active cells:       52
+DEAL::Number of degrees of freedom: 73   L0: 4   L1: 9   L2: 25   L3: 25   L4: 25   L5: 25
+DEAL::   Number of degrees of freedom: 73 (by level: 4, 9, 25, 25, 25, 25)
+DEAL:cg::Starting value 0.620015
+DEAL:cg::Convergence step 8 value 4.43441e-14
+DEAL::   8 CG iterations needed to obtain convergence.
+DEAL::Cycle 5:
+DEAL::   Number of active cells:       160
+DEAL::Number of degrees of freedom: 201   L0: 4   L1: 9   L2: 25   L3: 65   L4: 65   L5: 65   L6: 65
+DEAL::   Number of degrees of freedom: 201 (by level: 4, 9, 25, 65, 65, 65, 65)
+DEAL:cg::Starting value 0.415999
+DEAL:cg::Convergence step 9 value 4.46942e-14
+DEAL::   9 CG iterations needed to obtain convergence.
+DEAL::Cycle 6:
+DEAL::   Number of active cells:       304
+DEAL::Number of degrees of freedom: 357   L0: 4   L1: 9   L2: 25   L3: 81   L4: 81   L5: 81   L6: 81   L7: 153
+DEAL::   Number of degrees of freedom: 357 (by level: 4, 9, 25, 81, 81, 81, 81, 153)
+DEAL:cg::Starting value 0.377840
+DEAL:cg::Convergence step 8 value 1.00219e-13
+DEAL::   8 CG iterations needed to obtain convergence.
+DEAL::Cycle 7:
+DEAL::   Number of active cells:       676
+DEAL::Number of degrees of freedom: 761   L0: 4   L1: 9   L2: 25   L3: 81   L4: 81   L5: 81   L6: 121   L7: 209   L8: 465
+DEAL::   Number of degrees of freedom: 761 (by level: 4, 9, 25, 81, 81, 81, 121, 209, 465)
+DEAL:cg::Starting value 0.377713
+DEAL:cg::Convergence step 8 value 5.29691e-13
+DEAL::   8 CG iterations needed to obtain convergence.
index ab5c4b41dbccf9e7779e0ac7f71e837f4f199969..b975cb1807707738756d916f3fca38a97c0cf1d6 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-/*
- * Authors: Guido Kanschat, University of Heidelberg, 2003
- *          Baerbel Janssen, University of Heidelberg, 2010
- *          Wolfgang Bangerth, Texas A&M University, 2010
- */
+// Multigrid for continuous finite elements without MeshWorker
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>

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.