]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add MGInterfaceOperator and use it in a few tests 3605/head
authorDenis Davydov <davydden@gmail.com>
Thu, 17 Nov 2016 09:39:55 +0000 (10:39 +0100)
committerDenis Davydov <davydden@gmail.com>
Thu, 17 Nov 2016 09:39:55 +0000 (10:39 +0100)
include/deal.II/matrix_free/operators.h
tests/matrix_free/parallel_multigrid_adaptive_02.cc
tests/matrix_free/parallel_multigrid_adaptive_04.cc

index 87282f1ba8a23ae5fdd512a2551c65beb0292784..668e408a46ec0ad403f709a569f76e1ef96f1088 100644 (file)
@@ -245,6 +245,86 @@ namespace MatrixFreeOperators
 
 
 
+  /**
+   * Auxiliary class to provide interface vmult/Tvmult methods required in
+   * adaptive geometric multgrids. @p OperatorType class should be derived
+   * from MatrixFreeOperators::Base class.
+   *
+   * The adaptive multigrid realization in deal.II implements an approach
+   * called local smoothing. This means that the smoothing on the finest level
+   * only covers the local part of the mesh defined by the fixed (finest) grid
+   * level and ignores parts of the computational domain where the terminal
+   * cells are coarser than this level. As the method progresses to coarser
+   * levels, more and more of the global mesh will be covered. At some coarser
+   * level, the whole mesh will be covered. Since all level matrices in the
+   * multigrid method cover a single level in the mesh, no hanging nodes
+   * appear on the level matrices. At the interface between multigrid levels,
+   * homogeneous Dirichlet boundary conditions are set while smoothing. When
+   * the residual is transferred to the next coarser level, however, the
+   * coupling over the multigrid interface needs to be taken into account.
+   * This is done by the so-called interface (or edge) matrices that compute
+   * the part of the residual that is missed by the level matrix with
+   * homogeneous Dirichlet conditions. We refer to the @ref mg_paper
+   * "Multigrid paper by Janssen and Kanschat" for more details.
+   *
+   * For the implementation of those interface matrices, most infrastructure
+   * is already in place and provided by MatrixFreeOperators::Base through the
+   * two multiplication routines vmult_interface_down() and
+   * vmult_interface_up(). The only thing MGInterfaceOperator does is
+   * wrapping those operations and make them accessible via
+   * @p vmult() and @p Tvmult interface as expected by the multigrid routines
+   * (that were originally written for matrices, hence expecting those names).
+   * Note that the vmult_interface_down is used during the restriction phase of
+   * the multigrid V-cycle, whereas vmult_interface_up is used during the
+   * prolongation phase.
+   *
+   * @author Martin Kronbichler, 2016
+   */
+  template <typename OperatorType>
+  class MGInterfaceOperator : public Subscriptor
+  {
+  public:
+    /**
+     * Number typedef.
+     */
+    typedef typename OperatorType::value_type value_type;
+
+    /**
+     * Default constructor.
+     */
+    MGInterfaceOperator();
+
+    /**
+     * Clear the pointer to the OperatorType object.
+     */
+    void clear();
+
+    /**
+     * Initialize this class with an operator @p operator_in.
+     */
+    void initialize (const OperatorType &operator_in);
+
+    /**
+     * vmult operator, see class description for more info.
+     */
+    void vmult (LinearAlgebra::distributed::Vector<value_type> &dst,
+                const LinearAlgebra::distributed::Vector<value_type> &src) const;
+
+    /**
+     * Tvmult operator, see class description for more info.
+     */
+    void Tvmult (LinearAlgebra::distributed::Vector<value_type> &dst,
+                 const LinearAlgebra::distributed::Vector<value_type> &src) const;
+
+  private:
+    /**
+     * Const pointer to the operator class.
+     */
+    SmartPointer<const OperatorType> mf_base_operator;
+  };
+
+
+
   /**
    * This class implements the operation of the action of the inverse of a
    * mass matrix on an element for the special case of an evaluation object
@@ -965,6 +1045,62 @@ namespace MatrixFreeOperators
 
 
 
+  //------------------------- MGInterfaceOperator ------------------------------
+
+  template <typename OperatorType>
+  MGInterfaceOperator<OperatorType>::MGInterfaceOperator ()
+    :
+    Subscriptor(),
+    mf_base_operator(NULL)
+  {
+  }
+
+
+
+  template <typename OperatorType>
+  void
+  MGInterfaceOperator<OperatorType>::clear ()
+  {
+    mf_base_operator = NULL;
+  }
+
+
+
+  template <typename OperatorType>
+  void
+  MGInterfaceOperator<OperatorType>::initialize (const OperatorType &operator_in)
+  {
+    mf_base_operator = &operator_in;
+  }
+
+
+
+  template <typename OperatorType>
+  void
+  MGInterfaceOperator<OperatorType>::vmult (LinearAlgebra::distributed::Vector<typename MGInterfaceOperator<OperatorType>::value_type> &dst,
+                                            const LinearAlgebra::distributed::Vector<typename MGInterfaceOperator<OperatorType>::value_type> &src) const
+  {
+    Assert(mf_base_operator != NULL,
+           ExcNotInitialized());
+
+    mf_base_operator->vmult_interface_down(dst, src);
+  }
+
+
+
+  template <typename OperatorType>
+  void
+  MGInterfaceOperator<OperatorType>::Tvmult (LinearAlgebra::distributed::Vector<typename MGInterfaceOperator<OperatorType>::value_type> &dst,
+                                             const LinearAlgebra::distributed::Vector<typename MGInterfaceOperator<OperatorType>::value_type> &src) const
+  {
+    Assert(mf_base_operator != NULL,
+           ExcNotInitialized());
+
+    mf_base_operator->vmult_interface_up(dst, src);
+  }
+
+
+
   //-----------------------------MassOperator----------------------------------
 
   template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
index 1e99b8af183931cd22f812ef7edfa312c15d5be3..ccbf07a91d338d21afd2dcc07a3565f682d5c56a 100644 (file)
@@ -49,34 +49,6 @@ std::ofstream logfile("output");
 
 using namespace dealii::MatrixFreeOperators;
 
-
-template <typename LAPLACEOPERATOR>
-class MGInterfaceMatrix : public Subscriptor
-{
-public:
-  void initialize (const LAPLACEOPERATOR &laplace)
-  {
-    this->laplace = &laplace;
-  }
-
-  void vmult (LinearAlgebra::distributed::Vector<double> &dst,
-              const LinearAlgebra::distributed::Vector<double> &src) const
-  {
-    laplace->vmult_interface_down(dst, src);
-  }
-
-  void Tvmult (LinearAlgebra::distributed::Vector<double> &dst,
-               const LinearAlgebra::distributed::Vector<double> &src) const
-  {
-    laplace->vmult_interface_up(dst, src);
-  }
-
-private:
-  SmartPointer<const LAPLACEOPERATOR> laplace;
-};
-
-
-
 template <int dim, typename LAPLACEOPERATOR>
 class MGTransferMF : public MGTransferMatrixFree<dim, typename LAPLACEOPERATOR::value_type>
 {
@@ -239,7 +211,7 @@ void do_test (const DoFHandler<dim>  &dof)
                                     level);
       mg_matrices[level].compute_diagonal();
     }
-  MGLevelObject<MGInterfaceMatrix<LevelMatrixType> > mg_interface_matrices;
+  MGLevelObject<MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
   mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
   for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
     mg_interface_matrices[level].initialize(mg_matrices[level]);
index d23eec9f0f426debf1ec9c4f1d7d976f56b658e9..9068ad7b4bcc83563dfa548746c14f7ef05273fd 100644 (file)
@@ -49,33 +49,6 @@ std::ofstream logfile("output");
 
 using namespace dealii::MatrixFreeOperators;
 
-template <typename LAPLACEOPERATOR>
-class MGInterfaceMatrix : public Subscriptor
-{
-public:
-  void initialize (const LAPLACEOPERATOR &laplace)
-  {
-    this->laplace = &laplace;
-  }
-
-  void vmult (LinearAlgebra::distributed::Vector<double> &dst,
-              const LinearAlgebra::distributed::Vector<double> &src) const
-  {
-    laplace->vmult_interface_down(dst, src);
-  }
-
-  void Tvmult (LinearAlgebra::distributed::Vector<double> &dst,
-               const LinearAlgebra::distributed::Vector<double> &src) const
-  {
-    laplace->vmult_interface_up(dst, src);
-  }
-
-private:
-  SmartPointer<const LAPLACEOPERATOR> laplace;
-};
-
-
-
 template <int dim, typename LAPLACEOPERATOR>
 class MGTransferMF : public MGTransferMatrixFree<dim, typename LAPLACEOPERATOR::value_type>
 {
@@ -238,7 +211,7 @@ void do_test (const DoFHandler<dim>  &dof)
                                     level);
       mg_matrices[level].compute_diagonal();
     }
-  MGLevelObject<MGInterfaceMatrix<LevelMatrixType> > mg_interface_matrices;
+  MGLevelObject<MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
   mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
   for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
     mg_interface_matrices[level].initialize(mg_matrices[level]);

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.