]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Big clean-up (part 2).
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Nov 1999 15:32:34 +0000 (15:32 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Nov 1999 15:32:34 +0000 (15:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@1903 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/include/numerics/multigrid.h
deal.II/deal.II/include/numerics/multigrid.templates.h
deal.II/deal.II/source/multigrid/multigrid.cc
deal.II/deal.II/source/numerics/multigrid.cc

index 6086ecbdb4e407d5813af66b7096bb3613720e79..1c0b779a78d5a3ea9dc1dfa76da0daba42f50fd7 100644 (file)
@@ -14,6 +14,8 @@
 
 #include <vector>
 
+
+
 /**
  * Implementation of multigrid with matrices.
  * While #MGBase# was only an abstract framework for the v-cycle,
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
 template<int dim>
-class MG
-  :
-  public MGBase
+class MG : public MGBase
 {
   public:
-    MG(const MGDoFHandler<dim>&,
-       ConstraintMatrix& hanging_nodes,
-       const MGMatrix<SparseMatrix<double> >&,
-       const MGTransferBase& transfer,
-       unsigned int minlevel = 0, unsigned int maxlevel = 10000);
-    
-                                    /** Transfer from dVector to
-                                     * MGVector.
+                                    /**
+                                     * Constructor. Take the object
+                                     * describing the multilevel
+                                     * hierarchy of degrees of
+                                     * freedom, an object describing
+                                     * constraints to degrees of
+                                     * freedom on the global grid
+                                     * (for example for hanging
+                                     * nodes), a set of matrices
+                                     * which describe the equation
+                                     * discretized on each level, and
+                                     * an object mediating the
+                                     * transfer of solutions between
+                                     * different levels.
                                      *
-                                     * This function copies data from a
-                                     * dVector, that is, data on the
-                                     * locally finest level, into the
-                                     * corresponding levels of an
-                                     * MGVector.
+                                     * This function already
+                                     * initializes the vectors which
+                                     * will be used later on in the
+                                     * course of the
+                                     * computations. You should
+                                     * therefore create objects of
+                                     * this type as late as possible.
+                                     */
+//TODO: minlevel, maxlevel?
+    MG(const MGDoFHandler<dim>               &mg_dof_handler,
+       ConstraintMatrix                      &constraints,
+       const MGMatrix<SparseMatrix<double> > &level_matrices,
+       const MGTransferBase                  &transfer,
+       const unsigned int                     minlevel = 0,
+       const unsigned int                     maxlevel = 1000000);
+    
+                                    /**
+                                     * Transfer from a vector on the
+                                     * global grid to vectors defined
+                                     * on each of the levels
+                                     * separately, i.a. an #MGVector#.
                                      */
     template<typename number>
-    void copy_to_mg(const Vector<number>& src);
+    void copy_to_mg (const Vector<number> &src);
 
                                     /**
-                                     * Transfer from multi-grid vector to
+                                     * Transfer from multi-level vector to
                                      * normal vector.
                                      *
-                                     * Copies data from active portions
-                                     * of an MGVector into the
-                                     * respective positions of a
-                                     * Vector<double>. All other entries of
-                                     * #src# are zero.
+                                     * Copies data from active
+                                     * portions of an MGVector into
+                                     * the respective positions of a
+                                     * #Vector<double>#. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom are set to zero.
                                      */
     template<typename number>
-    void copy_from_mg(Vector<number>& dst) const;
+    void copy_from_mg (Vector<number> &dst) const;
 
                                     /**
                                      * Negative #vmult# on a level.
+//TODO: what does this function do?                                  
                                      * @see MGBase.
                                      */
-    virtual void level_vmult(unsigned int,
-                               Vector<double> &,
-                               const Vector<double> &,
-                               const Vector<double> &);
+    virtual void level_vmult (const unsigned int    level,
+                             Vector<double>       &dest,
+                             const Vector<double> &src,
+                             const Vector<double> &rhs);
 
                                     /**
                                      * Read-only access to level matrices.
                                      */
-    const SparseMatrix<double>& get_matrix(unsigned int level) const;
+    const SparseMatrix<double>& get_matrix (const unsigned int level) const;
 
   private:
                                     /**
                                      * Associated #MGDoFHandler#.
                                      */
-    SmartPointer<const MGDoFHandler<dim> > dofs;
+    SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
 
                                     /**
                                      * Matrices for each level.
@@ -85,11 +109,19 @@ class MG
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
-    ConstraintMatrix& hanging_nodes;
+    SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
+
+                                    /**
+                                     * Pointer to the object holding
+                                     * constraints.
+                                     */
+    SmartPointer<const ConstraintMatrix>                 constraints;
 };
 
 
+
+
+
 /**
  * Implementation of the #MGTransferBase# interface for which the transfer
  * operations are prebuilt upon construction of the object of this class as
@@ -102,13 +134,6 @@ class MG
 class MGTransferPrebuilt : public MGTransferBase 
 {
   public:
-                                    /**
-                                     * Preliminary constructor
-                                     */
-    MGTransferPrebuilt()
-      {}
-    
-    
                                     /**
                                      * Destructor.
                                      */
@@ -123,7 +148,10 @@ class MGTransferPrebuilt : public MGTransferBase
 
                                     /**
                                      * Prolongate a vector from level
-                                     * #to_level-1# to level #to_level#.
+                                     * #to_level-1# to level
+                                     * #to_level#. The previous
+                                     * content of #dst# is
+                                     * overwritten.
                                      *
                                      * #src# is assumed to be a vector with
                                      * as many elements as there are degrees
@@ -133,14 +161,24 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * are degrees of freedom on the finer
                                      * level.
                                      */
-    virtual void prolongate (const unsigned int   to_level,
+    virtual void prolongate (const unsigned int    to_level,
                             Vector<double>       &dst,
                             const Vector<double> &src) const;
 
                                     /**
                                      * Restrict a vector from level
                                      * #from_level# to level
-                                     * #from_level-1#.
+                                     * #from_level-1# and add this
+                                     * restriction to
+                                     * #dst#. Obviously, if the
+                                     * refined region on level
+                                     * #from_level# is smaller than
+                                     * that on level #from_level-1#,
+                                     * some degrees of freedom in
+                                     * #dst# are not covered and will
+                                     * not be altered. For the other
+                                     * degress of freedom, the result
+                                     * of the restriction is added.
                                      *
                                      * #src# is assumed to be a vector with
                                      * as many elements as there are degrees
@@ -150,9 +188,9 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * are degrees of freedom on the coarser
                                      * level.
                                      */
-    virtual void restrict (const unsigned int   from_level,
-                          Vector<double>       &dst,
-                          const Vector<double> &src) const;
+    virtual void restrict_and_add (const unsigned int    from_level,
+                                  Vector<double>       &dst,
+                                  const Vector<double> &src) const;
 
   private:
 
@@ -170,14 +208,20 @@ class MGTransferPrebuilt : public MGTransferBase
 };
 
 
+
+
+/* --------------------------- inline functions --------------------- */
+
+
 template<int dim>
 inline
 const SparseMatrix<double>&
 MG<dim>::get_matrix(unsigned int level) const
 {
-  Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
+  Assert((level>=minlevel) && (level<maxlevel),
+        ExcIndexRange(level, minlevel, maxlevel));
   
-  return (*matrices)[level];
+  return (*level_matrices)[level];
 }
 
 
index 5d01b9b16ae90a5e37e2ad227c2fb78dea2aa880..f63dc58d73922b3eecfa1bb80d687631ab9df581 100644 (file)
@@ -1,5 +1,6 @@
 // $Id$
 
+#include <grid/dof_constraints.h>
 #include <numerics/multigrid.h>
 #include <algorithm>
 #include <fstream>
 
 
 
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs,
-           ConstraintMatrix& hanging_nodes,
-           const MGMatrix<SparseMatrix<double> >& matrices,
-           const MGTransferBase& transfer,
-           unsigned int minl, unsigned int maxl)
-               :
-               MGBase(transfer, minl,
-                      min(dofs.get_tria().n_levels()-1, maxl)),
-               dofs(&dofs),
-               matrices(&matrices),
-               hanging_nodes(hanging_nodes)
-{
-  for (unsigned int level = minlevel; level<=maxlevel ; ++level)
-    {
-      solution[level].reinit(matrices[level].m());
-      defect[level].reinit(matrices[level].m());
-    };
-};
-
-
-
 template <int dim>
 template <typename number>
 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;
+  const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
+  const unsigned int face_dofs = mg_dof_handler->get_fe().dofs_per_face;
 
                                   // set the elements of the vectors
                                   // on all levels to zero
   defect.clear();
-//  hanging_nodes.condense(src);
+//  constraints->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");
-//     out.write_gnuplot(of,1);
-//   }
 
-  for (int level = maxlevel; level >= static_cast<int>(minlevel); --level)
+  for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
     {
-      DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
-      MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
-
+      DoFHandler<dim>::active_cell_iterator
+       dc = mg_dof_handler->DoFHandler<dim>::begin_active(level);
+      MGDoFHandler<dim>::active_cell_iterator
+       level_cell = mg_dof_handler->begin_active(level);
+      const MGDoFHandler<dim>::active_cell_iterator
+       level_end  = mg_dof_handler->end_active(level);
+       
 //TODO: Document on difference between #end# and #end_level#
 
-      for (; c != dofs->end_active(level) ; ++c, ++dc)
+      for (; level_cell!=level_end; ++level_cell, ++dc)
        {
          dc->get_dof_indices(index);
-         c->get_mg_dof_indices(mgindex);
-         for (unsigned int i=0;i<fe_dofs;++i)
-           defect[c->level()](mgindex[i]) = src(index[i]);
+         level_cell->get_mg_dof_indices (mgindex);
+         for (unsigned int i=0; i<fe_dofs; ++i)
+           defect[level](mgindex[i]) = src(index[i]);
 
+//TODO: what happens here?
                                           // Delete values on refinement edge
-         for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+         for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
            {
-             MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+             const MGDoFHandler<dim>::face_iterator face = level_cell->face(face_n);
              if (face->has_children())
                {
                  face->get_mg_dof_indices(mgfaceindex);
-                 for (unsigned int i=0;i<face_dofs;++i)
-                   defect[c->level()](mgfaceindex[i]) = 0.;
-               }
-           }
-       }
-//      if (level > (int) minlevel)
-//        transfer->restrict(level, d[level-1], d[level]);
-      if (level < (int) maxlevel)
-       transfer->restrict(level+1, defect[level], defect[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);
-//     }
-}
+                 for (unsigned int i=0; i<face_dofs; ++i)
+                   defect[level](mgfaceindex[i]) = 0.;
+               };
+           };
+       };
+
+                                      // for that part of the level
+                                      // which is further refined:
+                                      // get the defect by
+                                      // restriction of the defect on
+                                      // one level higher
+      if (static_cast<unsigned int>(level) < maxlevel)
+       transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+    };
+};
+
+
 
 template <int dim>
 template <typename number>
 void
 MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
-  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
 
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
@@ -113,16 +87,16 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //       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)
+  DoFHandler<dim>::active_cell_iterator dc = mg_dof_handler->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = mg_dof_handler->begin_active();
+       c != mg_dof_handler->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]) = solution[c->level()](mgindex[i]);
     } 
-  hanging_nodes.set_zero(dst);
+  constraints->set_zero(dst);
 //   {
 //     ofstream of("CF");
 //     DataOut<2> out;
@@ -132,25 +106,3 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //   }
 }
 
-
-
-template <int dim>
-void
-MG<dim>::level_vmult(unsigned int l,
-                       Vector<double> &r,
-                       const Vector<double> &u,
-                       const Vector<double> &)
-{
-  (*matrices)[l].vmult(r,u);
-  r.scale(-1.);
-}
-
-/*
-template <int dim>
-MG::
-template <int dim>
-MG::
-template <int dim>
-MG::
-
-*/
index 6086ecbdb4e407d5813af66b7096bb3613720e79..1c0b779a78d5a3ea9dc1dfa76da0daba42f50fd7 100644 (file)
@@ -14,6 +14,8 @@
 
 #include <vector>
 
+
+
 /**
  * Implementation of multigrid with matrices.
  * While #MGBase# was only an abstract framework for the v-cycle,
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
 template<int dim>
-class MG
-  :
-  public MGBase
+class MG : public MGBase
 {
   public:
-    MG(const MGDoFHandler<dim>&,
-       ConstraintMatrix& hanging_nodes,
-       const MGMatrix<SparseMatrix<double> >&,
-       const MGTransferBase& transfer,
-       unsigned int minlevel = 0, unsigned int maxlevel = 10000);
-    
-                                    /** Transfer from dVector to
-                                     * MGVector.
+                                    /**
+                                     * Constructor. Take the object
+                                     * describing the multilevel
+                                     * hierarchy of degrees of
+                                     * freedom, an object describing
+                                     * constraints to degrees of
+                                     * freedom on the global grid
+                                     * (for example for hanging
+                                     * nodes), a set of matrices
+                                     * which describe the equation
+                                     * discretized on each level, and
+                                     * an object mediating the
+                                     * transfer of solutions between
+                                     * different levels.
                                      *
-                                     * This function copies data from a
-                                     * dVector, that is, data on the
-                                     * locally finest level, into the
-                                     * corresponding levels of an
-                                     * MGVector.
+                                     * This function already
+                                     * initializes the vectors which
+                                     * will be used later on in the
+                                     * course of the
+                                     * computations. You should
+                                     * therefore create objects of
+                                     * this type as late as possible.
+                                     */
+//TODO: minlevel, maxlevel?
+    MG(const MGDoFHandler<dim>               &mg_dof_handler,
+       ConstraintMatrix                      &constraints,
+       const MGMatrix<SparseMatrix<double> > &level_matrices,
+       const MGTransferBase                  &transfer,
+       const unsigned int                     minlevel = 0,
+       const unsigned int                     maxlevel = 1000000);
+    
+                                    /**
+                                     * Transfer from a vector on the
+                                     * global grid to vectors defined
+                                     * on each of the levels
+                                     * separately, i.a. an #MGVector#.
                                      */
     template<typename number>
-    void copy_to_mg(const Vector<number>& src);
+    void copy_to_mg (const Vector<number> &src);
 
                                     /**
-                                     * Transfer from multi-grid vector to
+                                     * Transfer from multi-level vector to
                                      * normal vector.
                                      *
-                                     * Copies data from active portions
-                                     * of an MGVector into the
-                                     * respective positions of a
-                                     * Vector<double>. All other entries of
-                                     * #src# are zero.
+                                     * Copies data from active
+                                     * portions of an MGVector into
+                                     * the respective positions of a
+                                     * #Vector<double>#. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom are set to zero.
                                      */
     template<typename number>
-    void copy_from_mg(Vector<number>& dst) const;
+    void copy_from_mg (Vector<number> &dst) const;
 
                                     /**
                                      * Negative #vmult# on a level.
+//TODO: what does this function do?                                  
                                      * @see MGBase.
                                      */
-    virtual void level_vmult(unsigned int,
-                               Vector<double> &,
-                               const Vector<double> &,
-                               const Vector<double> &);
+    virtual void level_vmult (const unsigned int    level,
+                             Vector<double>       &dest,
+                             const Vector<double> &src,
+                             const Vector<double> &rhs);
 
                                     /**
                                      * Read-only access to level matrices.
                                      */
-    const SparseMatrix<double>& get_matrix(unsigned int level) const;
+    const SparseMatrix<double>& get_matrix (const unsigned int level) const;
 
   private:
                                     /**
                                      * Associated #MGDoFHandler#.
                                      */
-    SmartPointer<const MGDoFHandler<dim> > dofs;
+    SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
 
                                     /**
                                      * Matrices for each level.
@@ -85,11 +109,19 @@ class MG
                                      * the constructor of #MG# and can
                                      * be accessed for assembling.
                                      */
-    SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
-    ConstraintMatrix& hanging_nodes;
+    SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
+
+                                    /**
+                                     * Pointer to the object holding
+                                     * constraints.
+                                     */
+    SmartPointer<const ConstraintMatrix>                 constraints;
 };
 
 
+
+
+
 /**
  * Implementation of the #MGTransferBase# interface for which the transfer
  * operations are prebuilt upon construction of the object of this class as
@@ -102,13 +134,6 @@ class MG
 class MGTransferPrebuilt : public MGTransferBase 
 {
   public:
-                                    /**
-                                     * Preliminary constructor
-                                     */
-    MGTransferPrebuilt()
-      {}
-    
-    
                                     /**
                                      * Destructor.
                                      */
@@ -123,7 +148,10 @@ class MGTransferPrebuilt : public MGTransferBase
 
                                     /**
                                      * Prolongate a vector from level
-                                     * #to_level-1# to level #to_level#.
+                                     * #to_level-1# to level
+                                     * #to_level#. The previous
+                                     * content of #dst# is
+                                     * overwritten.
                                      *
                                      * #src# is assumed to be a vector with
                                      * as many elements as there are degrees
@@ -133,14 +161,24 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * are degrees of freedom on the finer
                                      * level.
                                      */
-    virtual void prolongate (const unsigned int   to_level,
+    virtual void prolongate (const unsigned int    to_level,
                             Vector<double>       &dst,
                             const Vector<double> &src) const;
 
                                     /**
                                      * Restrict a vector from level
                                      * #from_level# to level
-                                     * #from_level-1#.
+                                     * #from_level-1# and add this
+                                     * restriction to
+                                     * #dst#. Obviously, if the
+                                     * refined region on level
+                                     * #from_level# is smaller than
+                                     * that on level #from_level-1#,
+                                     * some degrees of freedom in
+                                     * #dst# are not covered and will
+                                     * not be altered. For the other
+                                     * degress of freedom, the result
+                                     * of the restriction is added.
                                      *
                                      * #src# is assumed to be a vector with
                                      * as many elements as there are degrees
@@ -150,9 +188,9 @@ class MGTransferPrebuilt : public MGTransferBase
                                      * are degrees of freedom on the coarser
                                      * level.
                                      */
-    virtual void restrict (const unsigned int   from_level,
-                          Vector<double>       &dst,
-                          const Vector<double> &src) const;
+    virtual void restrict_and_add (const unsigned int    from_level,
+                                  Vector<double>       &dst,
+                                  const Vector<double> &src) const;
 
   private:
 
@@ -170,14 +208,20 @@ class MGTransferPrebuilt : public MGTransferBase
 };
 
 
+
+
+/* --------------------------- inline functions --------------------- */
+
+
 template<int dim>
 inline
 const SparseMatrix<double>&
 MG<dim>::get_matrix(unsigned int level) const
 {
-  Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
+  Assert((level>=minlevel) && (level<maxlevel),
+        ExcIndexRange(level, minlevel, maxlevel));
   
-  return (*matrices)[level];
+  return (*level_matrices)[level];
 }
 
 
index 5d01b9b16ae90a5e37e2ad227c2fb78dea2aa880..f63dc58d73922b3eecfa1bb80d687631ab9df581 100644 (file)
@@ -1,5 +1,6 @@
 // $Id$
 
+#include <grid/dof_constraints.h>
 #include <numerics/multigrid.h>
 #include <algorithm>
 #include <fstream>
 
 
 
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs,
-           ConstraintMatrix& hanging_nodes,
-           const MGMatrix<SparseMatrix<double> >& matrices,
-           const MGTransferBase& transfer,
-           unsigned int minl, unsigned int maxl)
-               :
-               MGBase(transfer, minl,
-                      min(dofs.get_tria().n_levels()-1, maxl)),
-               dofs(&dofs),
-               matrices(&matrices),
-               hanging_nodes(hanging_nodes)
-{
-  for (unsigned int level = minlevel; level<=maxlevel ; ++level)
-    {
-      solution[level].reinit(matrices[level].m());
-      defect[level].reinit(matrices[level].m());
-    };
-};
-
-
-
 template <int dim>
 template <typename number>
 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;
+  const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
+  const unsigned int face_dofs = mg_dof_handler->get_fe().dofs_per_face;
 
                                   // set the elements of the vectors
                                   // on all levels to zero
   defect.clear();
-//  hanging_nodes.condense(src);
+//  constraints->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");
-//     out.write_gnuplot(of,1);
-//   }
 
-  for (int level = maxlevel; level >= static_cast<int>(minlevel); --level)
+  for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
     {
-      DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
-      MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
-
+      DoFHandler<dim>::active_cell_iterator
+       dc = mg_dof_handler->DoFHandler<dim>::begin_active(level);
+      MGDoFHandler<dim>::active_cell_iterator
+       level_cell = mg_dof_handler->begin_active(level);
+      const MGDoFHandler<dim>::active_cell_iterator
+       level_end  = mg_dof_handler->end_active(level);
+       
 //TODO: Document on difference between #end# and #end_level#
 
-      for (; c != dofs->end_active(level) ; ++c, ++dc)
+      for (; level_cell!=level_end; ++level_cell, ++dc)
        {
          dc->get_dof_indices(index);
-         c->get_mg_dof_indices(mgindex);
-         for (unsigned int i=0;i<fe_dofs;++i)
-           defect[c->level()](mgindex[i]) = src(index[i]);
+         level_cell->get_mg_dof_indices (mgindex);
+         for (unsigned int i=0; i<fe_dofs; ++i)
+           defect[level](mgindex[i]) = src(index[i]);
 
+//TODO: what happens here?
                                           // Delete values on refinement edge
-         for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+         for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
            {
-             MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+             const MGDoFHandler<dim>::face_iterator face = level_cell->face(face_n);
              if (face->has_children())
                {
                  face->get_mg_dof_indices(mgfaceindex);
-                 for (unsigned int i=0;i<face_dofs;++i)
-                   defect[c->level()](mgfaceindex[i]) = 0.;
-               }
-           }
-       }
-//      if (level > (int) minlevel)
-//        transfer->restrict(level, d[level-1], d[level]);
-      if (level < (int) maxlevel)
-       transfer->restrict(level+1, defect[level], defect[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);
-//     }
-}
+                 for (unsigned int i=0; i<face_dofs; ++i)
+                   defect[level](mgfaceindex[i]) = 0.;
+               };
+           };
+       };
+
+                                      // for that part of the level
+                                      // which is further refined:
+                                      // get the defect by
+                                      // restriction of the defect on
+                                      // one level higher
+      if (static_cast<unsigned int>(level) < maxlevel)
+       transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+    };
+};
+
+
 
 template <int dim>
 template <typename number>
 void
 MG<dim>::copy_from_mg(Vector<number>& dst) const
 {
-  const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+  const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
 
   vector<int> index(fe_dofs);
   vector<int> mgindex(fe_dofs);
@@ -113,16 +87,16 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //       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)
+  DoFHandler<dim>::active_cell_iterator dc = mg_dof_handler->DoFHandler<dim>::begin_active();
+  for (MGDoFHandler<dim>::active_cell_iterator c = mg_dof_handler->begin_active();
+       c != mg_dof_handler->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]) = solution[c->level()](mgindex[i]);
     } 
-  hanging_nodes.set_zero(dst);
+  constraints->set_zero(dst);
 //   {
 //     ofstream of("CF");
 //     DataOut<2> out;
@@ -132,25 +106,3 @@ MG<dim>::copy_from_mg(Vector<number>& dst) const
 //   }
 }
 
-
-
-template <int dim>
-void
-MG<dim>::level_vmult(unsigned int l,
-                       Vector<double> &r,
-                       const Vector<double> &u,
-                       const Vector<double> &)
-{
-  (*matrices)[l].vmult(r,u);
-  r.scale(-1.);
-}
-
-/*
-template <int dim>
-MG::
-template <int dim>
-MG::
-template <int dim>
-MG::
-
-*/
index 4b5c8b794670737b79fff8b85ecd22b2fdf64ad3..a593df249eb44a912c1bee9ce5580afa800d84cb 100644 (file)
 #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);
+
+
+/* --------------------------------- MG ------------------------------------ */
+
+template <int dim>
+MG<dim>::MG(const MGDoFHandler<dim>               &mg_dof_handler,
+           ConstraintMatrix                      &constraints,
+           const MGMatrix<SparseMatrix<double> > &level_matrices,
+           const MGTransferBase                  &transfer,
+           unsigned int minl, unsigned int maxl)
+               :
+               MGBase(transfer, minl,
+                      min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
+               mg_dof_handler(&mg_dof_handler),
+               level_matrices(&level_matrices),
+               constraints(&constraints)
+{
+                                  // initialize the objects with
+                                  // their correct size
+  for (unsigned int level=minlevel; level<=maxlevel ; ++level)
+    {
+      solution[level].reinit(level_matrices[level].m());
+      defect[level].reinit(level_matrices[level].m());
+    };
+};
+
+
+
+template <int dim>
+void
+MG<dim>::level_vmult(const unsigned int    level,
+                    Vector<double>       &result,
+                    const Vector<double> &u,
+                    const Vector<double> &/* rhs */)
+{
+  (*level_matrices)[level].vmult(result,u);
+  result.scale(-1.);
+}
+
+
+
+
 
 
 /* ----------------------------- MGTransferPrebuilt ------------------------ */
@@ -61,7 +101,7 @@ 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),
-// evil hack, must be corrected!!!
+//TODO: evil hack, must be corrected!!!
                                             dofs_per_cell+1);
       
       for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
@@ -148,9 +188,9 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
 
 
 
-void MGTransferPrebuilt::restrict (const unsigned int   from_level,
-                                  Vector<double>       &dst,
-                                  const Vector<double> &src) const 
+void MGTransferPrebuilt::restrict_and_add (const unsigned int   from_level,
+                                          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));
@@ -160,7 +200,7 @@ void MGTransferPrebuilt::restrict (const unsigned int   from_level,
 
 
 // explicit instatations
+template class MG<deal_II_dimension>;
 template
 void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
 
-template MG<deal_II_dimension>;
index 4b5c8b794670737b79fff8b85ecd22b2fdf64ad3..a593df249eb44a912c1bee9ce5580afa800d84cb 100644 (file)
 #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);
+
+
+/* --------------------------------- MG ------------------------------------ */
+
+template <int dim>
+MG<dim>::MG(const MGDoFHandler<dim>               &mg_dof_handler,
+           ConstraintMatrix                      &constraints,
+           const MGMatrix<SparseMatrix<double> > &level_matrices,
+           const MGTransferBase                  &transfer,
+           unsigned int minl, unsigned int maxl)
+               :
+               MGBase(transfer, minl,
+                      min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
+               mg_dof_handler(&mg_dof_handler),
+               level_matrices(&level_matrices),
+               constraints(&constraints)
+{
+                                  // initialize the objects with
+                                  // their correct size
+  for (unsigned int level=minlevel; level<=maxlevel ; ++level)
+    {
+      solution[level].reinit(level_matrices[level].m());
+      defect[level].reinit(level_matrices[level].m());
+    };
+};
+
+
+
+template <int dim>
+void
+MG<dim>::level_vmult(const unsigned int    level,
+                    Vector<double>       &result,
+                    const Vector<double> &u,
+                    const Vector<double> &/* rhs */)
+{
+  (*level_matrices)[level].vmult(result,u);
+  result.scale(-1.);
+}
+
+
+
+
 
 
 /* ----------------------------- MGTransferPrebuilt ------------------------ */
@@ -61,7 +101,7 @@ 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),
-// evil hack, must be corrected!!!
+//TODO: evil hack, must be corrected!!!
                                             dofs_per_cell+1);
       
       for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
@@ -148,9 +188,9 @@ void MGTransferPrebuilt::prolongate (const unsigned int   to_level,
 
 
 
-void MGTransferPrebuilt::restrict (const unsigned int   from_level,
-                                  Vector<double>       &dst,
-                                  const Vector<double> &src) const 
+void MGTransferPrebuilt::restrict_and_add (const unsigned int   from_level,
+                                          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));
@@ -160,7 +200,7 @@ void MGTransferPrebuilt::restrict (const unsigned int   from_level,
 
 
 // explicit instatations
+template class MG<deal_II_dimension>;
 template
 void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
 
-template MG<deal_II_dimension>;

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.