]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The First MultiGrid Test works
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 4 May 1999 15:29:14 +0000 (15:29 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 4 May 1999 15:29:14 +0000 (15:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@1265 0785d39b-7218-0410-832d-ea1e28bc413d

15 files changed:
deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/mg_smoother.h
deal.II/deal.II/include/numerics/mg_smoother.h
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/deal.II/source/multigrid/mg_smoother.cc
deal.II/deal.II/source/numerics/mg_smoother.cc
deal.II/lac/include/lac/mgbase.h
deal.II/lac/include/lac/solver_cg.h
deal.II/lac/source/mgbase.cc
deal.II/lac/source/solver_control.cc
tests/lac/mg.cc
tests/lac/mgbase.cc
tests/lac/solver.cc
tests/lac/testmatrix.cc
tests/lac/testmatrix.h

index e0e83596d5f06da041c91d03d6d6f17b4cec62f8..c30ef1930c08153207aa5d1fb825471f3f17a6dd 100644 (file)
@@ -11,6 +11,8 @@
 
 #include <vector>
 
+#define FLOAT float
+
 
 /**
  * Abstract base class for coarse grid solvers.
@@ -31,8 +33,8 @@ class MGCoarseGridSolver
                                      * about the matrix is removed to
                                      * that class.
                                      */
-    virtual void operator() (unsigned int level, Vector<float>& dst,
-                            const Vector<float>& src) const = 0;
+    virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src) const = 0;
 };
 
 
@@ -65,9 +67,27 @@ class MGSmootherBase
                                      * things.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const = 0;
+                        Vector<FLOAT>       &u,
+                        const Vector<FLOAT> &rhs) const = 0;
+
+};
 
+/**
+ * Smoother doing nothing.
+ * @author Guido Kanschat, 1999
+ */
+class MGSmootherIdentity
+  :
+  public MGSmootherBase
+{
+  public:
+                                    /**
+                                     * Implementation of the interface in #MGSmootherBase#.
+                                     * This function does nothing.
+                                     */
+    virtual void smooth (const unsigned int   level,
+                        Vector<FLOAT>       &u,
+                        const Vector<FLOAT> &rhs) const;
 };
 
 /**
@@ -100,8 +120,8 @@ class MGTransferBase  :  public Subscriptor
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const = 0;
+                            Vector<FLOAT>       &dst,
+                            const Vector<FLOAT> &src) const = 0;
 
                                     /**
                                      * Restrict a vector from level
@@ -117,10 +137,44 @@ class MGTransferBase  :  public Subscriptor
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const = 0;
+                          Vector<FLOAT>       &dst,
+                          const Vector<FLOAT> &src) const = 0;
+};
+
+/**
+ * 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.
+ * @author Guido Kanschat, 1999
+ */
+template<class VECTOR>
+class MGVector
+  :
+  public Subscriptor,
+  public vector<VECTOR>
+{
+  public:
+                                    /**
+                                     * Constructor allowing to initialize the number of levels.
+                                     */
+    MGVector(unsigned int minlevel, unsigned int maxlevel);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    VECTOR& operator[](unsigned int);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    const VECTOR& operator[](unsigned int) const;
+  private:
+                                    /**
+                                     * Level of first component.
+                                     */
+    unsigned int minlevel;
 };
 
+
+
 /**
  * An array of matrices for each level.
  * This class provides the functionality of #vector<MATRIX># combined
@@ -137,7 +191,21 @@ class MGMatrix
                                     /**
                                      * Constructor allowing to initialize the number of levels.
                                      */
-    MGMatrix(unsigned int n_levels);
+    MGMatrix(unsigned int minlevel, unsigned int maxlevel);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    MATRIX& operator[](unsigned int);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    const MATRIX& operator[](unsigned int) const;
+    
+  private:
+                                    /**
+                                     * Level of first component.
+                                     */
+    unsigned int minlevel;
 };
 
 
@@ -168,8 +236,8 @@ class MGCoarseGridLACIteration
                                      * matrix, vectors and
                                      * preconditioner.
                                      */
-    virtual void operator() (unsigned int level, Vector<float>& dst,
-                            const Vector<float>& src) const;
+    virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src) const;
   private:
                                     /**
                                      * Reference to the solver.
@@ -207,23 +275,35 @@ class MGCoarseGridLACIteration
  * */
 class MGBase
 {
+  private:
     MGBase(const MGBase&);
     const MGBase& operator=(const MGBase&);
-    
+
+  protected:
+                                    /**
+                                     * Highest level of cells.
+                                     */
+    unsigned int maxlevel;
+
+                                    /**
+                                     * Level for coarse grid solution.
+                                     */
+    unsigned int minlevel;
                                     /**
                                      * Auxiliary vector, defect.
                                      */
-    vector<Vector<float> > d;
+    MGVector<Vector<FLOAT> > d;
 
                                     /**
                                      * Auxiliary vector, solution.
                                      */
-    vector<Vector<float> > s;
+    MGVector<Vector<FLOAT> > s;
 
+  private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<float> t;
+    Vector<FLOAT> t;
     
                                     /**
                                      * Prolongation and restriction object.
@@ -231,15 +311,6 @@ class MGBase
     SmartPointer<const MGTransferBase> transfer;
     
   protected:
-                                    /**
-                                     * Highest level of cells.
-                                     */
-    unsigned int maxlevel;
-
-                                    /**
-                                     * Level for coarse grid solution.
-                                     */
-    unsigned int minlevel;
 
                                     /**
                                      * The actual v-cycle multigrid method.
@@ -262,9 +333,9 @@ class MGBase
                                      * derived class.
                                      */
     virtual void level_residual(unsigned int level,
-                            Vector<float>& dst,
-                            const Vector<float>& src,
-                            const Vector<float>& rhs) = 0;  
+                            Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src,
+                            const Vector<FLOAT>& rhs) = 0;  
 
   
   public:
@@ -273,8 +344,27 @@ class MGBase
                                      */
     MGBase(const MGTransferBase& transfer,
                  unsigned int maxlevel, unsigned int minlevel);
+                                    /**
+                                     * Virtual destructor.
+                                     */
     virtual ~MGBase();
-    
+                                    /**
+                                     * 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#.
+                                     */
+    void vcycle(const MGSmootherBase& pre_smooth,
+               const MGSmootherBase& post_smooth,
+               const MGCoarseGridSolver& cgs);
 };
 
 
@@ -288,16 +378,64 @@ MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
 template<class SOLVER, class MATRIX, class PRECOND>
 void
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<FLOAT>& dst, const Vector<FLOAT>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
 
+//////////////////////////////////////////////////////////////////////
+
+template<class VECTOR>
+MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+               :
+               vector<VECTOR>(max-min+1),
+               minlevel(min)
+{}
+
+
+template<class VECTOR>
+VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i)
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<VECTOR>::operator[](i-minlevel);
+}
+
+
+template<class VECTOR>
+const VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i) const
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<VECTOR>::operator[](i-minlevel);
+}
+
+//////////////////////////////////////////////////////////////////////
+
 template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int nlevels)
+MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
                :
-               vector<MATRIX>(nlevels)
+               vector<MATRIX>(max-min+1),
+               minlevel(min)
 {}
 
 
+template<class MATRIX>
+MATRIX&
+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>
+const MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i) const
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
 #endif
index 1016d092c56ed727e4fa77ef367b121a8bce0b46..fb8fa27b4a7f2d7bce75702adcc204a6270cd379 100644 (file)
@@ -110,29 +110,6 @@ class MGSmoother
     vector<vector<int> > interior_boundary_dofs;
 };
 
-/**
- * Smoother doing nothing.
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity
-  :
-  public MGSmoother
-{
-  public:
-                                    /**
-                                     * Constructor.
-                                     */
-    template<int dim>
-    MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof);
-                                    /**
-                                     * Implementation of the interface in #MGSmootherBase#.
-                                     * This function does nothing.
-                                     */
-    virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
-};
-
 /**
  * Implementation of an SOR smoother.
  *
index 1016d092c56ed727e4fa77ef367b121a8bce0b46..fb8fa27b4a7f2d7bce75702adcc204a6270cd379 100644 (file)
@@ -110,29 +110,6 @@ class MGSmoother
     vector<vector<int> > interior_boundary_dofs;
 };
 
-/**
- * Smoother doing nothing.
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity
-  :
-  public MGSmoother
-{
-  public:
-                                    /**
-                                     * Constructor.
-                                     */
-    template<int dim>
-    MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof);
-                                    /**
-                                     * Implementation of the interface in #MGSmootherBase#.
-                                     * This function does nothing.
-                                     */
-    virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
-};
-
 /**
  * Implementation of an SOR smoother.
  *
index c6252fa0d8d1dafe6c3f5ecff1b52a573581a262..73a80470fb19336e436d37da737be686f8918290 100644 (file)
@@ -1,6 +1,24 @@
 // $Id$
 
 #include <lac/mgbase.h>
+#include <iostream>
+#include <cmath>
+
+template<class VECTOR>
+void print_vector(ostream& s, const VECTOR& v)
+{
+  const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
+  unsigned int k=0;
+  
+  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 () 
 {};
@@ -10,11 +28,21 @@ MGBase::~MGBase ()
 MGBase::MGBase(const MGTransferBase& transfer,
                             unsigned maxlevel, unsigned minlevel)
                :
-               d(maxlevel-minlevel),
-               s(maxlevel-minlevel),
-               transfer(&transfer),
-               maxlevel(maxlevel), minlevel(minlevel)
-{}
+               maxlevel(maxlevel), minlevel(minlevel),
+               d(minlevel,maxlevel),
+               s(minlevel,maxlevel),
+               transfer(&transfer)
+{
+}
+
+void
+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,
@@ -22,31 +50,51 @@ MGBase::level_mgstep(unsigned level,
                     const MGSmootherBase& post_smooth,
                     const MGCoarseGridSolver& coarse_grid_solver)
 {
+  s[level] = 0.;
+  
   if (level == minlevel)
     {
       coarse_grid_solver(level, s[level], d[level]);
       return;
     }
-  
-                                  // smoothing of the residual by modifying s
-  pre_smooth.smooth(level, s[level], d[level]);
 
+                          // smoothing of the residual by modifying s
+  pre_smooth.smooth(level, s[level], d[level]);
                                   // t = d-As
+  t.reinit(s[level].size());
   level_residual(level, t, s[level], d[level]);
-  
+
                                   // make t rhs of lower level
-  transfer->restrict(level, t, d[level-1]);
+  transfer->restrict(level, d[level-1], t);
+  
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
                                   // do coarse grid correction
-  transfer->prolongate(level, s[level], s[level-1]);
-
+  t.reinit(s[level].size());
+  transfer->prolongate(level, t, s[level-1]);
+  s[level].add(t);
+  
                                   // smoothing (modify s again)
   post_smooth.smooth(level, s[level], d[level]);
 }
 
+
+//////////////////////////////////////////////////////////////////////
+
 MGTransferBase::~MGTransferBase()
 {}
 
 MGSmootherBase::~MGSmootherBase()
 {}
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+                           Vector<FLOAT>       &,
+                           const Vector<FLOAT> &) const
+{}
+
+
index cc16654c94a5eadcf65a2d91d75389e3b4859b82..71b80cb4b3c4af2b7454baf00a08a7b011c6eea2 100644 (file)
@@ -107,22 +107,6 @@ MGSmoother::set_zero_interior_boundary (const unsigned int  level,
       u(*p) = 0;
 };
 
-//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof)
-               :
-               MGSmoother(mg_dof, 0)
-{}
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
-                      Vector<float>       &,
-                           const Vector<float> &) const
-{}
-
-
 //////////////////////////////////////////////////////////////////////
 
 template<int dim>
@@ -153,5 +137,3 @@ template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_
                                      const MGMatrix<SparseMatrix<float> >& matrix,
                                      unsigned int steps);
 
-template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<deal_II_dimension> &mg_dof);
-
index cc16654c94a5eadcf65a2d91d75389e3b4859b82..71b80cb4b3c4af2b7454baf00a08a7b011c6eea2 100644 (file)
@@ -107,22 +107,6 @@ MGSmoother::set_zero_interior_boundary (const unsigned int  level,
       u(*p) = 0;
 };
 
-//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof)
-               :
-               MGSmoother(mg_dof, 0)
-{}
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
-                      Vector<float>       &,
-                           const Vector<float> &) const
-{}
-
-
 //////////////////////////////////////////////////////////////////////
 
 template<int dim>
@@ -153,5 +137,3 @@ template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_
                                      const MGMatrix<SparseMatrix<float> >& matrix,
                                      unsigned int steps);
 
-template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<deal_II_dimension> &mg_dof);
-
index 3cb773a6b2a738a97771f6dc8304a26bb013aa97..54c15ce614ae563eca818dd2b2a92836309f0852 100644 (file)
@@ -11,6 +11,8 @@
 
 #include <vector>
 
+#define FLOAT float
+
 
 /**
  * Abstract base class for coarse grid solvers.
@@ -31,8 +33,8 @@ class MGCoarseGridSolver
                                      * about the matrix is removed to
                                      * that class.
                                      */
-    virtual void operator() (unsigned int level, Vector<float>& dst,
-                            const Vector<float>& src) const = 0;
+    virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src) const = 0;
 };
 
 
@@ -65,9 +67,27 @@ class MGSmootherBase
                                      * things.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const = 0;
+                        Vector<FLOAT>       &u,
+                        const Vector<FLOAT> &rhs) const = 0;
+
+};
 
+/**
+ * Smoother doing nothing.
+ * @author Guido Kanschat, 1999
+ */
+class MGSmootherIdentity
+  :
+  public MGSmootherBase
+{
+  public:
+                                    /**
+                                     * Implementation of the interface in #MGSmootherBase#.
+                                     * This function does nothing.
+                                     */
+    virtual void smooth (const unsigned int   level,
+                        Vector<FLOAT>       &u,
+                        const Vector<FLOAT> &rhs) const;
 };
 
 /**
@@ -100,8 +120,8 @@ class MGTransferBase  :  public Subscriptor
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const = 0;
+                            Vector<FLOAT>       &dst,
+                            const Vector<FLOAT> &src) const = 0;
 
                                     /**
                                      * Restrict a vector from level
@@ -117,10 +137,44 @@ class MGTransferBase  :  public Subscriptor
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const = 0;
+                          Vector<FLOAT>       &dst,
+                          const Vector<FLOAT> &src) const = 0;
+};
+
+/**
+ * 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.
+ * @author Guido Kanschat, 1999
+ */
+template<class VECTOR>
+class MGVector
+  :
+  public Subscriptor,
+  public vector<VECTOR>
+{
+  public:
+                                    /**
+                                     * Constructor allowing to initialize the number of levels.
+                                     */
+    MGVector(unsigned int minlevel, unsigned int maxlevel);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    VECTOR& operator[](unsigned int);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    const VECTOR& operator[](unsigned int) const;
+  private:
+                                    /**
+                                     * Level of first component.
+                                     */
+    unsigned int minlevel;
 };
 
+
+
 /**
  * An array of matrices for each level.
  * This class provides the functionality of #vector<MATRIX># combined
@@ -137,7 +191,21 @@ class MGMatrix
                                     /**
                                      * Constructor allowing to initialize the number of levels.
                                      */
-    MGMatrix(unsigned int n_levels);
+    MGMatrix(unsigned int minlevel, unsigned int maxlevel);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    MATRIX& operator[](unsigned int);
+                                    /**
+                                     * Safe access operator.
+                                     */
+    const MATRIX& operator[](unsigned int) const;
+    
+  private:
+                                    /**
+                                     * Level of first component.
+                                     */
+    unsigned int minlevel;
 };
 
 
@@ -168,8 +236,8 @@ class MGCoarseGridLACIteration
                                      * matrix, vectors and
                                      * preconditioner.
                                      */
-    virtual void operator() (unsigned int level, Vector<float>& dst,
-                            const Vector<float>& src) const;
+    virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src) const;
   private:
                                     /**
                                      * Reference to the solver.
@@ -207,23 +275,35 @@ class MGCoarseGridLACIteration
  * */
 class MGBase
 {
+  private:
     MGBase(const MGBase&);
     const MGBase& operator=(const MGBase&);
-    
+
+  protected:
+                                    /**
+                                     * Highest level of cells.
+                                     */
+    unsigned int maxlevel;
+
+                                    /**
+                                     * Level for coarse grid solution.
+                                     */
+    unsigned int minlevel;
                                     /**
                                      * Auxiliary vector, defect.
                                      */
-    vector<Vector<float> > d;
+    MGVector<Vector<FLOAT> > d;
 
                                     /**
                                      * Auxiliary vector, solution.
                                      */
-    vector<Vector<float> > s;
+    MGVector<Vector<FLOAT> > s;
 
+  private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<float> t;
+    Vector<FLOAT> t;
     
                                     /**
                                      * Prolongation and restriction object.
@@ -231,15 +311,6 @@ class MGBase
     SmartPointer<const MGTransferBase> transfer;
     
   protected:
-                                    /**
-                                     * Highest level of cells.
-                                     */
-    unsigned int maxlevel;
-
-                                    /**
-                                     * Level for coarse grid solution.
-                                     */
-    unsigned int minlevel;
 
                                     /**
                                      * The actual v-cycle multigrid method.
@@ -262,9 +333,9 @@ class MGBase
                                      * derived class.
                                      */
     virtual void level_residual(unsigned int level,
-                            Vector<float>& dst,
-                            const Vector<float>& src,
-                            const Vector<float>& rhs) = 0;  
+                            Vector<FLOAT>& dst,
+                            const Vector<FLOAT>& src,
+                            const Vector<FLOAT>& rhs) = 0;  
 
   
   public:
@@ -273,8 +344,27 @@ class MGBase
                                      */
     MGBase(const MGTransferBase& transfer,
                  unsigned int maxlevel, unsigned int minlevel);
+                                    /**
+                                     * Virtual destructor.
+                                     */
     virtual ~MGBase();
-    
+                                    /**
+                                     * 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#.
+                                     */
+    void vcycle(const MGSmootherBase& pre_smooth,
+               const MGSmootherBase& post_smooth,
+               const MGCoarseGridSolver& cgs);
 };
 
 
@@ -288,16 +378,64 @@ MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
 template<class SOLVER, class MATRIX, class PRECOND>
 void
 MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<FLOAT>& dst, const Vector<FLOAT>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
 
+//////////////////////////////////////////////////////////////////////
+
+template<class VECTOR>
+MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+               :
+               vector<VECTOR>(max-min+1),
+               minlevel(min)
+{}
+
+
+template<class VECTOR>
+VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i)
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<VECTOR>::operator[](i-minlevel);
+}
+
+
+template<class VECTOR>
+const VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i) const
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<VECTOR>::operator[](i-minlevel);
+}
+
+//////////////////////////////////////////////////////////////////////
+
 template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int nlevels)
+MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
                :
-               vector<MATRIX>(nlevels)
+               vector<MATRIX>(max-min+1),
+               minlevel(min)
 {}
 
 
+template<class MATRIX>
+MATRIX&
+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>
+const MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i) const
+{
+  Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+  return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
 #endif
index 8b3c58a7b0f49bbb36371002956e5b5620c233f3..a0000885a65a6759f058c793318a52915d043eb9 100644 (file)
@@ -158,7 +158,7 @@ SolverCG<Matrix,Vector>::solve (const Matrix &A,
       
       g.add(alpha,Ad);
       x.add(alpha,d );
-      res = sqrt(g*g);
+      res = g.l2_norm();
 
       conv = control().check(it,res);
       if (conv) break;
index c6252fa0d8d1dafe6c3f5ecff1b52a573581a262..73a80470fb19336e436d37da737be686f8918290 100644 (file)
@@ -1,6 +1,24 @@
 // $Id$
 
 #include <lac/mgbase.h>
+#include <iostream>
+#include <cmath>
+
+template<class VECTOR>
+void print_vector(ostream& s, const VECTOR& v)
+{
+  const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
+  unsigned int k=0;
+  
+  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 () 
 {};
@@ -10,11 +28,21 @@ MGBase::~MGBase ()
 MGBase::MGBase(const MGTransferBase& transfer,
                             unsigned maxlevel, unsigned minlevel)
                :
-               d(maxlevel-minlevel),
-               s(maxlevel-minlevel),
-               transfer(&transfer),
-               maxlevel(maxlevel), minlevel(minlevel)
-{}
+               maxlevel(maxlevel), minlevel(minlevel),
+               d(minlevel,maxlevel),
+               s(minlevel,maxlevel),
+               transfer(&transfer)
+{
+}
+
+void
+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,
@@ -22,31 +50,51 @@ MGBase::level_mgstep(unsigned level,
                     const MGSmootherBase& post_smooth,
                     const MGCoarseGridSolver& coarse_grid_solver)
 {
+  s[level] = 0.;
+  
   if (level == minlevel)
     {
       coarse_grid_solver(level, s[level], d[level]);
       return;
     }
-  
-                                  // smoothing of the residual by modifying s
-  pre_smooth.smooth(level, s[level], d[level]);
 
+                          // smoothing of the residual by modifying s
+  pre_smooth.smooth(level, s[level], d[level]);
                                   // t = d-As
+  t.reinit(s[level].size());
   level_residual(level, t, s[level], d[level]);
-  
+
                                   // make t rhs of lower level
-  transfer->restrict(level, t, d[level-1]);
+  transfer->restrict(level, d[level-1], t);
+  
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
                                   // do coarse grid correction
-  transfer->prolongate(level, s[level], s[level-1]);
-
+  t.reinit(s[level].size());
+  transfer->prolongate(level, t, s[level-1]);
+  s[level].add(t);
+  
                                   // smoothing (modify s again)
   post_smooth.smooth(level, s[level], d[level]);
 }
 
+
+//////////////////////////////////////////////////////////////////////
+
 MGTransferBase::~MGTransferBase()
 {}
 
 MGSmootherBase::~MGSmootherBase()
 {}
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+                           Vector<FLOAT>       &,
+                           const Vector<FLOAT> &) const
+{}
+
+
index 5fc7b74c64965d3dce4483f51c182e4eba2427c3..85fa0e5b349b19897cf31b8991a3e2a1078af772 100644 (file)
@@ -107,7 +107,12 @@ ReductionControl::check (const unsigned int step,
     };
 
   if (check_value < reduced_tol)
-    return success;
+    {
+      if (log_result)
+       deallog << "Convergence step " << step
+               << " value " << check_value << endl;
+      return success;
+    }
   else
     return SolverControl::check(step, check_value);
 };
index 439d47f20df529bb89ba719475d8c3e4c33ec4f1..46995507745c84c537785a3894b1632a5223efaf 100644 (file)
 #include <lac/precondition.h>
 #include <lac/mgbase.h>
 
+template<class VECTOR>
+void print_vector(ostream& s, const VECTOR& v)
+{
+  const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
+  unsigned int k=0;
+  
+  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;
+}
+
 class FDMG
   :
   public MGBase
 {
   private:
-    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+                                    /**
+                                     * Pointer to the level matrices.
+                                     */
+    SmartPointer<MGMatrix<SparseMatrix<FLOAT> > >matrices;
   public:
-    FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
-        FDMGTransfer& transfer)
-                   :
-                   MGBase(transfer, maxlevel, 0),
-                   matrices(&matrices)
-      {}
+    FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<FLOAT> >& matrices,
+        FDMGTransfer& transfer);
     
     
     virtual void level_residual (unsigned int level,
-                                Vector<float>& dst,
-                                const Vector<float>& src,
-                                const Vector<float>& rhs);
+                                Vector<FLOAT>& dst,
+                                const Vector<FLOAT>& src,
+                                const Vector<FLOAT>& rhs);
+
+    void copy_to_mg(const Vector<double>& rhs);
+    void copy_from_mg(Vector<double>& lsg);
 };
 
 class MGSmootherLAC
   :
-  MGSmootherBase
+  public MGSmootherBase
 {
   private:
-    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+    SmartPointer<MGMatrix<SparseMatrix<FLOAT> > >matrices;
   public:
-    MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+    MGSmootherLAC(MGMatrix<SparseMatrix<FLOAT> >&);
     
     virtual void smooth (const unsigned int level,
-                        Vector<float> &u,
-                        const Vector<float> &rhs) const;
+                        Vector<FLOAT> &u,
+                        const Vector<FLOAT> &rhs) const;
     
 };
 
-typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float>  >,
-SparseMatrix<float>, PreconditionRelaxation<SparseMatrix<float> , Vector<float> > >
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<FLOAT> , Vector<FLOAT>  >,
+SparseMatrix<FLOAT>, /*PreconditionRelaxation<SparseMatrix<FLOAT> ,*/
+  PreconditionIdentity<Vector<FLOAT> > >
 Coarse;
 
 class MGPrecondition
 {
   private:
     FDMG& mg;
-    MGSmootherLAC& smooth;
+    MGSmootherBase& smooth;
     Coarse& coarse;
   public:
-    MGPrecondition(FDMG& m, MGSmootherLAC& s, Coarse& c)
+    MGPrecondition(FDMG& m, MGSmootherBase& s, Coarse& c)
                    :
                    mg(m), smooth(s), coarse(c)
       {}
@@ -72,7 +91,11 @@ class MGPrecondition
     
     void operator () (Vector<double>& dst, const Vector<double>& src) const
       {
+//     print_vector(cout, src);
        
+       mg.copy_to_mg(src);
+       mg.vcycle(smooth,smooth,coarse);
+       mg.copy_from_mg(dst);
       }
     
     
@@ -88,7 +111,7 @@ main()
   SolverControl control(100, 1.e-5, true);
 
   const unsigned int base = 3;
-  const unsigned int maxlevel = 6;
+  const unsigned int maxlevel = 8;
   
   const unsigned int maxsize = base * (1 << maxlevel);
   
@@ -96,50 +119,102 @@ main()
   FDMGTransfer transfer(maxsize, maxsize, maxlevel);
 
                                   // coarse grid solver
-  PrimitiveVectorMemory<Vector<float> > cgmem;
-  SolverControl cgcontrol(100, 1.e-12);
-  SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+  PrimitiveVectorMemory<Vector<FLOAT> > cgmem;
+  ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false);
+  SolverCG<SparseMatrix<FLOAT> , Vector<FLOAT> > cgcg(cgcontrol,cgmem);
 
-  vector<SparseMatrixStruct >  structure(maxlevel+1);
-  MGMatrix<SparseMatrix<float> > A(maxlevel+1);
   
   for (unsigned int level = 0; level <= maxlevel; ++level)
     {
+      const unsigned int minlevel = 0;
+
       const unsigned int size = base * (1 << level);
 
-      const unsigned int dim = (size+1)*(size+1);
+      const unsigned int dim = (size-1)*(size-1);
+      deallog << "Level " << level << " size " << size << endl;
 
                                       // Make matrix
+      vector<SparseMatrixStruct >  structure(maxlevel+1);
+      MGMatrix<SparseMatrix<FLOAT> > A(minlevel,maxlevel);
+
       FDMatrix testproblem(size, size);
-      structure[level].reinit(dim, dim, 5);
-      testproblem.build_structure(structure[level]);
-      structure[level].compress();
-      A[level].reinit(structure[level]);
-      testproblem.laplacian(A[level]);
 
+      for(unsigned ll = minlevel; ll <= level; ++ll)
+       {
+         const unsigned int size = base * (1 << ll);
+         const unsigned int dim = (size-1)*(size-1);
+         FDMatrix testproblem(size, size);
+         structure[ll].reinit(dim, dim, 5);
+         testproblem.build_structure(structure[ll]);
+         structure[ll].compress();
+         A[ll].reinit(structure[ll]);
+         testproblem.laplacian(A[ll]);
+       }
+      
       FDMG multigrid(level, A, transfer);
 
-      PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
-       cgprec(A[0], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+//      PreconditionRelaxation<SparseMatrix<FLOAT> , Vector<FLOAT> >
+      PreconditionIdentity<Vector<FLOAT> >
+       cgprec;//(A[minlevel], &SparseMatrix<FLOAT> ::template precondition_SSOR<FLOAT>, 1.2);
       
-      Coarse coarsegrid(cgcg, A[0], cgprec);
+      Coarse coarsegrid(cgcg, A[minlevel], cgprec);
       
       MGSmootherLAC smoother(A);
-      
-//      MGPrecondition prec(multigrid);
+//      MGSmootherIdentity smoother;
+
+      MGPrecondition precondition(multigrid, smoother, coarsegrid);
+
+//      SolverRichardson<SparseMatrix<FLOAT> , Vector<double> > solver(control, mem);
+      SolverCG<SparseMatrix<FLOAT> , Vector<double> > solver(control, mem);
+
+      Vector<double> u(dim);
+      Vector<double> f(dim);
+      u = 0.;
+      f = 1.;//(size*size);
+
+      solver.solve(A[level], u, f, precondition);
+      ofstream out("T");
+      testproblem.gnuplot_print(out,u);
     }
 }
 
+FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<FLOAT> >& matrices,
+          FDMGTransfer& transfer)
+               :
+               MGBase(transfer, maxlevel, 0),
+               matrices(&matrices)
+{
+  for (unsigned int level = minlevel; level<=maxlevel ; ++level)
+    {
+      s[level].reinit(matrices[level].m());
+      d[level].reinit(matrices[level].m());
+    }
+}
+
+
 void
 FDMG::level_residual (unsigned int level,
-                     Vector<float>& dst,
-                     const Vector<float>& src,
-                     const Vector<float>& rhs)
+                     Vector<FLOAT>& dst,
+                     const Vector<FLOAT>& src,
+                     const Vector<FLOAT>& rhs)
 {
   (*matrices)[level].residual(dst, src, rhs);
 }
 
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+void
+FDMG::copy_to_mg(const Vector<double>& v)
+{
+  d[maxlevel] = v;
+}
+
+void
+FDMG::copy_from_mg(Vector<double>& v)
+{
+  v = s[maxlevel];
+}
+
+
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<FLOAT> >& matrix)
                :
                matrices(&matrix)
 {}
@@ -147,14 +222,14 @@ MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
 
 void
 MGSmootherLAC::smooth (const unsigned int level,
-                      Vector<float> &u,
-                      const Vector<float> &rhs) const
+                      Vector<FLOAT> &u,
+                      const Vector<FLOAT> &rhs) const
 {
-  SolverControl control(0.,1);
-  PrimitiveVectorMemory<Vector<float> > mem;
-  SolverRichardson<SparseMatrix<float> , Vector<float>  > rich(control, mem);
-  PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
-    prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+  SolverControl control(1,1.e-300,false,false);
+  PrimitiveVectorMemory<Vector<FLOAT> > mem;
+  SolverRichardson<SparseMatrix<FLOAT> , Vector<FLOAT>  > rich(control, mem);
+  PreconditionRelaxation<SparseMatrix<FLOAT> , Vector<FLOAT> >
+    prec((*matrices)[level], &SparseMatrix<FLOAT> ::template precondition_SSOR<FLOAT>, 1.);
 
   rich.solve((*matrices)[level], u, rhs, prec);
 }
index acd703a6c856047893e31e44c12211220f92dc47..59d7287ba2ccbac9fa139da07bb6b1bf683c1073 100644 (file)
@@ -47,11 +47,17 @@ class MGTest
     MGCoarseGridSolver& solver;
 
   public:
-    MGTest(TransferTest& t, SmoothTest& s, MGCoarseGridTest& c)
+    MGTest(TransferTest& t, SmoothTest& sm, MGCoarseGridTest& cs)
                    : MGBase(t, 9, 3),
-                     smoother(s),
-                     solver(c)
-      {}
+                     smoother(sm),
+                     solver(cs)
+      {
+       for (unsigned int l = minlevel; l <= maxlevel; ++l)
+         {
+           d[l].reinit(1);
+           s[l].reinit(1);
+         }
+      }
     
     virtual void level_residual(unsigned level,
                                Vector<float>& dst,
index acaa7ed9b83accea4285fe16e31780f08d0bc36c..62d6bcd8420f5c23ab886780b3b6e5a0dfc72954 100644 (file)
@@ -42,7 +42,7 @@ main()
     {
       deallog << "Size " << size << endl;
       
-      unsigned int dim = (size+1)*(size+1);
+      unsigned int dim = (size-1)*(size-1);
 
                                       // Make matrix
       FDMatrix testproblem(size, size);
@@ -73,6 +73,7 @@ main()
 
       check_method(rich,A,u,f,prec_ssor);
       check_method(cg,A,u,f,prec_ssor);
+//      testproblem.gnuplot_print(cout, u);
       check_method(bicgstab,A,u,f,prec_ssor);
       check_method(gmres,A,u,f,prec_ssor);
 
index 8b22d664e731ffc21bde3c2c1117d18d490e6120..53eeee5063c91473a8e2da683759613b821fc786 100644 (file)
@@ -11,32 +11,32 @@ FDMatrix::FDMatrix(unsigned int nx, unsigned int ny)
 void
 FDMatrix::build_structure(SparseMatrixStruct& structure) const
 {
-  for(unsigned int i=0;i<=ny;i++)
+  for(unsigned int i=0;i<=ny-2;i++)
     {
-      for(unsigned int j=0;j<=nx; j++)
+      for(unsigned int j=0;j<=nx-2; j++)
        {
                                           // Number of the row to be entered
-         unsigned int row = j+(nx+1)*i;
+         unsigned int row = j+(nx-1)*i;
          structure.add(row, row);
          if (j>0)
            {
              structure.add(row-1, row);
              structure.add(row, row-1);
            }
-         if (j<nx)
+         if (j<nx-2)
            {
              structure.add(row+1, row);
              structure.add(row, row+1);
            }
          if (i>0)
            {
-             structure.add(row-nx, row);
-             structure.add(row, row-nx);
+             structure.add(row-(nx-1), row);
+             structure.add(row, row-(nx-1));
            }
-         if (i<ny)
+         if (i<ny-2)
            {
-             structure.add(row+nx, row);
-             structure.add(row, row+nx);
+             structure.add(row+(nx-1), row);
+             structure.add(row, row+(nx-1));
            }
        }
     }
@@ -46,12 +46,12 @@ template<typename number>
 void
 FDMatrix::laplacian(SparseMatrix<number>& A) const
 {
-   for(unsigned int i=0;i<=ny;i++)
+   for(unsigned int i=0;i<=ny-2;i++)
     {
-      for(unsigned int j=0;j<=nx; j++)
+      for(unsigned int j=0;j<=nx-2; j++)
        {
                                           // Number of the row to be entered
-         unsigned int row = j+(nx+1)*i;
+         unsigned int row = j+(nx-1)*i;
          
          A.set(row, row, 4.);
          if (j>0)
@@ -59,29 +59,46 @@ FDMatrix::laplacian(SparseMatrix<number>& A) const
              A.set(row-1, row, -1.);
              A.set(row, row-1, -1.);
            }
-         if (j<nx)
+         if (j<nx-2)
            {
              A.set(row+1, row, -1.);
              A.set(row, row+1, -1.);
            }
          if (i>0)
            {
-             A.set(row-nx, row, -1.);
-             A.set(row, row-nx, -1.);
+             A.set(row-(nx-1), row, -1.);
+             A.set(row, row-(nx-1), -1.);
            }
-         if (i<ny)
+         if (i<ny-2)
            {
-             A.set(row+nx, row, -1.);
-             A.set(row, row+nx, -1.);
+             A.set(row+(nx-1), row, -1.);
+             A.set(row, row+(nx-1), -1.);
            }
        }
     } 
 }
 
+template<typename number>
+void
+FDMatrix::gnuplot_print(ostream& s, const Vector<number>& V) const
+{
+   for(unsigned int i=0;i<=ny-2;i++)
+    {
+      for(unsigned int j=0;j<=nx-2; j++)
+       {
+                                          // Number of the row to be entered
+         unsigned int row = j+(nx-1)*i;
+         s << (j+1)/(float)nx << '\t' << (i+1)/(float)ny << '\t' << V(row) << endl;
+       }
+      s << endl;
+    } 
+   s << endl;
+}
+
 FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
                           unsigned int nlevels)
                :
-               structures(nlevels-1), matrices(nlevels-1)
+               structures(nlevels), matrices(nlevels)
 {
   unsigned int power = 1 << nlevels;
   
@@ -91,7 +108,7 @@ FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
   nx /= power;
   ny /= power;
   
-  for (unsigned int level = 0; level < nlevels-1; ++ level)
+  for (unsigned int level = 0; level < nlevels; ++ level)
     {
       build_matrix(nx,ny,structures[level],matrices[level]);
       nx *= 2;
@@ -101,18 +118,18 @@ FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
 
 void
 FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
-                          SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+                          SparseMatrixStruct& structure, SparseMatrix<FLOAT>& matrix)
 {
-  structure.reinit((nx+1)*(ny+1),(2*nx+1)*(2*ny+1),9);
+  structure.reinit((nx-1)*(ny-1),(2*nx-1)*(2*ny-1),9);
   
                                   // Traverse all points of coarse grid
-  for (unsigned int i=0 ; i<=ny; ++i)
-    for (unsigned int j=0 ; j<=nx; ++j)
+  for (unsigned int i=1 ; i<ny; ++i)
+    for (unsigned int j=1 ; j<nx; ++j)
       {
                                         // coarse grid point number
-       unsigned int ncoarse =j+(nx+1)*i;
+       unsigned int ncoarse =j+(nx-1)*i -(nx-1) -1;
                                         // same point on fine grid
-       unsigned int nfine = 2*j+(4*nx+2)*i;
+       unsigned int nfine = 2*j+(4*nx-2)*i -(2*nx-1) -1;
        
        structure.add(ncoarse,nfine);
                                         // left
@@ -121,10 +138,10 @@ FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
            structure.add(ncoarse,nfine-1);
                                             // lower left
            if (i>0)
-             structure.add(ncoarse,nfine-(2*nx+1)-1);
+             structure.add(ncoarse,nfine-(2*nx-1)-1);
                                             // upper left
            if (i<ny)
-             structure.add(ncoarse,nfine+(2*nx+1)-1);
+             structure.add(ncoarse,nfine+(2*nx-1)-1);
          }
                                         // right
        if (j<nx)
@@ -132,30 +149,30 @@ FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
            structure.add(ncoarse, nfine+1);
                                             // lower right
            if (i>0)
-             structure.add(ncoarse,nfine-(2*nx+1)+1);
+             structure.add(ncoarse,nfine-(2*nx-1)+1);
                                             // upper right
            if (i<ny)
-             structure.add(ncoarse,nfine+(2*nx+1)+1);
+             structure.add(ncoarse,nfine+(2*nx-1)+1);
          }
 
                                             // lower
            if (i>0)
-             structure.add(ncoarse,nfine-(2*nx+1));
+             structure.add(ncoarse,nfine-(2*nx-1));
                                             // upper
            if (i<ny)
-             structure.add(ncoarse,nfine+(2*nx+1));
+             structure.add(ncoarse,nfine+(2*nx-1));
       }
 
   structure.compress();
   matrix.reinit(structure);
   
-  for (unsigned int i=0 ; i<=ny; ++i)
-    for (unsigned int j=0 ; j<=nx; ++j)
+  for (unsigned int i=1 ; i<ny; ++i)
+    for (unsigned int j=1 ; j<nx; ++j)
       {
                                         // coarse grid point number
-       unsigned int ncoarse =j+(nx+1)*i;
+       unsigned int ncoarse =j+(nx-1)*i -(nx-1) -1;
                                         // same point on fine grid
-       unsigned int nfine = 2*j+(4*nx+2)*i;
+       unsigned int nfine =2*j+(4*nx-2)*i -(2*nx-1) -1;
        
        matrix.set(ncoarse,nfine,1.);
                                         // left
@@ -164,10 +181,10 @@ FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
            matrix.set(ncoarse,nfine-1,.5);
                                             // lower left
            if (i>0)
-             matrix.set(ncoarse,nfine-(2*nx+1)-1,.25);
+             matrix.set(ncoarse,nfine-(2*nx-1)-1,.25);
                                             // upper left
            if (i<ny)
-             matrix.set(ncoarse,nfine+(2*nx+1)-1,.25);
+             matrix.set(ncoarse,nfine+(2*nx-1)-1,.25);
          }
                                         // right
        if (j<nx)
@@ -175,26 +192,26 @@ FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
            matrix.set(ncoarse,nfine+1,.5);
                                             // lower right
            if (i>0)
-             matrix.set(ncoarse,nfine-(2*nx+1)+1,.25);
+             matrix.set(ncoarse,nfine-(2*nx-1)+1,.25);
                                             // upper right
            if (i<ny)
-             matrix.set(ncoarse,nfine+(2*nx+1)+1,.25);
+             matrix.set(ncoarse,nfine+(2*nx-1)+1,.25);
          }
 
                                             // lower
            if (i>0)
-             matrix.set(ncoarse,nfine-(2*nx+1),.5);
+             matrix.set(ncoarse,nfine-(2*nx-1),.5);
                                             // upper
            if (i<ny)
-             matrix.set(ncoarse,nfine+(2*nx+1),.5);
+             matrix.set(ncoarse,nfine+(2*nx-1),.5);
       }
   
 }
 
 void
 FDMGTransfer::prolongate (const unsigned int   to_level,
-                         Vector<float>       &dst,
-                         const Vector<float> &src) const
+                         Vector<FLOAT>       &dst,
+                         const Vector<FLOAT> &src) const
 {
   Assert((to_level>0) && (to_level<=matrices.size()),
         ExcIndexRange(to_level, 0, matrices.size()+1));
@@ -205,8 +222,8 @@ FDMGTransfer::prolongate (const unsigned int   to_level,
 
 void
 FDMGTransfer::restrict (const unsigned int   from_level,
-                       Vector<float>       &dst,
-                       const Vector<float> &src) const
+                       Vector<FLOAT>       &dst,
+                       const Vector<FLOAT> &src) const
 {
   Assert((from_level>0) && (from_level<=matrices.size()),
         ExcIndexRange(from_level, 0, matrices.size()+1));
@@ -217,3 +234,6 @@ FDMGTransfer::restrict (const unsigned int   from_level,
 
 template void FDMatrix::laplacian(SparseMatrix<double>&) const;
 template void FDMatrix::laplacian(SparseMatrix<float>&) const;
+template void FDMatrix::gnuplot_print(ostream& s, const Vector<double>& V) const;
+template void FDMatrix::gnuplot_print(ostream& s, const Vector<float>& V) const;
+
index 641f949e0fa0438d04f9b50a5b84d703dad99a70..659ab397527cf4f76e913ddf45aeab806909aca9 100644 (file)
@@ -28,6 +28,9 @@ class FDMatrix
     template <typename number>
     void laplacian(SparseMatrix<number>&) const;
 
+    template <typename number>
+    void gnuplot_print(ostream&, const Vector<number>&) const;
+    
   private:
                                     /**
                                      * Number of gridpoints in x-direction.
@@ -63,16 +66,16 @@ class FDMGTransfer
                                      * function in #MGTranferBase#.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const;
+                            Vector<FLOAT>       &dst,
+                            const Vector<FLOAT> &src) const;
 
                                     /**
                                      * Implementation of abstract
                                      * function in #MGTranferBase#.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const;
+                          Vector<FLOAT>       &dst,
+                          const Vector<FLOAT> &src) const;
 
                                     /**
                                      * Exception.
@@ -88,7 +91,7 @@ class FDMGTransfer
                                     /**
                                      * Prolongation matrices.
                                      */
-    vector<SparseMatrix<float> > matrices;
+    vector<SparseMatrix<FLOAT> > matrices;
 
                                     /**
                                      * Matrix generator.
@@ -98,5 +101,5 @@ class FDMGTransfer
                                      * fine to coarse (#vmult#).
                                      */
     void build_matrix(unsigned int nx, unsigned int ny,
-                     SparseMatrixStruct&, SparseMatrix<float>&);
+                     SparseMatrixStruct&, SparseMatrix<FLOAT>&);
 };

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.