]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Large clean-ups of multigrid things (part 1).
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Nov 1999 08:12:28 +0000 (08:12 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Nov 1999 08:12:28 +0000 (08:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@1893 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/include/numerics/multigrid.templates.h
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/lac/include/lac/mgbase.h
deal.II/lac/source/mgbase.cc

index 1a5d7134433fb8da3399d319ad44ec078a69e19e..4214e4269b9a33bfe811941a8f855a6db19fc4be 100644 (file)
 #include <vector>
 #include <base/smartpointer.h>
 
+
+
 /**
- * Abstract base class for coarse grid solvers.
- * The interface of a function call operator is defined to execute coarse
- * grid solution in a derived class.
+ * Abstract base class for coarse grid solvers.  The interface of a
+ * function call operator is defined to execute coarse grid solution
+ * in a derived class.
+ *
+ * @author Guido Kanschat, 1999
  */
 class MGCoarseGridSolver : public Subscriptor
 {
   public:
                                     /**
-                                     * Virtual destructor.
+                                     * Virtual destructor. Does
+                                     * nothing in particular, but
+                                     * needs to be declared anyhow.
                                      */
     virtual ~MGCoarseGridSolver();
     
@@ -57,24 +63,33 @@ class MGSmootherBase :  public Subscriptor
 {  
   public:
                                     /**
-                                     * Virtual destructor.
+                                     * Virtual destructor. Does
+                                     * nothing in particular, but
+                                     * needs to be declared anyhow.
                                      */
     virtual ~MGSmootherBase();
     
                                     /**
-                                     * Smooth the residual of #u# on the given
-                                     * level. This function should keep the interior
-                                     * level boundary values, so you may want
-                                     * to call #set_zero_interior_boundary#
-                                     * in the derived class #MGSmoother#
-                                     * somewhere in your derived function,
-                                     * or another function doing similar
-                                     * things.
+                                     * Smooth the residual of #u# on
+                                     * the given level. If $S$ is the
+                                     * smoothing operator, then this
+                                     * function should do the
+                                     * following operation:
+                                     * $u += S (rhs - Au)$, where #u# and
+                                     * #rhs# are the input parameters.
+                                     *
+                                     * This function should keep the
+                                     * interior level boundary
+                                     * values, so you may want to
+                                     * call #set_zero_interior_boundary#
+                                     * in the derived class
+                                     * #MGSmoother# somewhere in your
+                                     * derived function, or another
+                                     * function doing similar things.
                                      */
     virtual void smooth (const unsigned int    level,
                         Vector<double>       &u,
                         const Vector<double> &rhs) const = 0;
-
 };
 
 
@@ -93,7 +108,8 @@ class MGSmootherIdentity : public MGSmootherBase
 {
   public:
                                     /**
-                                     * Implementation of the interface in #MGSmootherBase#.
+                                     * Implementation of the
+                                     * interface in #MGSmootherBase#.
                                      * This function does nothing.
                                      */
     virtual void smooth (const unsigned int   level,
@@ -103,6 +119,7 @@ class MGSmootherIdentity : public MGSmootherBase
 
 
 
+
 /**
  * Base class used to declare the operations needed by a concrete class
  * implementing prolongation and restriction of vectors in the multigrid
@@ -132,7 +149,7 @@ class MGTransferBase : public Subscriptor
                                      * are degrees of freedom on the finer
                                      * level.
                                      */
-    virtual void prolongate (const unsigned int   to_level,
+    virtual void prolongate (const unsigned int    to_level,
                             Vector<double>       &dst,
                             const Vector<double> &src) const = 0;
 
@@ -149,7 +166,7 @@ class MGTransferBase : public Subscriptor
                                      * are degrees of freedom on the coarser
                                      * level.
                                      */
-    virtual void restrict (const unsigned int   from_level,
+    virtual void restrict (const unsigned int    from_level,
                           Vector<double>       &dst,
                           const Vector<double> &src) const = 0;
 };
@@ -157,37 +174,61 @@ class MGTransferBase : public Subscriptor
 
 
 /**
- * An array with a vector for each level.
- * The purpose of this class is mostly to provide range checking that is missing in
- * those ultra-fashionable STL vectors.
+ * An array with a vector for each level.  The purpose of this class
+ * is mostly to allow access by level number, even if the lower levels
+ * are not used and therefore have no vector at all; this is done by
+ * simply shifting the given index by the minimum level we have
+ * stored.
+ *
  * @author Guido Kanschat, 1999
  */
-template<class VECTOR>
-class MGVector : public Subscriptor,
-                public vector<VECTOR>
+template<class VECTOR = Vector<double> >
+class MGVector : public Subscriptor
 {
   public:
                                     /**
-                                     * Constructor allowing to initialize the number of levels.
+                                     * Constructor allowing to
+                                     * initialize the number of
+                                     * levels.
                                      */
-    MGVector(unsigned int minlevel, unsigned int maxlevel);
+    MGVector (const unsigned int minlevel,
+             const unsigned int maxlevel);
+    
                                     /**
-                                     * Safe access operator.
+                                     * Access vector on level #level#.
                                      */
-    VECTOR& operator[](unsigned int);
+    VECTOR & operator[] (const unsigned int level);
+    
                                     /**
-                                     * Safe access operator.
+                                     * Access vector on level
+                                     * #level#. Constant version.
                                      */
-    const VECTOR& operator[](unsigned int) const;
+    const VECTOR & operator[] (const unsigned int level) const;
+    
                                     /**
-                                     * Delete all entries.
+                                     * Call #clear# on all vectors
+                                     * stored by this object. Note
+                                     * that if #VECTOR==Vector<T>#,
+                                     * #clear# will set all entries
+                                     * to zero, while if
+                                     * #VECTOR==std::vector<T>#,
+                                     * #clear# deletes the elements
+                                     * of the vectors. This class
+                                     * might therefore not be useful
+                                     * for STL vectors.
                                      */
     void clear();
+    
   private:
                                     /**
                                      * Level of first component.
                                      */
-    unsigned int minlevel;
+    const unsigned int minlevel;
+
+                                    /**
+                                     * Array of the vectors to be held.
+                                     */
+    vector<VECTOR> vectors;
 };
 
 
@@ -234,6 +275,10 @@ class MGMatrix : public Subscriptor,
  * This is a little wrapper, transforming a triplet of iterative
  * solver, matrix and preconditioner into a coarse grid solver.
  *
+ * The type of the matrix (i.e. the template parameter #MATRIX#)
+ * should be derived from #Subscriptor# to allow for the use of a
+ * smart pointer to it.
+ *
  * @author Guido Kanschat, 1999
  */
 template<class SOLVER, class MATRIX, class PRECOND>
@@ -246,7 +291,9 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
                                      * preconditioning method for later
                                      * use.
                                      */
-    MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&);
+    MGCoarseGridLACIteration (SOLVER        &,
+                             const MATRIX  &,
+                             const PRECOND &);
     
                                     /**
                                      * Implementation of the abstract
@@ -267,7 +314,7 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
                                     /**
                                      * Reference to the matrix.
                                      */
-    const MATRIX& matrix;
+    const SmartPointer<const MATRIX> matrix;
     
                                     /**
                                      * Reference to the preconditioner.
@@ -277,16 +324,16 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
 
 
 
+
 /**
- * Basic class for multigrid preconditioning.
+ * Basic class for preconditioning by multigrid.
  *
  * The functionality of the multigrid method is restricted to defect
- * correction. It is <B>not</B> iterative and the start solution is
- * always zero. Since by this <I>u<SUP>E</SUP><SUB>l</SUB></I> and
- * <I>u<SUP>A</SUP><SUB>l</SUB></I> (see report on multigrid) are
- * always zero, restriction is simplified a lot and maybe even the
- * seam condition on grids is oblivious. Still, I am not sure that
- * these restrictions on the class might cause numerical
+ * correction. It is {\bf not} iterative and the start solution is
+ * always zero. Since by this $u^E_l$ and $u^A_l$ (see report on
+ * multigrid) are always zero, restriction is simplified a lot and
+ * maybe even the seam condition on grids is oblivious. Still, I am
+ * not sure that these restrictions on the class might cause numerical
  * inefficiencies.
  *
  * The function #precondition# is the actual multigrid method and
@@ -295,65 +342,109 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
  * preconditioned defect in #dst#.
  *
  * @author Guido Kanschat, 1999
- * */
+ */
 class MGBase : public Subscriptor
 {
   private:
+                                    /**
+                                     * Copy constructor. Made private
+                                     * to prevent use.
+                                     */
     MGBase(const MGBase&);
+
+                                    /**
+                                     * Copy operator. Made private
+                                     * to prevent use.
+                                     */
     const MGBase& operator=(const MGBase&);
 
-  protected:
+  public:
                                     /**
-                                     * Highest level of cells.
+                                     * Constructor, subject to change.
                                      */
-    unsigned int maxlevel;
+    MGBase (const MGTransferBase &transfer,
+           const unsigned int    minlevel,
+           const unsigned int    maxlevel);
+    
+                                    /**
+                                     * Virtual destructor.
+                                     */
+    virtual ~MGBase();
+    
+                                    /**
+                                     * Execute one step of the
+                                     * v-cycle algorithm.  This
+                                     * function assumes, that the
+                                     * vector #d# is properly filled
+                                     * with the residual in the outer
+                                     * defect correction
+                                     * scheme. After execution of
+                                     * #vcycle()#, the result is in
+                                     * the vector #s#. We propose to
+                                     * write functions #copy_from_mg#
+                                     * and #copy_to_mg# in derived
+                                     * classes of #MGBase# to access
+                                     * #d# and #s#.
+                                     *
+                                     * The actual work for this
+                                     * function is done in
+                                     * #level_mgstep#.
+                                     */
+    void vcycle(const MGSmootherBase     &pre_smooth,
+               const MGSmootherBase     &post_smooth,
+               const MGCoarseGridSolver &cgs);
 
                                     /**
-                                     * Level for coarse grid solution.
+                                     * Exception.
                                      */
-    unsigned int minlevel;
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
+    
+  protected:
                                     /**
-                                     * Auxiliary vector, defect.
+                                     * Highest level of cells.
                                      */
-    MGVector<Vector<double> > d;
+    unsigned int maxlevel;
 
                                     /**
-                                     * Auxiliary vector, solution.
+                                     * Level for coarse grid solution.
                                      */
-    MGVector<Vector<double> > s;
+    unsigned int minlevel;
+    
                                     /**
-                                     * Prolongation and restriction object.
+                                     * Auxiliary vector.
                                      */
-    SmartPointer<const MGTransferBase> transfer;
+    MGVector<Vector<double> > defect;
 
-  private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<double> t;
+    MGVector<Vector<double> > solution;
     
                                     /**
-                                     * Exception.
+                                     * Prolongation and restriction object.
                                      */
-    DeclException2(ExcSwitchedLevels, int, int,
-                  << "minlevel and maxlevel switched, should be: "
-                  << arg1 << "<=" << arg2);
-    
-  protected:
+    SmartPointer<const MGTransferBase> transfer;
 
                                     /**
                                      * The actual v-cycle multigrid method.
-                                     * This function is called on the
-                                     * highest level and recursively
-                                     * invokes itself down to the
-                                     * coarsest. There, it calls
-                                     * #coarse_grid_solver# and
-                                     * proceeds back up.
-                                     */
-    void level_mgstep(unsigned int level,
-                     const MGSmootherBase& pre_smooth,
-                     const MGSmootherBase& post_smooth,
-                     const MGCoarseGridSolver& cgs);  
+                                     * #level# is the level the
+                                     * function should work on. It
+                                     * will usually be called for the
+                                     * highest level from outside,
+                                     * but will then call itself
+                                     * recursively for #level-1#,
+                                     * unless we are on #minlevel#
+                                     * where instead of recursing
+                                     * deeper, the coarse grid solver
+                                     * is used to solve the matrix of
+                                     * this level.
+                                     */
+    void level_mgstep (const unsigned int        level,
+                      const MGSmootherBase     &pre_smooth,
+                      const MGSmootherBase     &post_smooth,
+                      const MGCoarseGridSolver &cgs);  
 
                                     /**
                                      * Apply negative #vmult# on all
@@ -361,39 +452,16 @@ class MGBase : public Subscriptor
                                      * This is implemented in a
                                      * derived class.
                                      */
-    virtual void level_vmult(unsigned int level,
-                            Vector<double>& dst,
-                            const Vector<double>& src,
-                            const Vector<double>& rhs) = 0;  
-
+    virtual void level_vmult (const unsigned int    level,
+                             Vector<double>       &dst,
+                             const Vector<double> &src,
+                             const Vector<double> &rhs) = 0;
   
-  public:
-                                    /**
-                                     * Constructor, subject to change.
-                                     */
-    MGBase(const MGTransferBase& transfer,
-                 unsigned int minlevel, unsigned int maxlevel);
-                                    /**
-                                     * Virtual destructor.
-                                     */
-    virtual ~MGBase();
+  private:
                                     /**
-                                     * Execute v-cycle algorithm.
-                                     * This function assumes, that
-                                     * the vector #d# is properly
-                                     * filled with the residual in
-                                     * the outer defect correction
-                                     * scheme. After execution of
-                                     * #vcycle()#, the result is in
-                                     * the vector #s#. We propose to
-                                     * write functions #copy_from_mg#
-                                     * and #copy_to_mg# in derived
-                                     * classes of #MGBase# to access
-                                     * #d# and #s#.
+                                     * Auxiliary vector.
                                      */
-    void vcycle(const MGSmootherBase& pre_smooth,
-               const MGSmootherBase& post_smooth,
-               const MGCoarseGridSolver& cgs);
+    Vector<double> t;    
 };
 
 
@@ -407,9 +475,10 @@ class MGBase : public Subscriptor
  * Furthermore, it needs functions #void copy_to_mg(const VECTOR&)#
  * to store #src# in the right hand side of the multi-level method and
  * #void copy_from_mg(VECTOR&)# to store the result of the v-cycle in #dst#.
+ *
  * @author Guido Kanschat, 1999
  */
-template<class MG, class VECTOR>
+template<class MG, class VECTOR = Vector<double> >
 class PreconditionMG
 {
   public:
@@ -419,29 +488,39 @@ class PreconditionMG
                                      * pre-smoother, post-smoother and
                                      * coarse grid solver.
                                      */
-    PreconditionMG(MG& mg,
-                  const MGSmootherBase& pre,
-                  const MGSmootherBase& post,
-                  const MGCoarseGridSolver& coarse);
+    PreconditionMG(MG                       &mg,
+                  const MGSmootherBase     &pre,
+                  const MGSmootherBase     &post,
+                  const MGCoarseGridSolver &coarse);
+    
                                     /**
-                                     * Preconditioning operator.
-                                     * This is the operator used by LAC
-                                     * iterative solvers.
+                                     * Preconditioning operator,
+                                     * calling the #vcycle# function
+                                     * of the #MG# object passed to
+                                     * the constructor.
+                                     *
+                                     * This is the operator used by
+                                     * LAC iterative solvers.
                                      */
-    void operator() (VECTOR& dst, const VECTOR& src) const;
+    void operator() (VECTOR       &dst,
+                    const VECTOR &src) const;
+    
   private:
                                     /**
-                                     * The mulrigrid object.
+                                     * The multigrid object.
                                      */
     SmartPointer<MG> multigrid;
+    
                                     /**
                                      * The pre-smoothing object.
                                      */
     SmartPointer<const MGSmootherBase> pre;
+
                                     /**
                                      * The post-smoothing object.
                                      */
     SmartPointer<const MGSmootherBase> post;
+    
                                     /**
                                      * The coarse grid solver.
                                      */
@@ -451,58 +530,44 @@ class PreconditionMG
 
 
 
-
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
-MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
-               :
-               solver(s), matrix(m), precondition(p)
-{}
-
-template<class SOLVER, class MATRIX, class PRECOND>
-void
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<double>& dst, const Vector<double>& src) const
-{
-  solver.solve(matrix, dst, src, precondition);
-}
-
-
-
-
 /* ------------------------------------------------------------------- */
 
 template<class VECTOR>
-MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+MGVector<VECTOR>::MGVector(const unsigned int min,
+                          const unsigned int max)
                :
-               vector<VECTOR>(max-min+1),
-               minlevel(min)
-{}
+               minlevel(min),
+               vectors(max-min+1)
+{};
+
 
 
 template<class VECTOR>
-VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i)
+VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i)
 {
-  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
-  return vector<VECTOR>::operator[](i-minlevel);
+  Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+        ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+  return vectors[i-minlevel];
 }
 
 
 template<class VECTOR>
-const VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i) const
+const VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i) const
 {
-  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
-  return vector<VECTOR>::operator[](i-minlevel);
+  Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+        ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+  return vectors[i-minlevel];
 }
 
+
 template<class VECTOR>
 void
 MGVector<VECTOR>::clear()
 {
   typename vector<VECTOR>::iterator v;
-  for (v = begin(); v != end(); ++v)
+  for (v = vectors.begin(); v != vectors.end(); ++v)
     v->clear();
   
 }
@@ -525,7 +590,8 @@ MGMatrix<MATRIX>::operator[](unsigned int i)
 {
   Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
   return vector<MATRIX>::operator[](i-minlevel);
-}
+};
+
 
 
 template<class MATRIX>
@@ -534,18 +600,49 @@ MGMatrix<MATRIX>::operator[](unsigned int i) const
 {
   Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
   return vector<MATRIX>::operator[](i-minlevel);
+};
+
+
+
+/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+MGCoarseGridLACIteration(SOLVER        &s,
+                        const MATRIX  &m,
+                        const PRECOND &p)
+               :
+               solver(s),
+               matrix(&m),
+               precondition(p)
+{};
+
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+void
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+operator() (const unsigned int    /* level */,
+           Vector<double>       &dst,
+           const Vector<double> &src) const
+{
+  solver.solve(*matrix, dst, src, precondition);
 }
 
 
 
+
+
+
 /* ------------------------------------------------------------------- */
 
 
 template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MGmg,
-                                          const MGSmootherBasepre,
-                                          const MGSmootherBasepost,
-                                          const MGCoarseGridSolvercoarse)
+PreconditionMG<MG, VECTOR>::PreconditionMG(MG                       &mg,
+                                          const MGSmootherBase     &pre,
+                                          const MGSmootherBase     &post,
+                                          const MGCoarseGridSolver &coarse)
                :
                multigrid(&mg),
                pre(&pre),
@@ -557,8 +654,8 @@ PreconditionMG<MG, VECTOR>::PreconditionMG(MG& mg,
 
 template<class MG, class VECTOR>
 void
-PreconditionMG<MG, VECTOR>::operator() (VECTOR& dst,
-                                       const VECTOR& src) const
+PreconditionMG<MG,VECTOR>::operator() (VECTOR       &dst,
+                                      const VECTOR &src) const
 {
   multigrid->copy_to_mg(src);
   multigrid->vcycle(*pre, *post, *coarse);
index 7caabae233640f59cb63d58acc0a485410e23e36..5d01b9b16ae90a5e37e2ad227c2fb78dea2aa880 100644 (file)
@@ -1,12 +1,11 @@
 // $Id$
 
 #include <numerics/multigrid.h>
+#include <algorithm>
 #include <fstream>
 
-inline int minimum(int a, int b)
-{
-  return ((a<b) ? a : b);
-}
+
+
 
 template <int dim>
 MG<dim>::MG(const MGDoFHandler<dim>& dofs,
@@ -16,17 +15,19 @@ MG<dim>::MG(const MGDoFHandler<dim>& dofs,
            unsigned int minl, unsigned int maxl)
                :
                MGBase(transfer, minl,
-                      minimum(dofs.get_tria().n_levels()-1, maxl)),
+                      min(dofs.get_tria().n_levels()-1, maxl)),
                dofs(&dofs),
                matrices(&matrices),
                hanging_nodes(hanging_nodes)
 {
   for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
-      s[level].reinit(matrices[level].m());
-      d[level].reinit(matrices[level].m());
-    }
-}
+      solution[level].reinit(matrices[level].m());
+      defect[level].reinit(matrices[level].m());
+    };
+};
+
+
 
 template <int dim>
 template <typename number>
@@ -38,7 +39,7 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
 
                                   // set the elements of the vectors
                                   // on all levels to zero
-  d.clear();
+  defect.clear();
 //  hanging_nodes.condense(src);
   
   vector<int> index(fe_dofs);
@@ -65,7 +66,7 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
          dc->get_dof_indices(index);
          c->get_mg_dof_indices(mgindex);
          for (unsigned int i=0;i<fe_dofs;++i)
-           d[c->level()](mgindex[i]) = src(index[i]);
+           defect[c->level()](mgindex[i]) = src(index[i]);
 
                                           // Delete values on refinement edge
          for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
@@ -75,14 +76,14 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
                {
                  face->get_mg_dof_indices(mgfaceindex);
                  for (unsigned int i=0;i<face_dofs;++i)
-                   d[c->level()](mgfaceindex[i]) = 0.;
+                   defect[c->level()](mgfaceindex[i]) = 0.;
                }
            }
        }
 //      if (level > (int) minlevel)
 //        transfer->restrict(level, d[level-1], d[level]);
       if (level < (int) maxlevel)
-       transfer->restrict(level+1, d[level], d[level+1]);
+       transfer->restrict(level+1, defect[level], defect[level+1]);
     }
   
 //   for (unsigned int i=minlevel; i<= maxlevel; ++i)
@@ -94,7 +95,8 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
 //     }
 }
 
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
 void
 MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
@@ -112,13 +114,13 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //     }
 
   DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
-  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
-                                           ; c != dofs->end() ; ++c, ++dc)
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active();
+       c != dofs->end(); ++c, ++dc)
     {
       dc->get_dof_indices(index);
       c->get_mg_dof_indices(mgindex);
       for (unsigned int i=0;i<fe_dofs;++i)
-       dst(index[i]) = s[c->level()](mgindex[i]);
+       dst(index[i]) = solution[c->level()](mgindex[i]);
     } 
   hanging_nodes.set_zero(dst);
 //   {
@@ -130,6 +132,8 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //   }
 }
 
+
+
 template <int dim>
 void
 MG<dim>::level_vmult(unsigned int l,
index 7caabae233640f59cb63d58acc0a485410e23e36..5d01b9b16ae90a5e37e2ad227c2fb78dea2aa880 100644 (file)
@@ -1,12 +1,11 @@
 // $Id$
 
 #include <numerics/multigrid.h>
+#include <algorithm>
 #include <fstream>
 
-inline int minimum(int a, int b)
-{
-  return ((a<b) ? a : b);
-}
+
+
 
 template <int dim>
 MG<dim>::MG(const MGDoFHandler<dim>& dofs,
@@ -16,17 +15,19 @@ MG<dim>::MG(const MGDoFHandler<dim>& dofs,
            unsigned int minl, unsigned int maxl)
                :
                MGBase(transfer, minl,
-                      minimum(dofs.get_tria().n_levels()-1, maxl)),
+                      min(dofs.get_tria().n_levels()-1, maxl)),
                dofs(&dofs),
                matrices(&matrices),
                hanging_nodes(hanging_nodes)
 {
   for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
-      s[level].reinit(matrices[level].m());
-      d[level].reinit(matrices[level].m());
-    }
-}
+      solution[level].reinit(matrices[level].m());
+      defect[level].reinit(matrices[level].m());
+    };
+};
+
+
 
 template <int dim>
 template <typename number>
@@ -38,7 +39,7 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
 
                                   // set the elements of the vectors
                                   // on all levels to zero
-  d.clear();
+  defect.clear();
 //  hanging_nodes.condense(src);
   
   vector<int> index(fe_dofs);
@@ -65,7 +66,7 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
          dc->get_dof_indices(index);
          c->get_mg_dof_indices(mgindex);
          for (unsigned int i=0;i<fe_dofs;++i)
-           d[c->level()](mgindex[i]) = src(index[i]);
+           defect[c->level()](mgindex[i]) = src(index[i]);
 
                                           // Delete values on refinement edge
          for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
@@ -75,14 +76,14 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
                {
                  face->get_mg_dof_indices(mgfaceindex);
                  for (unsigned int i=0;i<face_dofs;++i)
-                   d[c->level()](mgfaceindex[i]) = 0.;
+                   defect[c->level()](mgfaceindex[i]) = 0.;
                }
            }
        }
 //      if (level > (int) minlevel)
 //        transfer->restrict(level, d[level-1], d[level]);
       if (level < (int) maxlevel)
-       transfer->restrict(level+1, d[level], d[level+1]);
+       transfer->restrict(level+1, defect[level], defect[level+1]);
     }
   
 //   for (unsigned int i=minlevel; i<= maxlevel; ++i)
@@ -94,7 +95,8 @@ MG<dim>::copy_to_mg(const Vector<number>& src)
 //     }
 }
 
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
 void
 MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
@@ -112,13 +114,13 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //     }
 
   DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
-  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
-                                           ; c != dofs->end() ; ++c, ++dc)
+  for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active();
+       c != dofs->end(); ++c, ++dc)
     {
       dc->get_dof_indices(index);
       c->get_mg_dof_indices(mgindex);
       for (unsigned int i=0;i<fe_dofs;++i)
-       dst(index[i]) = s[c->level()](mgindex[i]);
+       dst(index[i]) = solution[c->level()](mgindex[i]);
     } 
   hanging_nodes.set_zero(dst);
 //   {
@@ -130,6 +132,8 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //   }
 }
 
+
+
 template <int dim>
 void
 MG<dim>::level_vmult(unsigned int l,
index 9df47ff5949ba4b6e4ce8af46e129b7cc60caba5..3662347101cb84a7124a0ff1dc5036512899aea4 100644 (file)
 #include <iostream>
 #include <cmath>
 
+
+
+//TODO: this function is only for debugging purposes and should be removed sometimes
 template<class VECTOR>
+static
 void print_vector(ostream& s, const VECTOR& v)
 {
   const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
   unsigned int k=0;
-  
+
+                                  // write the vector entries in a
+                                  // kind of square
   for (unsigned int i=0;i<n;++i)
     {
       for (unsigned int j=0;j<n;++j)
-//     s << ((float)i)/n << '\t' << ((float)j)/n << '\t' << v(k++) << endl;
        s << '\t' << v(k++);
       s << endl;
     }
   s << endl;
 }
 
+
+
 MGBase::~MGBase () 
 {};
 
 
+//////////////////////////////////////////////////////////////////////
+
+
+MGSmootherBase::~MGSmootherBase()
+{};
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+                           Vector<double>       &,
+                           const Vector<double> &) const
+{};
+
 
-MGBase::MGBase(const MGTransferBase& transfer,
-                            unsigned minlevel, unsigned maxlevel)
+
+//////////////////////////////////////////////////////////////////////
+
+
+
+MGBase::MGBase(const MGTransferBase &transfer,
+              const unsigned        minlevel,
+              const unsigned        maxlevel)
                :
-               maxlevel(maxlevel), minlevel(minlevel),
-               d(minlevel,maxlevel),
-               s(minlevel,maxlevel),
+               maxlevel(maxlevel),
+               minlevel(minlevel),
+               defect(minlevel,maxlevel),
+               solution(minlevel,maxlevel),
                transfer(&transfer)
 {
-  Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
-}
+  Assert(minlevel <= maxlevel,
+        ExcSwitchedLevels(minlevel, maxlevel));
+};
+
+
 
 void
-MGBase::vcycle(const MGSmootherBasepre_smooth,
-              const MGSmootherBasepost_smooth,
-              const MGCoarseGridSolvercoarse_grid_solver)
+MGBase::vcycle(const MGSmootherBase     &pre_smooth,
+              const MGSmootherBase     &post_smooth,
+              const MGCoarseGridSolver &coarse_grid_solver)
 {
   level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-}
+};
+
 
 
 void
-MGBase::level_mgstep(unsigned level,
-                    const MGSmootherBasepre_smooth,
-                    const MGSmootherBasepost_smooth,
-                    const MGCoarseGridSolvercoarse_grid_solver)
+MGBase::level_mgstep(const unsigned int        level,
+                    const MGSmootherBase     &pre_smooth,
+                    const MGSmootherBase     &post_smooth,
+                    const MGCoarseGridSolver &coarse_grid_solver)
 {
-  s[level] = 0.;
+  solution[level] = 0.;
   
   if (level == minlevel)
     {
-      coarse_grid_solver(level, s[level], d[level]);
+      coarse_grid_solver(level, solution[level], defect[level]);
       return;
     }
 
                           // smoothing of the residual by modifying s
-  pre_smooth.smooth(level, s[level], d[level]);
+  pre_smooth.smooth(level, solution[level], defect[level]);
                                   // t = d-As
-  t.reinit(s[level].size());
-  level_vmult(level, t, s[level], d[level]);
+  t.reinit(solution[level].size());
+  level_vmult(level, t, solution[level], defect[level]);
 
                                   // make t rhs of lower level
-  transfer->restrict(level, d[level-1], t);
+  transfer->restrict(level, defect[level-1], t);
   
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
                                   // do coarse grid correction
-  t.reinit(s[level].size());
-  transfer->prolongate(level, t, s[level-1]);
-  s[level].add(t);
+  t.reinit(solution[level].size());
+  transfer->prolongate(level, t, solution[level-1]);
+  solution[level].add(t);
   
                                   // smoothing (modify s again)
-  post_smooth.smooth(level, s[level], d[level]);
+  post_smooth.smooth(level, solution[level], defect[level]);
 }
 
 
 //////////////////////////////////////////////////////////////////////
 
 MGCoarseGridSolver::~MGCoarseGridSolver()
-{}
+{};
+
 
 
 //////////////////////////////////////////////////////////////////////
 
 MGTransferBase::~MGTransferBase()
-{}
-
-MGSmootherBase::~MGSmootherBase()
-{}
-
+{};
 
-//////////////////////////////////////////////////////////////////////
 
 
-void
-MGSmootherIdentity::smooth (const unsigned int,
-                           Vector<double>       &,
-                           const Vector<double> &) const
-{}
 
 
index d1fd1fe0c3a9c9ff90f06c8aa950fe14a3800155..9456174df4d5552d83484f370441ffa0a1ea6ea4 100644 (file)
 #include <vector>
 #include <base/smartpointer.h>
 
+
+
 /**
- * Abstract base class for coarse grid solvers.
- * The interface of a function call operator is defined to execute coarse
- * grid solution in a derived class.
+ * Abstract base class for coarse grid solvers.  The interface of a
+ * function call operator is defined to execute coarse grid solution
+ * in a derived class.
+ *
+ * @author Guido Kanschat, 1999
  */
 class MGCoarseGridSolver : public Subscriptor
 {
   public:
                                     /**
-                                     * Virtual destructor.
+                                     * Virtual destructor. Does
+                                     * nothing in particular, but
+                                     * needs to be declared anyhow.
                                      */
     virtual ~MGCoarseGridSolver();
     
@@ -57,24 +63,33 @@ class MGSmootherBase :  public Subscriptor
 {  
   public:
                                     /**
-                                     * Virtual destructor.
+                                     * Virtual destructor. Does
+                                     * nothing in particular, but
+                                     * needs to be declared anyhow.
                                      */
     virtual ~MGSmootherBase();
     
                                     /**
-                                     * Smooth the residual of #u# on the given
-                                     * level. This function should keep the interior
-                                     * level boundary values, so you may want
-                                     * to call #set_zero_interior_boundary#
-                                     * in the derived class #MGSmoother#
-                                     * somewhere in your derived function,
-                                     * or another function doing similar
-                                     * things.
+                                     * Smooth the residual of #u# on
+                                     * the given level. If $S$ is the
+                                     * smoothing operator, then this
+                                     * function should do the
+                                     * following operation:
+                                     * $u += S (rhs - Au)$, where #u# and
+                                     * #rhs# are the input parameters.
+                                     *
+                                     * This function should keep the
+                                     * interior level boundary
+                                     * values, so you may want to
+                                     * call #set_zero_interior_boundary#
+                                     * in the derived class
+                                     * #MGSmoother# somewhere in your
+                                     * derived function, or another
+                                     * function doing similar things.
                                      */
     virtual void smooth (const unsigned int    level,
                         Vector<double>       &u,
                         const Vector<double> &rhs) const = 0;
-
 };
 
 
@@ -93,7 +108,8 @@ class MGSmootherIdentity : public MGSmootherBase
 {
   public:
                                     /**
-                                     * Implementation of the interface in #MGSmootherBase#.
+                                     * Implementation of the
+                                     * interface in #MGSmootherBase#.
                                      * This function does nothing.
                                      */
     virtual void smooth (const unsigned int   level,
@@ -103,6 +119,7 @@ class MGSmootherIdentity : public MGSmootherBase
 
 
 
+
 /**
  * Base class used to declare the operations needed by a concrete class
  * implementing prolongation and restriction of vectors in the multigrid
@@ -132,7 +149,7 @@ class MGTransferBase : public Subscriptor
                                      * are degrees of freedom on the finer
                                      * level.
                                      */
-    virtual void prolongate (const unsigned int   to_level,
+    virtual void prolongate (const unsigned int    to_level,
                             Vector<double>       &dst,
                             const Vector<double> &src) const = 0;
 
@@ -149,7 +166,7 @@ class MGTransferBase : public Subscriptor
                                      * are degrees of freedom on the coarser
                                      * level.
                                      */
-    virtual void restrict (const unsigned int   from_level,
+    virtual void restrict (const unsigned int    from_level,
                           Vector<double>       &dst,
                           const Vector<double> &src) const = 0;
 };
@@ -157,37 +174,61 @@ class MGTransferBase : public Subscriptor
 
 
 /**
- * An array with a vector for each level.
- * The purpose of this class is mostly to provide range checking that is missing in
- * those ultra-fashionable STL vectors.
+ * An array with a vector for each level.  The purpose of this class
+ * is mostly to allow access by level number, even if the lower levels
+ * are not used and therefore have no vector at all; this is done by
+ * simply shifting the given index by the minimum level we have
+ * stored.
+ *
  * @author Guido Kanschat, 1999
  */
-template<class VECTOR>
-class MGVector : public Subscriptor,
-                public vector<VECTOR>
+template<class VECTOR = Vector<double> >
+class MGVector : public Subscriptor
 {
   public:
                                     /**
-                                     * Constructor allowing to initialize the number of levels.
+                                     * Constructor allowing to
+                                     * initialize the number of
+                                     * levels.
                                      */
-    MGVector(unsigned int minlevel, unsigned int maxlevel);
+    MGVector (const unsigned int minlevel,
+             const unsigned int maxlevel);
+    
                                     /**
-                                     * Safe access operator.
+                                     * Access vector on level #level#.
                                      */
-    VECTOR& operator[](unsigned int);
+    VECTOR & operator[] (const unsigned int level);
+    
                                     /**
-                                     * Safe access operator.
+                                     * Access vector on level
+                                     * #level#. Constant version.
                                      */
-    const VECTOR& operator[](unsigned int) const;
+    const VECTOR & operator[] (const unsigned int level) const;
+    
                                     /**
-                                     * Delete all entries.
+                                     * Call #clear# on all vectors
+                                     * stored by this object. Note
+                                     * that if #VECTOR==Vector<T>#,
+                                     * #clear# will set all entries
+                                     * to zero, while if
+                                     * #VECTOR==std::vector<T>#,
+                                     * #clear# deletes the elements
+                                     * of the vectors. This class
+                                     * might therefore not be useful
+                                     * for STL vectors.
                                      */
     void clear();
+    
   private:
                                     /**
                                      * Level of first component.
                                      */
-    unsigned int minlevel;
+    const unsigned int minlevel;
+
+                                    /**
+                                     * Array of the vectors to be held.
+                                     */
+    vector<VECTOR> vectors;
 };
 
 
@@ -234,6 +275,10 @@ class MGMatrix : public Subscriptor,
  * This is a little wrapper, transforming a triplet of iterative
  * solver, matrix and preconditioner into a coarse grid solver.
  *
+ * The type of the matrix (i.e. the template parameter #MATRIX#)
+ * should be derived from #Subscriptor# to allow for the use of a
+ * smart pointer to it.
+ *
  * @author Guido Kanschat, 1999
  */
 template<class SOLVER, class MATRIX, class PRECOND>
@@ -246,7 +291,9 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
                                      * preconditioning method for later
                                      * use.
                                      */
-    MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&);
+    MGCoarseGridLACIteration (SOLVER        &,
+                             const MATRIX  &,
+                             const PRECOND &);
     
                                     /**
                                      * Implementation of the abstract
@@ -267,7 +314,7 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
                                     /**
                                      * Reference to the matrix.
                                      */
-    const MATRIX& matrix;
+    const SmartPointer<const MATRIX> matrix;
     
                                     /**
                                      * Reference to the preconditioner.
@@ -277,16 +324,16 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
 
 
 
+
 /**
- * Basic class for multigrid preconditioning.
+ * Basic class for preconditioning by multigrid.
  *
  * The functionality of the multigrid method is restricted to defect
- * correction. It is <B>not</B> iterative and the start solution is
- * always zero. Since by this <I>u<SUP>E</SUP><SUB>l</SUB></I> and
- * <I>u<SUP>A</SUP><SUB>l</SUB></I> (see report on multigrid) are
- * always zero, restriction is simplified a lot and maybe even the
- * seam condition on grids is oblivious. Still, I am not sure that
- * these restrictions on the class might cause numerical
+ * correction. It is {\bf not} iterative and the start solution is
+ * always zero. Since by this $u^E_l$ and $u^A_l$ (see report on
+ * multigrid) are always zero, restriction is simplified a lot and
+ * maybe even the seam condition on grids is oblivious. Still, I am
+ * not sure that these restrictions on the class might cause numerical
  * inefficiencies.
  *
  * The function #precondition# is the actual multigrid method and
@@ -295,65 +342,109 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
  * preconditioned defect in #dst#.
  *
  * @author Guido Kanschat, 1999
- * */
+ */
 class MGBase : public Subscriptor
 {
   private:
+                                    /**
+                                     * Copy constructor. Made private
+                                     * to prevent use.
+                                     */
     MGBase(const MGBase&);
+
+                                    /**
+                                     * Copy operator. Made private
+                                     * to prevent use.
+                                     */
     const MGBase& operator=(const MGBase&);
 
-  protected:
+  public:
                                     /**
-                                     * Highest level of cells.
+                                     * Constructor, subject to change.
                                      */
-    unsigned int maxlevel;
+    MGBase (const MGTransferBase &transfer,
+           const unsigned int    minlevel,
+           const unsigned int    maxlevel);
+    
+                                    /**
+                                     * Virtual destructor.
+                                     */
+    virtual ~MGBase();
+    
+                                    /**
+                                     * Execute one step of the
+                                     * v-cycle algorithm.  This
+                                     * function assumes, that the
+                                     * vector #d# is properly filled
+                                     * with the residual in the outer
+                                     * defect correction
+                                     * scheme. After execution of
+                                     * #vcycle()#, the result is in
+                                     * the vector #s#. We propose to
+                                     * write functions #copy_from_mg#
+                                     * and #copy_to_mg# in derived
+                                     * classes of #MGBase# to access
+                                     * #d# and #s#.
+                                     *
+                                     * The actual work for this
+                                     * function is done in
+                                     * #level_mgstep#.
+                                     */
+    void vcycle(const MGSmootherBase     &pre_smooth,
+               const MGSmootherBase     &post_smooth,
+               const MGCoarseGridSolver &cgs);
 
                                     /**
-                                     * Level for coarse grid solution.
+                                     * Exception.
                                      */
-    unsigned int minlevel;
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
+    
+  protected:
                                     /**
-                                     * Auxiliary vector, defect.
+                                     * Highest level of cells.
                                      */
-    MGVector<Vector<double> > d;
+    unsigned int maxlevel;
 
                                     /**
-                                     * Auxiliary vector, solution.
+                                     * Level for coarse grid solution.
                                      */
-    MGVector<Vector<double> > s;
+    unsigned int minlevel;
+    
                                     /**
-                                     * Prolongation and restriction object.
+                                     * Auxiliary vector.
                                      */
-    SmartPointer<const MGTransferBase> transfer;
+    MGVector<Vector<double> > defect;
 
-  private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<double> t;
+    MGVector<Vector<double> > solution;
     
                                     /**
-                                     * Exception.
+                                     * Prolongation and restriction object.
                                      */
-    DeclException2(ExcSwitchedLevels, int, int,
-                  << "minlevel and maxlevel switched, should be: "
-                  << arg1 << "<=" << arg2);
-    
-  protected:
+    SmartPointer<const MGTransferBase> transfer;
 
                                     /**
                                      * The actual v-cycle multigrid method.
-                                     * This function is called on the
-                                     * highest level and recursively
-                                     * invokes itself down to the
-                                     * coarsest. There, it calls
-                                     * #coarse_grid_solver# and
-                                     * proceeds back up.
-                                     */
-    void level_mgstep(unsigned int level,
-                     const MGSmootherBase& pre_smooth,
-                     const MGSmootherBase& post_smooth,
-                     const MGCoarseGridSolver& cgs);  
+                                     * #level# is the level the
+                                     * function should work on. It
+                                     * will usually be called for the
+                                     * highest level from outside,
+                                     * but will then call itself
+                                     * recursively for #level-1#,
+                                     * unless we are on #minlevel#
+                                     * where instead of recursing
+                                     * deeper, the coarse grid solver
+                                     * is used to solve the matrix of
+                                     * this level.
+                                     */
+    void level_mgstep (const unsigned int        level,
+                      const MGSmootherBase     &pre_smooth,
+                      const MGSmootherBase     &post_smooth,
+                      const MGCoarseGridSolver &cgs);  
 
                                     /**
                                      * Apply negative #vmult# on all
@@ -361,39 +452,16 @@ class MGBase : public Subscriptor
                                      * This is implemented in a
                                      * derived class.
                                      */
-    virtual void level_vmult(unsigned int level,
-                            Vector<double>& dst,
-                            const Vector<double>& src,
-                            const Vector<double>& rhs) = 0;  
-
+    virtual void level_vmult (const unsigned int    level,
+                             Vector<double>       &dst,
+                             const Vector<double> &src,
+                             const Vector<double> &rhs) = 0;
   
-  public:
-                                    /**
-                                     * Constructor, subject to change.
-                                     */
-    MGBase(const MGTransferBase& transfer,
-                 unsigned int minlevel, unsigned int maxlevel);
-                                    /**
-                                     * Virtual destructor.
-                                     */
-    virtual ~MGBase();
+  private:
                                     /**
-                                     * Execute v-cycle algorithm.
-                                     * This function assumes, that
-                                     * the vector #d# is properly
-                                     * filled with the residual in
-                                     * the outer defect correction
-                                     * scheme. After execution of
-                                     * #vcycle()#, the result is in
-                                     * the vector #s#. We propose to
-                                     * write functions #copy_from_mg#
-                                     * and #copy_to_mg# in derived
-                                     * classes of #MGBase# to access
-                                     * #d# and #s#.
+                                     * Auxiliary vector.
                                      */
-    void vcycle(const MGSmootherBase& pre_smooth,
-               const MGSmootherBase& post_smooth,
-               const MGCoarseGridSolver& cgs);
+    Vector<double> t;    
 };
 
 
@@ -407,9 +475,10 @@ class MGBase : public Subscriptor
  * Furthermore, it needs functions #void copy_to_mg(const VECTOR&)#
  * to store #src# in the right hand side of the multi-level method and
  * #void copy_from_mg(VECTOR&)# to store the result of the v-cycle in #dst#.
+ *
  * @author Guido Kanschat, 1999
  */
-template<class MG, class VECTOR>
+template<class MG, class VECTOR = Vector<double> >
 class PreconditionMG
 {
   public:
@@ -419,29 +488,39 @@ class PreconditionMG
                                      * pre-smoother, post-smoother and
                                      * coarse grid solver.
                                      */
-    PreconditionMG(MG& mg,
-                  const MGSmootherBase& pre,
-                  const MGSmootherBase& post,
-                  const MGCoarseGridSolver& coarse);
+    PreconditionMG(MG                       &mg,
+                  const MGSmootherBase     &pre,
+                  const MGSmootherBase     &post,
+                  const MGCoarseGridSolver &coarse);
+    
                                     /**
-                                     * Preconditioning operator.
-                                     * This is the operator used by LAC
-                                     * iterative solvers.
+                                     * Preconditioning operator,
+                                     * calling the #vcycle# function
+                                     * of the #MG# object passed to
+                                     * the constructor.
+                                     *
+                                     * This is the operator used by
+                                     * LAC iterative solvers.
                                      */
-    void operator() (VECTOR& dst, const VECTOR& src) const;
+    void operator() (VECTOR       &dst,
+                    const VECTOR &src) const;
+    
   private:
                                     /**
-                                     * The mulrigrid object.
+                                     * The multigrid object.
                                      */
     SmartPointer<MG> multigrid;
+    
                                     /**
                                      * The pre-smoothing object.
                                      */
     SmartPointer<const MGSmootherBase> pre;
+
                                     /**
                                      * The post-smoothing object.
                                      */
     SmartPointer<const MGSmootherBase> post;
+    
                                     /**
                                      * The coarse grid solver.
                                      */
@@ -451,58 +530,44 @@ class PreconditionMG
 
 
 
-
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
-MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
-               :
-               solver(s), matrix(m), precondition(p)
-{}
-
-template<class SOLVER, class MATRIX, class PRECOND>
-void
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<double>& dst, const Vector<double>& src) const
-{
-  solver.solve(matrix, dst, src, precondition);
-}
-
-
-
-
 /* ------------------------------------------------------------------- */
 
 template<class VECTOR>
-MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+MGVector<VECTOR>::MGVector(const unsigned int min,
+                          const unsigned int max)
                :
-               vector<VECTOR>(max-min+1),
-               minlevel(min)
-{}
+               minlevel(min),
+               vectors(max-min+1)
+{};
+
 
 
 template<class VECTOR>
-VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i)
+VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i)
 {
-  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
-  return vector<VECTOR>::operator[](i-minlevel);
+  Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+        ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+  return vectors[i-minlevel];
 }
 
 
 template<class VECTOR>
-const VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i) const
+const VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i) const
 {
-  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
-  return vector<VECTOR>::operator[](i-minlevel);
+  Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+        ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+  return vectors[i-minlevel];
 }
 
+
 template<class VECTOR>
 void
 MGVector<VECTOR>::clear()
 {
   typename vector<VECTOR>::iterator v;
-  for (v = begin(); v != end(); ++v)
+  for (v = vectors.begin(); v != vectors.end(); ++v)
     v->clear();
   
 }
@@ -525,7 +590,8 @@ MGMatrix<MATRIX>::operator[](unsigned int i)
 {
   Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
   return vector<MATRIX>::operator[](i-minlevel);
-}
+};
+
 
 
 template<class MATRIX>
@@ -534,18 +600,49 @@ MGMatrix<MATRIX>::operator[](unsigned int i) const
 {
   Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
   return vector<MATRIX>::operator[](i-minlevel);
+};
+
+
+
+/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+MGCoarseGridLACIteration(SOLVER        &s,
+                        const MATRIX  &m,
+                        const PRECOND &p)
+               :
+               solver(s),
+               matrix(&m),
+               precondition(p)
+{};
+
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+void
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+operator() (const unsigned int    /* level */,
+           Vector<double>       &dst,
+           const Vector<double> &src) const
+{
+  solver.solve(*matrix, dst, src, precondition);
 }
 
 
 
+
+
+
 /* ------------------------------------------------------------------- */
 
 
 template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MGmg,
-                                          const MGSmootherBasepre,
-                                          const MGSmootherBasepost,
-                                          const MGCoarseGridSolvercoarse)
+PreconditionMG<MG, VECTOR>::PreconditionMG(MG                       &mg,
+                                          const MGSmootherBase     &pre,
+                                          const MGSmootherBase     &post,
+                                          const MGCoarseGridSolver &coarse)
                :
                multigrid(&mg),
                pre(&pre),
@@ -557,8 +654,8 @@ PreconditionMG<MG, VECTOR>::PreconditionMG(MG& mg,
 
 template<class MG, class VECTOR>
 void
-PreconditionMG<MG, VECTOR>::operator() (VECTOR& dst,
-                                       const VECTOR& src) const
+PreconditionMG<MG,VECTOR>::operator() (VECTOR       &dst,
+                                      const VECTOR &src) const
 {
   multigrid->copy_to_mg(src);
   multigrid->vcycle(*pre, *post, *coarse);
index 9df47ff5949ba4b6e4ce8af46e129b7cc60caba5..3662347101cb84a7124a0ff1dc5036512899aea4 100644 (file)
 #include <iostream>
 #include <cmath>
 
+
+
+//TODO: this function is only for debugging purposes and should be removed sometimes
 template<class VECTOR>
+static
 void print_vector(ostream& s, const VECTOR& v)
 {
   const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
   unsigned int k=0;
-  
+
+                                  // write the vector entries in a
+                                  // kind of square
   for (unsigned int i=0;i<n;++i)
     {
       for (unsigned int j=0;j<n;++j)
-//     s << ((float)i)/n << '\t' << ((float)j)/n << '\t' << v(k++) << endl;
        s << '\t' << v(k++);
       s << endl;
     }
   s << endl;
 }
 
+
+
 MGBase::~MGBase () 
 {};
 
 
+//////////////////////////////////////////////////////////////////////
+
+
+MGSmootherBase::~MGSmootherBase()
+{};
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+                           Vector<double>       &,
+                           const Vector<double> &) const
+{};
+
 
-MGBase::MGBase(const MGTransferBase& transfer,
-                            unsigned minlevel, unsigned maxlevel)
+
+//////////////////////////////////////////////////////////////////////
+
+
+
+MGBase::MGBase(const MGTransferBase &transfer,
+              const unsigned        minlevel,
+              const unsigned        maxlevel)
                :
-               maxlevel(maxlevel), minlevel(minlevel),
-               d(minlevel,maxlevel),
-               s(minlevel,maxlevel),
+               maxlevel(maxlevel),
+               minlevel(minlevel),
+               defect(minlevel,maxlevel),
+               solution(minlevel,maxlevel),
                transfer(&transfer)
 {
-  Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
-}
+  Assert(minlevel <= maxlevel,
+        ExcSwitchedLevels(minlevel, maxlevel));
+};
+
+
 
 void
-MGBase::vcycle(const MGSmootherBasepre_smooth,
-              const MGSmootherBasepost_smooth,
-              const MGCoarseGridSolvercoarse_grid_solver)
+MGBase::vcycle(const MGSmootherBase     &pre_smooth,
+              const MGSmootherBase     &post_smooth,
+              const MGCoarseGridSolver &coarse_grid_solver)
 {
   level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-}
+};
+
 
 
 void
-MGBase::level_mgstep(unsigned level,
-                    const MGSmootherBasepre_smooth,
-                    const MGSmootherBasepost_smooth,
-                    const MGCoarseGridSolvercoarse_grid_solver)
+MGBase::level_mgstep(const unsigned int        level,
+                    const MGSmootherBase     &pre_smooth,
+                    const MGSmootherBase     &post_smooth,
+                    const MGCoarseGridSolver &coarse_grid_solver)
 {
-  s[level] = 0.;
+  solution[level] = 0.;
   
   if (level == minlevel)
     {
-      coarse_grid_solver(level, s[level], d[level]);
+      coarse_grid_solver(level, solution[level], defect[level]);
       return;
     }
 
                           // smoothing of the residual by modifying s
-  pre_smooth.smooth(level, s[level], d[level]);
+  pre_smooth.smooth(level, solution[level], defect[level]);
                                   // t = d-As
-  t.reinit(s[level].size());
-  level_vmult(level, t, s[level], d[level]);
+  t.reinit(solution[level].size());
+  level_vmult(level, t, solution[level], defect[level]);
 
                                   // make t rhs of lower level
-  transfer->restrict(level, d[level-1], t);
+  transfer->restrict(level, defect[level-1], t);
   
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
                                   // do coarse grid correction
-  t.reinit(s[level].size());
-  transfer->prolongate(level, t, s[level-1]);
-  s[level].add(t);
+  t.reinit(solution[level].size());
+  transfer->prolongate(level, t, solution[level-1]);
+  solution[level].add(t);
   
                                   // smoothing (modify s again)
-  post_smooth.smooth(level, s[level], d[level]);
+  post_smooth.smooth(level, solution[level], defect[level]);
 }
 
 
 //////////////////////////////////////////////////////////////////////
 
 MGCoarseGridSolver::~MGCoarseGridSolver()
-{}
+{};
+
 
 
 //////////////////////////////////////////////////////////////////////
 
 MGTransferBase::~MGTransferBase()
-{}
-
-MGSmootherBase::~MGSmootherBase()
-{}
-
+{};
 
-//////////////////////////////////////////////////////////////////////
 
 
-void
-MGSmootherIdentity::smooth (const unsigned int,
-                           Vector<double>       &,
-                           const Vector<double> &) const
-{}
 
 

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.