]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
MG works on locally refined meshes, needs improvement and cleaning
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 28 Jun 1999 15:19:09 +0000 (15:19 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 28 Jun 1999 15:19:09 +0000 (15:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@1486 0785d39b-7218-0410-832d-ea1e28bc413d

21 files changed:
deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_constraints.templates.h
deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/mg_smoother.h
deal.II/deal.II/include/multigrid/mg_smoother.templates.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/include/numerics/data_io.h
deal.II/deal.II/include/numerics/mg_smoother.h
deal.II/deal.II/include/numerics/mg_smoother.templates.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/multigrid.h
deal.II/deal.II/include/numerics/multigrid.templates.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/deal.II/source/multigrid/mg_smoother.cc
deal.II/deal.II/source/multigrid/multigrid.cc
deal.II/deal.II/source/numerics/mg_smoother.cc
deal.II/deal.II/source/numerics/multigrid.cc
deal.II/lac/include/lac/mgbase.h
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/source/mgbase.cc

index 92a14f7bc83f9574e0afca5ea0642687bff1b117..34fb1ab7d81e0f803b640ac4e8151190bdee8226 100644 (file)
@@ -247,6 +247,13 @@ class ConstraintMatrix {
     template<typename number>
     void distribute (Vector<number> &vec) const;
     
+                                    /**
+                                     * Delete hanging nodes in a vector.
+                                     * Sets all hanging node values to zero.
+                                     */
+    template<typename number>
+    void set_zero (Vector<number> &vec) const;
+
     
                                     /**
                                      * Print the constraint lines. Mainly for
index 4e964898b8bcd733cd5e57e1b20dd40674b5cbf2..5c2f975c32fc1e7a85ced7ee0b66061e60a31cc0 100644 (file)
@@ -348,6 +348,32 @@ ConstraintMatrix::condense (Vector<number> &vec) const
   
 
 
+template<typename number>
+void
+ConstraintMatrix::set_zero (Vector<number> &vec) const
+{
+  Assert (sorted == true, ExcMatrixNotClosed());
+
+  if (lines.size() == 0)
+                                    // nothing to do
+    return;
+  
+  vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+  for (unsigned int row=0; row<vec.size(); ++row)
+    if (row == next_constraint->line)
+      {
+                                        // set entry to zero
+       vec(row) = 0.;
+       
+       ++next_constraint;
+       if (next_constraint == lines.end())
+                                          // nothing more to do
+         break;
+      };
+};
+  
+
+
 template<typename number>
 void
 ConstraintMatrix::distribute (const Vector<number> &condensed,
index 18704a4799f8bc75d76087f4a215c7a17da1a7d7..fd6a230c29eea6fa93cf872cfb394dcb5fc19d5a 100644 (file)
@@ -40,8 +40,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<double>& dst,
+                            const Vector<double>& src) const = 0;
 };
 
 
@@ -74,8 +74,8 @@ class MGSmootherBase
                                      * things.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const = 0;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const = 0;
 
 };
 
@@ -93,8 +93,8 @@ class MGSmootherIdentity
                                      * This function does nothing.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const;
 };
 
 /**
@@ -129,8 +129,8 @@ class MGTransferBase
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const = 0;
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const = 0;
 
                                     /**
                                      * Restrict a vector from level
@@ -146,8 +146,8 @@ class MGTransferBase
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const = 0;
+                          Vector<double>       &dst,
+                          const Vector<double> &src) const = 0;
 };
 
 /**
@@ -175,6 +175,10 @@ class MGVector
                                      * Safe access operator.
                                      */
     const VECTOR& operator[](unsigned int) const;
+                                    /**
+                                     * Delete all entries.
+                                     */
+    void clear();
   private:
                                     /**
                                      * Level of first component.
@@ -245,8 +249,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<double>& dst,
+                            const Vector<double>& src) const;
   private:
                                     /**
                                      * Reference to the solver.
@@ -303,23 +307,23 @@ class MGBase
                                     /**
                                      * Auxiliary vector, defect.
                                      */
-    MGVector<Vector<float> > d;
+    MGVector<Vector<double> > d;
 
                                     /**
                                      * Auxiliary vector, solution.
                                      */
-    MGVector<Vector<float> > s;
+    MGVector<Vector<double> > s;
+                                    /**
+                                     * Prolongation and restriction object.
+                                     */
+    SmartPointer<const MGTransferBase> transfer;
 
   private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<float> t;
+    Vector<double> t;
     
-                                    /**
-                                     * Prolongation and restriction object.
-                                     */
-    SmartPointer<const MGTransferBase> transfer;
                                     /**
                                      * Exception.
                                      */
@@ -344,15 +348,15 @@ class MGBase
                      const MGCoarseGridSolver& cgs);  
 
                                     /**
-                                     * Apply residual operator on all
+                                     * Apply negative #vmult# on all
                                      * cells of a level.
                                      * This is implemented in a
                                      * derived class.
                                      */
-    virtual void level_residual(unsigned int level,
-                            Vector<float>& dst,
-                            const Vector<float>& src,
-                            const Vector<float>& rhs) = 0;  
+    virtual void level_vmult(unsigned int level,
+                            Vector<double>& dst,
+                            const Vector<double>& src,
+                            const Vector<double>& rhs) = 0;  
 
   
   public:
@@ -450,7 +454,7 @@ MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
 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<double>& dst, const Vector<double>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
@@ -485,6 +489,15 @@ MGVector<VECTOR>::operator[](unsigned int i) const
   return vector<VECTOR>::operator[](i-minlevel);
 }
 
+template<class VECTOR>
+void
+MGVector<VECTOR>::clear()
+{
+  vector<VECTOR>::iterator v;
+  for (v = begin(); v != end(); ++v)
+    v->clear();
+  
+}
 
 
 /* ------------------------------------------------------------------- */
index 6db42e8039943910de355993d7893efca366fe1f..110cca5eb1cdbe4b2348274abc98f3c6292b9118 100644 (file)
@@ -22,8 +22,7 @@
  *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGSmoother
-  :  public MGSmootherBase
+class MGSmoother : public MGSmootherBase
 {
   private:
                                     /**
@@ -82,7 +81,7 @@ class MGSmoother
                                      * function does nothing in this case.
                                      */
     void set_zero_interior_boundary (const unsigned int  level,
-                                    Vector<float>      &u) const;
+                                    Vector<double>      &u) const;
                                     /**
                                      * Modify the number of smoothing steps.
                                      */
@@ -113,31 +112,45 @@ class MGSmoother
 };
 
 /**
- * Implementation of an SOR smoother.
- *
- * This class is an implementation of a smootin algorithm for multigrid. It
- * knows the actual instance of the multigrid matrix and uses the SOR method of the level
- * matrices.
+ * Implementation of a smoother using matrix builtin relaxation methods.
+ * 
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGSmootherSOR
-  :
-  public MGSmoother
+template<typename number>
+class MGSmootherRelaxation : public MGSmoother
 {
   public:
+                                    /**
+                                     * Type of the smoothing
+                                     * function of the matrix.
+                                     */
+    typedef void
+    (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
+                                          typename SparseMatrix<number>::value_type) const;
+    
                                     /**
                                      * Constructor.
                                      * The constructor uses an #MGDoFHandler# to initialize
                                      * data structures for interior level boundary handling
-                                     * in #MGSmoother#. Furthermore, it takes a pointer to the
-                                     * matrices and the number of SOR steps required in each
-                                     * smoothing operation.
+                                     * in #MGSmoother#.
+                                     *
+                                     * Furthermore, it takes a pointer to the
+                                     * level matrices and their smoothing function.
+                                     * This function must perform one relaxation step
+                                     * like #SparseMatrix<number>::SOR_step# does. Do not
+                                     * use the preconditioning methods because they apply
+                                     * a preconditioning matrix to the residual.
+                                     *
+                                     * The final two parameters are the number of relaxation
+                                     * steps and the relaxation parameter.
                                      */
 
     template<int dim>
-    MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
-                 const MGMatrix<SparseMatrix<float> >& matrix,
-                 unsigned int steps);
+    MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+                        const MGMatrix<SparseMatrix<number> >& matrix,
+                        function_ptr function,
+                        unsigned int steps,
+                        double omega = 1.);
     
                                     /**
                                      * Implementation of the interface in #MGSmootherBase#.
@@ -145,13 +158,25 @@ class MGSmootherSOR
                                      * and find a way to keep the boundary values.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const;
   private:
                                     /**
                                      * Pointer to the matrices.
                                      */
-    SmartPointer< const MGMatrix< SparseMatrix<float> > > matrix;
+    SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
+                                    /**
+                                     * Pointer to the relaxation function.
+                                     */
+    function_ptr relaxation;
+                                    /**
+                                     * Relaxation parameter.
+                                     */
+    double omega;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    mutable Vector<double> h1, h2;
 };
 
 inline
diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h
new file mode 100644 (file)
index 0000000..f72af10
--- /dev/null
@@ -0,0 +1,34 @@
+// $Id$
+
+template<typename number>
+template<int dim>
+MGSmootherRelaxation<number>::
+MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+                    const MGMatrix<SparseMatrix<number> >& matrix,
+                    function_ptr relaxation,
+                    unsigned int steps,
+                    double omega)
+               :
+               MGSmoother(mg_dof, steps),
+               matrix(&matrix),
+               relaxation(relaxation),
+               omega(omega)
+{}
+
+template<typename number>
+void
+MGSmootherRelaxation<number>::smooth (const unsigned int   level,
+                      Vector<double>       &u,
+                      const Vector<double> &rhs) const
+{
+  h1.reinit(u);
+  h2.reinit(u);
+  for(unsigned i=0;i<get_steps();++i)
+    {
+      (*matrix)[level].residual(h1, u, rhs);
+      ((*matrix)[level].*relaxation)(h2,h1,omega);
+      set_zero_interior_boundary(level,h2);
+      u.add(h2);
+    }
+}
+
index d9587dadfee635d179fb3f52c54777bf3503b440..74a92453c0dc105bf15918c66dd6fd14fab89c7c 100644 (file)
@@ -29,8 +29,10 @@ class MG
 {
   public:
     MG(const MGDoFHandler<dim>&,
-       const MGMatrix<SparseMatrix<float> >&,
-       const MGTransferBase& transfer);
+       ConstraintMatrix& hanging_nodes,
+       const MGMatrix<SparseMatrix<double> >&,
+       const MGTransferBase& transfer,
+       unsigned int minlevel = 0, unsigned int maxlevel = 10000);
     
                                     /** Transfer from dVector to
                                      * MGVector.
@@ -58,17 +60,18 @@ class MG
     void copy_from_mg(Vector<number>& dst) const;
 
                                     /**
-                                     * Residual on a level.
+                                     * Negative #vmult# on a level.
+                                     * @see MGBase.
                                      */
-    virtual void level_residual(unsigned int,
-                               Vector<float> &,
-                               const Vector<float> &,
-                               const Vector<float> &);
+    virtual void level_vmult(unsigned int,
+                               Vector<double> &,
+                               const Vector<double> &,
+                               const Vector<double> &);
 
                                     /**
                                      * Read-only access to level matrices.
                                      */
-    const SparseMatrix<float>& get_matrix(unsigned int level) const;
+    const SparseMatrix<double>& get_matrix(unsigned int level) const;
 
   private:
                                     /**
@@ -82,7 +85,8 @@ class MG
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
+    SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
+    ConstraintMatrix& hanging_nodes;
 };
 
 
@@ -97,8 +101,16 @@ class MG
  */
 class MGTransferPrebuilt : public MGTransferBase 
 {
+    const MGSmoother& smoother;
   public:
-
+                                    /**
+                                     * Preliminary constructor
+                                     */
+    MGTransferPrebuilt(const MGSmoother& smoother) :
+                   smoother(smoother)
+      {}
+    
+    
                                     /**
                                      * Destructor.
                                      */
@@ -124,8 +136,8 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const;
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const;
 
                                     /**
                                      * Restrict a vector from level
@@ -141,8 +153,8 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const;
+                          Vector<double>       &dst,
+                          const Vector<double> &src) const;
 
   private:
 
@@ -162,7 +174,7 @@ class MGTransferPrebuilt : public MGTransferBase
 
 template<int dim>
 inline
-const SparseMatrix<float>&
+const SparseMatrix<double>&
 MG<dim>::get_matrix(unsigned int level) const
 {
   Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
index 9e170e3d5b2553474f24f346d41d49daca682d4d..55ae154eee936b9c28811391b491879530498427 100644 (file)
@@ -1,15 +1,25 @@
 // $Id$
 
 #include <numerics/multigrid.h>
+#include <fstream>
+
+inline int minimum(int a, int b)
+{
+  return ((a<b) ? a : b);
+}
 
 template <int dim>
 MG<dim>::MG(const MGDoFHandler<dim>& dofs,
-           const MGMatrix<SparseMatrix<float> >& matrices,
-           const MGTransferBase& transfer)
+           ConstraintMatrix& hanging_nodes,
+           const MGMatrix<SparseMatrix<double> >& matrices,
+           const MGTransferBase& transfer,
+           unsigned int minl, unsigned int maxl)
                :
-               MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
+               MGBase(transfer, minl,
+                      minimum(dofs.get_tria().n_levels()-1, maxl)),
                dofs(&dofs),
-               matrices(&matrices)
+               matrices(&matrices),
+               hanging_nodes(hanging_nodes)
 {
   for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
@@ -24,19 +34,59 @@ void
 MG<dim>::copy_to_mg(const Vector<number>& src)
 {
   const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  const unsigned int face_dofs = dofs->get_fe().dofs_per_face;
+  
+  d.clear();
+  hanging_nodes.condense(src);
   
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
+  vector<int> mgfaceindex(face_dofs);
+  
+  {
+    ofstream of("CT");
+    DataOut<2> out;
+    out.attach_dof_handler(*dofs);
+    out.add_data_vector(src,"solution","m");
+    out.write_gnuplot(of,1);
+  }
 
-  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 (int level = maxlevel ; level >= (int) minlevel ; --level)
     {
-      dc->get_dof_indices(index);
-      c->get_mg_dof_indices(mgindex);
-      for (unsigned int i=0;i<fe_dofs;++i)
-       d[maxlevel](mgindex[i]) = src(index[i]);
-    } 
+      DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
+      MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
+      for (; c != dofs->end(level) ; ++c, ++dc)
+       {
+         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]);
+
+                                          // Delete values on refinement edge
+         for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+           {
+             MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+             if (face->has_children())
+               {
+                 face->get_mg_dof_indices(mgfaceindex);
+                 for (unsigned int i=0;i<face_dofs;++i)
+                   d[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]);
+    }
+  
+  for (unsigned int i=minlevel; i<= maxlevel; ++i)
+    {
+      char name[10];
+      sprintf(name,"CT%d",i);
+      ofstream of(name);
+      write_gnuplot(*dofs, d[i], i, of);
+    }
 }
 
 template <int dim> template <typename number>
@@ -48,25 +98,42 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
 
+  for (unsigned int i=minlevel; i<= maxlevel; ++i)
+    {
+      char name[10];
+      sprintf(name,"CF%d",i);
+      ofstream of(name);
+      write_gnuplot(*dofs, s[i], i, of);
+    }
+
   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)
+                                           ; 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[maxlevel](mgindex[i]);
-    }
+       dst(index[i]) = s[c->level()](mgindex[i]);
+    } 
+  hanging_nodes.set_zero(dst);
+  {
+    ofstream of("CF");
+    DataOut<2> out;
+    out.attach_dof_handler(*dofs);
+    out.add_data_vector(dst,"solution","m");
+    out.write_gnuplot(of,1);
+  }
 }
 
 template <int dim>
 void
-MG<dim>::level_residual(unsigned int l,
-                       Vector<float> &r,
-                       const Vector<float> &u,
-                       const Vector<float> &b)
+MG<dim>::level_vmult(unsigned int l,
+                       Vector<double> &r,
+                       const Vector<double> &u,
+                       const Vector<double> &)
 {
-  (*matrices)[l].residual(r,u,b);
+  (*matrices)[l].vmult(r,u);
+  r.scale(-1.);
 }
 
 /*
index 538f4428afa579b005bac1d6c4e1502ecc3bb1e4..75b6ee8c1340199e0aa3d802baec0c8cff1e1133 100644 (file)
@@ -809,14 +809,14 @@ class EpsOutputData{
                                      * used for the height information, within
                                      * the list of DoF data vectors.
                                      */
-    unsigned height_vector;
+    unsigned int height_vector;
 
                                     /**
                                      * Number of the vector which is to be
                                      * used for the cell shading values, within
                                      * the list of cell data vectors.
                                      */
-    unsigned cell_vector;
+    unsigned int cell_vector;
     
                                      /**
                                      * Azimuth of the spectators position.
index 6db42e8039943910de355993d7893efca366fe1f..110cca5eb1cdbe4b2348274abc98f3c6292b9118 100644 (file)
@@ -22,8 +22,7 @@
  *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGSmoother
-  :  public MGSmootherBase
+class MGSmoother : public MGSmootherBase
 {
   private:
                                     /**
@@ -82,7 +81,7 @@ class MGSmoother
                                      * function does nothing in this case.
                                      */
     void set_zero_interior_boundary (const unsigned int  level,
-                                    Vector<float>      &u) const;
+                                    Vector<double>      &u) const;
                                     /**
                                      * Modify the number of smoothing steps.
                                      */
@@ -113,31 +112,45 @@ class MGSmoother
 };
 
 /**
- * Implementation of an SOR smoother.
- *
- * This class is an implementation of a smootin algorithm for multigrid. It
- * knows the actual instance of the multigrid matrix and uses the SOR method of the level
- * matrices.
+ * Implementation of a smoother using matrix builtin relaxation methods.
+ * 
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGSmootherSOR
-  :
-  public MGSmoother
+template<typename number>
+class MGSmootherRelaxation : public MGSmoother
 {
   public:
+                                    /**
+                                     * Type of the smoothing
+                                     * function of the matrix.
+                                     */
+    typedef void
+    (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
+                                          typename SparseMatrix<number>::value_type) const;
+    
                                     /**
                                      * Constructor.
                                      * The constructor uses an #MGDoFHandler# to initialize
                                      * data structures for interior level boundary handling
-                                     * in #MGSmoother#. Furthermore, it takes a pointer to the
-                                     * matrices and the number of SOR steps required in each
-                                     * smoothing operation.
+                                     * in #MGSmoother#.
+                                     *
+                                     * Furthermore, it takes a pointer to the
+                                     * level matrices and their smoothing function.
+                                     * This function must perform one relaxation step
+                                     * like #SparseMatrix<number>::SOR_step# does. Do not
+                                     * use the preconditioning methods because they apply
+                                     * a preconditioning matrix to the residual.
+                                     *
+                                     * The final two parameters are the number of relaxation
+                                     * steps and the relaxation parameter.
                                      */
 
     template<int dim>
-    MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
-                 const MGMatrix<SparseMatrix<float> >& matrix,
-                 unsigned int steps);
+    MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+                        const MGMatrix<SparseMatrix<number> >& matrix,
+                        function_ptr function,
+                        unsigned int steps,
+                        double omega = 1.);
     
                                     /**
                                      * Implementation of the interface in #MGSmootherBase#.
@@ -145,13 +158,25 @@ class MGSmootherSOR
                                      * and find a way to keep the boundary values.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const;
   private:
                                     /**
                                      * Pointer to the matrices.
                                      */
-    SmartPointer< const MGMatrix< SparseMatrix<float> > > matrix;
+    SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
+                                    /**
+                                     * Pointer to the relaxation function.
+                                     */
+    function_ptr relaxation;
+                                    /**
+                                     * Relaxation parameter.
+                                     */
+    double omega;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    mutable Vector<double> h1, h2;
 };
 
 inline
diff --git a/deal.II/deal.II/include/numerics/mg_smoother.templates.h b/deal.II/deal.II/include/numerics/mg_smoother.templates.h
new file mode 100644 (file)
index 0000000..f72af10
--- /dev/null
@@ -0,0 +1,34 @@
+// $Id$
+
+template<typename number>
+template<int dim>
+MGSmootherRelaxation<number>::
+MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+                    const MGMatrix<SparseMatrix<number> >& matrix,
+                    function_ptr relaxation,
+                    unsigned int steps,
+                    double omega)
+               :
+               MGSmoother(mg_dof, steps),
+               matrix(&matrix),
+               relaxation(relaxation),
+               omega(omega)
+{}
+
+template<typename number>
+void
+MGSmootherRelaxation<number>::smooth (const unsigned int   level,
+                      Vector<double>       &u,
+                      const Vector<double> &rhs) const
+{
+  h1.reinit(u);
+  h2.reinit(u);
+  for(unsigned i=0;i<get_steps();++i)
+    {
+      (*matrix)[level].residual(h1, u, rhs);
+      ((*matrix)[level].*relaxation)(h2,h1,omega);
+      set_zero_interior_boundary(level,h2);
+      u.add(h2);
+    }
+}
+
index d9587dadfee635d179fb3f52c54777bf3503b440..74a92453c0dc105bf15918c66dd6fd14fab89c7c 100644 (file)
@@ -29,8 +29,10 @@ class MG
 {
   public:
     MG(const MGDoFHandler<dim>&,
-       const MGMatrix<SparseMatrix<float> >&,
-       const MGTransferBase& transfer);
+       ConstraintMatrix& hanging_nodes,
+       const MGMatrix<SparseMatrix<double> >&,
+       const MGTransferBase& transfer,
+       unsigned int minlevel = 0, unsigned int maxlevel = 10000);
     
                                     /** Transfer from dVector to
                                      * MGVector.
@@ -58,17 +60,18 @@ class MG
     void copy_from_mg(Vector<number>& dst) const;
 
                                     /**
-                                     * Residual on a level.
+                                     * Negative #vmult# on a level.
+                                     * @see MGBase.
                                      */
-    virtual void level_residual(unsigned int,
-                               Vector<float> &,
-                               const Vector<float> &,
-                               const Vector<float> &);
+    virtual void level_vmult(unsigned int,
+                               Vector<double> &,
+                               const Vector<double> &,
+                               const Vector<double> &);
 
                                     /**
                                      * Read-only access to level matrices.
                                      */
-    const SparseMatrix<float>& get_matrix(unsigned int level) const;
+    const SparseMatrix<double>& get_matrix(unsigned int level) const;
 
   private:
                                     /**
@@ -82,7 +85,8 @@ class MG
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
+    SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
+    ConstraintMatrix& hanging_nodes;
 };
 
 
@@ -97,8 +101,16 @@ class MG
  */
 class MGTransferPrebuilt : public MGTransferBase 
 {
+    const MGSmoother& smoother;
   public:
-
+                                    /**
+                                     * Preliminary constructor
+                                     */
+    MGTransferPrebuilt(const MGSmoother& smoother) :
+                   smoother(smoother)
+      {}
+    
+    
                                     /**
                                      * Destructor.
                                      */
@@ -124,8 +136,8 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const;
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const;
 
                                     /**
                                      * Restrict a vector from level
@@ -141,8 +153,8 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const;
+                          Vector<double>       &dst,
+                          const Vector<double> &src) const;
 
   private:
 
@@ -162,7 +174,7 @@ class MGTransferPrebuilt : public MGTransferBase
 
 template<int dim>
 inline
-const SparseMatrix<float>&
+const SparseMatrix<double>&
 MG<dim>::get_matrix(unsigned int level) const
 {
   Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
index 9e170e3d5b2553474f24f346d41d49daca682d4d..55ae154eee936b9c28811391b491879530498427 100644 (file)
@@ -1,15 +1,25 @@
 // $Id$
 
 #include <numerics/multigrid.h>
+#include <fstream>
+
+inline int minimum(int a, int b)
+{
+  return ((a<b) ? a : b);
+}
 
 template <int dim>
 MG<dim>::MG(const MGDoFHandler<dim>& dofs,
-           const MGMatrix<SparseMatrix<float> >& matrices,
-           const MGTransferBase& transfer)
+           ConstraintMatrix& hanging_nodes,
+           const MGMatrix<SparseMatrix<double> >& matrices,
+           const MGTransferBase& transfer,
+           unsigned int minl, unsigned int maxl)
                :
-               MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
+               MGBase(transfer, minl,
+                      minimum(dofs.get_tria().n_levels()-1, maxl)),
                dofs(&dofs),
-               matrices(&matrices)
+               matrices(&matrices),
+               hanging_nodes(hanging_nodes)
 {
   for (unsigned int level = minlevel; level<=maxlevel ; ++level)
     {
@@ -24,19 +34,59 @@ void
 MG<dim>::copy_to_mg(const Vector<number>& src)
 {
   const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  const unsigned int face_dofs = dofs->get_fe().dofs_per_face;
+  
+  d.clear();
+  hanging_nodes.condense(src);
   
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
+  vector<int> mgfaceindex(face_dofs);
+  
+  {
+    ofstream of("CT");
+    DataOut<2> out;
+    out.attach_dof_handler(*dofs);
+    out.add_data_vector(src,"solution","m");
+    out.write_gnuplot(of,1);
+  }
 
-  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 (int level = maxlevel ; level >= (int) minlevel ; --level)
     {
-      dc->get_dof_indices(index);
-      c->get_mg_dof_indices(mgindex);
-      for (unsigned int i=0;i<fe_dofs;++i)
-       d[maxlevel](mgindex[i]) = src(index[i]);
-    } 
+      DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
+      MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
+      for (; c != dofs->end(level) ; ++c, ++dc)
+       {
+         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]);
+
+                                          // Delete values on refinement edge
+         for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+           {
+             MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+             if (face->has_children())
+               {
+                 face->get_mg_dof_indices(mgfaceindex);
+                 for (unsigned int i=0;i<face_dofs;++i)
+                   d[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]);
+    }
+  
+  for (unsigned int i=minlevel; i<= maxlevel; ++i)
+    {
+      char name[10];
+      sprintf(name,"CT%d",i);
+      ofstream of(name);
+      write_gnuplot(*dofs, d[i], i, of);
+    }
 }
 
 template <int dim> template <typename number>
@@ -48,25 +98,42 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
 
+  for (unsigned int i=minlevel; i<= maxlevel; ++i)
+    {
+      char name[10];
+      sprintf(name,"CF%d",i);
+      ofstream of(name);
+      write_gnuplot(*dofs, s[i], i, of);
+    }
+
   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)
+                                           ; 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[maxlevel](mgindex[i]);
-    }
+       dst(index[i]) = s[c->level()](mgindex[i]);
+    } 
+  hanging_nodes.set_zero(dst);
+  {
+    ofstream of("CF");
+    DataOut<2> out;
+    out.attach_dof_handler(*dofs);
+    out.add_data_vector(dst,"solution","m");
+    out.write_gnuplot(of,1);
+  }
 }
 
 template <int dim>
 void
-MG<dim>::level_residual(unsigned int l,
-                       Vector<float> &r,
-                       const Vector<float> &u,
-                       const Vector<float> &b)
+MG<dim>::level_vmult(unsigned int l,
+                       Vector<double> &r,
+                       const Vector<double> &u,
+                       const Vector<double> &)
 {
-  (*matrices)[l].residual(r,u,b);
+  (*matrices)[l].vmult(r,u);
+  r.scale(-1.);
 }
 
 /*
index 1af3e7944b3cbbd44477a43f9d50ffbbd81f4ac1..e66b41e8cf92fdc8179dcd712213c6e5c4ee0cf9 100644 (file)
@@ -332,6 +332,7 @@ template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) cons
 template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
                            Vector<number>       &condensed) const;
 template void ConstraintMatrix::condense(Vector<number> &vec) const;
+template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
 template void ConstraintMatrix::distribute(const Vector<number> &condensed,
                                      Vector<number>       &uncondensed) const;
 template void ConstraintMatrix::distribute(Vector<number> &vec) const;
@@ -345,6 +346,9 @@ template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) cons
 template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
                            Vector<number>       &condensed) const;
 template void ConstraintMatrix::condense(Vector<number> &vec) const;
+template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
 template void ConstraintMatrix::distribute(const Vector<number> &condensed,
                                      Vector<number>       &uncondensed) const;
 template void ConstraintMatrix::distribute(Vector<number> &vec) const;
+
+#undef number
index e6340762e8d44d4eb285bd9bfb3952abbd780e4f..9df47ff5949ba4b6e4ce8af46e129b7cc60caba5 100644 (file)
@@ -63,7 +63,7 @@ MGBase::level_mgstep(unsigned level,
   pre_smooth.smooth(level, s[level], d[level]);
                                   // t = d-As
   t.reinit(s[level].size());
-  level_residual(level, t, s[level], d[level]);
+  level_vmult(level, t, s[level], d[level]);
 
                                   // make t rhs of lower level
   transfer->restrict(level, d[level-1], t);
@@ -100,8 +100,8 @@ MGSmootherBase::~MGSmootherBase()
 
 void
 MGSmootherIdentity::smooth (const unsigned int,
-                           Vector<float>       &,
-                           const Vector<float> &) const
+                           Vector<double>       &,
+                           const Vector<double> &) const
 {}
 
 
index 23bea7acf92eb1e14eb0ed1197fe88f1e87a8427..9fad33f9c8cb3c2fd6c8439cadf0175842a6d490 100644 (file)
@@ -92,7 +92,7 @@ MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
 
 void
 MGSmoother::set_zero_interior_boundary (const unsigned int  level,
-                                       Vector<float>      &u) const
+                                       Vector<double>      &u) const
 {
   if (level==0)
     return;
@@ -104,22 +104,7 @@ MGSmoother::set_zero_interior_boundary (const unsigned int  level,
 
 //////////////////////////////////////////////////////////////////////
 
-template<int dim>
-MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
-                            const MGMatrix<SparseMatrix<float> >& matrix,
-                            unsigned int steps)
-               :
-               MGSmoother(mg_dof, steps),
-               matrix(&matrix)
-{}
-
-void
-MGSmootherSOR::smooth (const unsigned int   /*level*/,
-                      Vector<float>       &/*u*/,
-                      const Vector<float> &/*rhs*/) const
-{
-
-}
+#include <numerics/mg_smoother.templates.h>
 
 // explicit instantiations
 // don't do the following instantiation in 1d, since there is a specialized
@@ -128,7 +113,18 @@ MGSmootherSOR::smooth (const unsigned int   /*level*/,
 template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
 #endif
 
-template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_dof,
-                                     const MGMatrix<SparseMatrix<float> >& matrix,
-                                     unsigned int steps);
+template
+MGSmootherRelaxation<float>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+                      const MGMatrix<SparseMatrix<float> >& matrix,
+                      function_ptr relaxation,
+                      unsigned int steps,
+                      double omega);
+template
+MGSmootherRelaxation<double>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+                      const MGMatrix<SparseMatrix<double> >& matrix,
+                      function_ptr relaxation,
+                      unsigned int steps,
+                      double omega);
 
index 429f929615908d12effdc73a45f973ac2543c2ce..2785d302c6da1ed839a30550cfb2177f5ab841af 100644 (file)
@@ -11,6 +11,8 @@
 #include <numerics/mg_smoother.h>
 #include <lac/vector.h>
 
+extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector<double>& v, unsigned int level,
+                          ostream &out, unsigned int accuracy);
 
 
 /* ----------------------------- MGTransferPrebuilt ------------------------ */
@@ -59,7 +61,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
                                       // cell
       prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
                                             mg_dof.n_dofs(level),
-                                            dofs_per_cell);
+// evil hack, must be corrected!!!
+                                            dofs_per_cell+1);
       
       for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
           cell != mg_dof.end(level); ++cell)
@@ -134,8 +137,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
 
 
 void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
-                                    Vector<float>       &dst,
-                                    const Vector<float> &src) const 
+                                    Vector<double>       &dst,
+                                    const Vector<double> &src) const 
 {
   Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
          ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
@@ -146,13 +149,13 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
 
 
 void MGTransferPrebuilt::restrict (const unsigned int   from_level,
-                                  Vector<float>       &dst,
-                                  const Vector<float> &src) const 
+                                  Vector<double>       &dst,
+                                  const Vector<double> &src) const 
 {
   Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
          ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
-
-  prolongation_matrices[from_level-1].Tvmult (dst, src);
+//  smoother.set_zero_interior_boundary(from_level, src);
+  prolongation_matrices[from_level-1].Tvmult_add (dst, src);
 };
 
 
index 23bea7acf92eb1e14eb0ed1197fe88f1e87a8427..9fad33f9c8cb3c2fd6c8439cadf0175842a6d490 100644 (file)
@@ -92,7 +92,7 @@ MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
 
 void
 MGSmoother::set_zero_interior_boundary (const unsigned int  level,
-                                       Vector<float>      &u) const
+                                       Vector<double>      &u) const
 {
   if (level==0)
     return;
@@ -104,22 +104,7 @@ MGSmoother::set_zero_interior_boundary (const unsigned int  level,
 
 //////////////////////////////////////////////////////////////////////
 
-template<int dim>
-MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
-                            const MGMatrix<SparseMatrix<float> >& matrix,
-                            unsigned int steps)
-               :
-               MGSmoother(mg_dof, steps),
-               matrix(&matrix)
-{}
-
-void
-MGSmootherSOR::smooth (const unsigned int   /*level*/,
-                      Vector<float>       &/*u*/,
-                      const Vector<float> &/*rhs*/) const
-{
-
-}
+#include <numerics/mg_smoother.templates.h>
 
 // explicit instantiations
 // don't do the following instantiation in 1d, since there is a specialized
@@ -128,7 +113,18 @@ MGSmootherSOR::smooth (const unsigned int   /*level*/,
 template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
 #endif
 
-template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_dof,
-                                     const MGMatrix<SparseMatrix<float> >& matrix,
-                                     unsigned int steps);
+template
+MGSmootherRelaxation<float>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+                      const MGMatrix<SparseMatrix<float> >& matrix,
+                      function_ptr relaxation,
+                      unsigned int steps,
+                      double omega);
+template
+MGSmootherRelaxation<double>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+                      const MGMatrix<SparseMatrix<double> >& matrix,
+                      function_ptr relaxation,
+                      unsigned int steps,
+                      double omega);
 
index 429f929615908d12effdc73a45f973ac2543c2ce..2785d302c6da1ed839a30550cfb2177f5ab841af 100644 (file)
@@ -11,6 +11,8 @@
 #include <numerics/mg_smoother.h>
 #include <lac/vector.h>
 
+extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector<double>& v, unsigned int level,
+                          ostream &out, unsigned int accuracy);
 
 
 /* ----------------------------- MGTransferPrebuilt ------------------------ */
@@ -59,7 +61,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
                                       // cell
       prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
                                             mg_dof.n_dofs(level),
-                                            dofs_per_cell);
+// evil hack, must be corrected!!!
+                                            dofs_per_cell+1);
       
       for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
           cell != mg_dof.end(level); ++cell)
@@ -134,8 +137,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
 
 
 void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
-                                    Vector<float>       &dst,
-                                    const Vector<float> &src) const 
+                                    Vector<double>       &dst,
+                                    const Vector<double> &src) const 
 {
   Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
          ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
@@ -146,13 +149,13 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
 
 
 void MGTransferPrebuilt::restrict (const unsigned int   from_level,
-                                  Vector<float>       &dst,
-                                  const Vector<float> &src) const 
+                                  Vector<double>       &dst,
+                                  const Vector<double> &src) const 
 {
   Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
          ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
-
-  prolongation_matrices[from_level-1].Tvmult (dst, src);
+//  smoother.set_zero_interior_boundary(from_level, src);
+  prolongation_matrices[from_level-1].Tvmult_add (dst, src);
 };
 
 
index d464b68f083909d0a588ccd91bb1ae620b6da329..8e860139e7c2721d42e495489972d7f93d2b6afc 100644 (file)
@@ -40,8 +40,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<double>& dst,
+                            const Vector<double>& src) const = 0;
 };
 
 
@@ -74,8 +74,8 @@ class MGSmootherBase
                                      * things.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const = 0;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const = 0;
 
 };
 
@@ -93,8 +93,8 @@ class MGSmootherIdentity
                                      * This function does nothing.
                                      */
     virtual void smooth (const unsigned int   level,
-                        Vector<float>       &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double>       &u,
+                        const Vector<double> &rhs) const;
 };
 
 /**
@@ -129,8 +129,8 @@ class MGTransferBase
                                      * level.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const = 0;
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const = 0;
 
                                     /**
                                      * Restrict a vector from level
@@ -146,8 +146,8 @@ class MGTransferBase
                                      * level.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const = 0;
+                          Vector<double>       &dst,
+                          const Vector<double> &src) const = 0;
 };
 
 /**
@@ -175,6 +175,10 @@ class MGVector
                                      * Safe access operator.
                                      */
     const VECTOR& operator[](unsigned int) const;
+                                    /**
+                                     * Delete all entries.
+                                     */
+    void clear();
   private:
                                     /**
                                      * Level of first component.
@@ -245,8 +249,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<double>& dst,
+                            const Vector<double>& src) const;
   private:
                                     /**
                                      * Reference to the solver.
@@ -303,23 +307,23 @@ class MGBase
                                     /**
                                      * Auxiliary vector, defect.
                                      */
-    MGVector<Vector<float> > d;
+    MGVector<Vector<double> > d;
 
                                     /**
                                      * Auxiliary vector, solution.
                                      */
-    MGVector<Vector<float> > s;
+    MGVector<Vector<double> > s;
+                                    /**
+                                     * Prolongation and restriction object.
+                                     */
+    SmartPointer<const MGTransferBase> transfer;
 
   private:
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector<float> t;
+    Vector<double> t;
     
-                                    /**
-                                     * Prolongation and restriction object.
-                                     */
-    SmartPointer<const MGTransferBase> transfer;
                                     /**
                                      * Exception.
                                      */
@@ -344,15 +348,15 @@ class MGBase
                      const MGCoarseGridSolver& cgs);  
 
                                     /**
-                                     * Apply residual operator on all
+                                     * Apply negative #vmult# on all
                                      * cells of a level.
                                      * This is implemented in a
                                      * derived class.
                                      */
-    virtual void level_residual(unsigned int level,
-                            Vector<float>& dst,
-                            const Vector<float>& src,
-                            const Vector<float>& rhs) = 0;  
+    virtual void level_vmult(unsigned int level,
+                            Vector<double>& dst,
+                            const Vector<double>& src,
+                            const Vector<double>& rhs) = 0;  
 
   
   public:
@@ -450,7 +454,7 @@ MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
 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<double>& dst, const Vector<double>& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
@@ -485,6 +489,15 @@ MGVector<VECTOR>::operator[](unsigned int i) const
   return vector<VECTOR>::operator[](i-minlevel);
 }
 
+template<class VECTOR>
+void
+MGVector<VECTOR>::clear()
+{
+  vector<VECTOR>::iterator v;
+  for (v = begin(); v != end(); ++v)
+    v->clear();
+  
+}
 
 
 /* ------------------------------------------------------------------- */
index 198697afb1f69aa7a4d814034de0351707d25b90..4f0acff54a6f99b98ac5766b62bbd884702cad21 100644 (file)
@@ -139,7 +139,8 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                       // and since it makes the backward
                                       // step when applying the decomposition
                                       // significantly faster
-      Assert((global_entry(rowstart_indices[row-1]) !=0), ExcDivideByZero());
+      AssertThrow((global_entry(rowstart_indices[row-1]) !=0),
+                 ExcDivideByZero());
       
       global_entry (rowstart_indices[row-1])
        = 1./global_entry (rowstart_indices[row-1]);
index e6340762e8d44d4eb285bd9bfb3952abbd780e4f..9df47ff5949ba4b6e4ce8af46e129b7cc60caba5 100644 (file)
@@ -63,7 +63,7 @@ MGBase::level_mgstep(unsigned level,
   pre_smooth.smooth(level, s[level], d[level]);
                                   // t = d-As
   t.reinit(s[level].size());
-  level_residual(level, t, s[level], d[level]);
+  level_vmult(level, t, s[level], d[level]);
 
                                   // make t rhs of lower level
   transfer->restrict(level, d[level-1], t);
@@ -100,8 +100,8 @@ MGSmootherBase::~MGSmootherBase()
 
 void
 MGSmootherIdentity::smooth (const unsigned int,
-                           Vector<float>       &,
-                           const Vector<float> &) const
+                           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.