]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new multigrid
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 3 Jan 2010 15:52:10 +0000 (15:52 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 3 Jan 2010 15:52:10 +0000 (15:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@20275 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/deal.II/include/multigrid/mg_smoother.h
deal.II/deal.II/include/multigrid/mg_tools.h
deal.II/deal.II/include/multigrid/mg_transfer.h
deal.II/deal.II/include/multigrid/mg_transfer.templates.h
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/source/multigrid/mg_dof_tools.cc
deal.II/deal.II/source/multigrid/mg_smoother.cc [deleted file]
deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc
deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
deal.II/examples/step-39/step-39.cc

index ff47759beff803d17604ab6d486e331c1e211fb1..ef498eed7ad39535b9e121b0e262a712e0e2f311 100644 (file)
@@ -166,115 +166,6 @@ class MGSmootherIdentity : public MGSmootherBase<VECTOR>
 };
 
 
-/**
- * Base class for multigrid smoothers. It implements some
- * functionality for setting the values of vectors at interior
- * boundaries (i.e. boundaries between differing levels of the
- * triangulation) to zero, by building a list of these degrees of
- * freedom's indices at construction time.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002
- */
-template <class VECTOR>
-class MGSmootherContinuous : public MGSmootherBase<VECTOR>
-{
-  private:
-                                    /**
-                                     * Default constructor. Made private
-                                     * to prevent it being called, which
-                                     * is necessary since this could
-                                     * circumvent the set-up of the list
-                                     * if interior boundary dofs.
-                                     */
-    MGSmootherContinuous ();
-    
-  public:
-
-                                    /**
-                                     * Constructor. This one collects
-                                     * the indices of the degrees of freedom
-                                     * on the interior boundaries between
-                                     * the different levels, which are
-                                     * needed by the function
-                                     * @p set_zero_interior_boundaries.
-                                     *
-                                     * Since this function is
-                                     * implemented a bit different in
-                                     * 1d (there are no faces of
-                                     * cells, just vertices), there
-                                     * are actually two sets of
-                                     * constructors, namely this one
-                                     * for 1d and the following one
-                                     * for all other
-                                     * dimensions. Really amusing
-                                     * about this text is, that there
-                                     * is no 1d implementation.
-                                     */
-    MGSmootherContinuous (const MGDoFHandler<1,1> &mg_dof,
-                         const unsigned int     steps);
-
-                                    /**
-                                     * Constructor. This one collects
-                                     * the indices of the degrees of
-                                     * freedom on the interior
-                                     * boundaries between the
-                                     * different levels, which are
-                                     * needed by the function
-                                     * @p set_zero_interior_boundaries.
-                                     *
-                                     * The parameter steps indicates
-                                     * the number of smoothing steps
-                                     * to be executed by @p smooth.
-                                     */
-    template <int dim, int spacedim>
-    MGSmootherContinuous (const MGDoFHandler<dim,spacedim> &mg_dof,
-                         const unsigned int       steps);    
-
-                                    /**
-                                     * Reset the values of the degrees of
-                                     * freedom on interior boundaries between
-                                     * different levels to zero in the given
-                                     * data vector @p u.
-                                     *
-                                     * Since the coarsest level (<tt>level==0</tt>)
-                                     * has no interior boundaries, this
-                                     * function does nothing in this case.
-                                     */
-    void set_zero_interior_boundary (const unsigned int level,
-                                    VECTOR&            u) const;
-    
-                                    /**
-                                     * Modify the number of smoothing steps.
-                                     */
-    void set_steps (const unsigned int steps);
-    
-                                    /**
-                                     * How many steps should be used?
-                                     */
-    unsigned int get_steps() const;
-
-  private:
-                                    /**
-                                     * Number of smoothing steps.
-                                     */
-    unsigned int steps;
-    
-                                    /**
-                                     * For each level, we store a list of
-                                     * degree of freedom indices which are
-                                     * located on interior boundaries between
-                                     * differing levels of the triangulation.
-                                     * Since the coarsest level has no
-                                     * interior boundary dofs, the first
-                                     * entry refers to the second level.
-                                     *
-                                     * These arrays are set by the constructor.
-                                     * The entries for each level are sorted ascendingly.
-                                     */
-    std::vector<std::vector<unsigned int> > interior_boundary_dofs;
-};
-
-
 /**
  * Smoother using relaxation classes.
  *
@@ -568,24 +459,6 @@ MGSmootherIdentity<VECTOR>::smooth (
   const VECTOR&) const
 {}
 
-/*----------------------------------------------------------------------*/
-
-template <class VECTOR>
-inline
-void
-MGSmootherContinuous<VECTOR>::set_steps (const unsigned int i)
-{
-  steps = i;
-}
-
-template <class VECTOR>
-inline
-unsigned int
-MGSmootherContinuous<VECTOR>::get_steps() const
-{
-  return steps;
-}
-
 //---------------------------------------------------------------------------
 
 template <class VECTOR>
index e750f812333283f783a627a433216a1638ee11ab..046decf5daef9dfd946531e8911ba2ea5a2ff63b 100644 (file)
@@ -126,7 +126,6 @@ class MGTools
     make_flux_sparsity_pattern_edge (const MGDoFHandler<dim,spacedim> &dof_handler,
                                     SparsityPattern         &sparsity,
                                     const unsigned int       level);
-
                                     /**
                                      * This function does the same as
                                      * the other with the same name,
@@ -243,12 +242,15 @@ class MGTools
       const typename FunctionMap<dim>::type& function_map,
       std::vector<std::set<unsigned int> >& boundary_indices,
       const std::vector<bool>& component_mask = std::vector<bool>());
+                                      /**
+                                       */
 
     template <typename number>
     static void apply_boundary_values (
       const std::set<unsigned int> &boundary_dofs,
       SparseMatrix<number>& matrix,
-      const bool preserve_symmetry);
+      const bool preserve_symmetry,
+      const bool ignore_zeros = false);
     
     template <typename number>
     static void apply_boundary_values (
index 871d3646033cd42c1b8398b51c49325648918cf7..600af6db4148c052456c573ba181e8584eb7009a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -16,6 +16,7 @@
 #include <base/config.h>
 
 #include <lac/block_vector.h>
+#include <lac/constraint_matrix.h>
 #ifdef DEAL_PREFER_MATRIX_EZ
 #  include <lac/sparse_matrix_ez.h>
 #  include <lac/block_sparse_matrix_ez.h>
@@ -63,11 +64,23 @@ template <class VECTOR>
 class MGTransferPrebuilt : public MGTransferBase<VECTOR> 
 {
   public:
+                                    /**
+                                     * Constructor without constraint
+                                     * matrices. Use this constructor
+                                     * only with discontinuous finite
+                                     * elements or with no local
+                                     * refinement.
+                                     */
+    MGTransferPrebuilt ();
+
+                                    /**
+                                     * Constructor with constraint matrices.
+                                     */
+    MGTransferPrebuilt (const ConstraintMatrix& constraints);
                                     /**
                                      * Destructor.
                                      */
     virtual ~MGTransferPrebuilt ();
-    
                                     /**
                                      * Actually build the prolongation
                                      * matrices for each level.
@@ -176,6 +189,42 @@ class MGTransferPrebuilt : public MGTransferBase<VECTOR>
                                     */
     std::vector<std::map<unsigned int, unsigned int> >
     copy_indices;
+
+                                    /**
+                                     * Fill the two boolean vectors
+                                     * #dofs_on_refinement_edge and
+                                     * #dofs_on_refinement_boundary.
+                                     */
+    template <int dim, int spacedim>
+    void find_dofs_on_refinement_edges (
+        const MGDoFHandler<dim,spacedim>& mg_dof);
+
+                                    /**
+                                     * Degrees of freedom on the
+                                     * refinement edge excluding
+                                     * those on the boundary.
+                                     *
+                                     * @todo Clean up names
+                                     */
+    std::vector<std::vector<bool> > dofs_on_refinement_edge;
+                                    /**
+                                     * The degrees of freedom on the
+                                     * the refinement edges. For each
+                                     * level (outer vector) and each
+                                     * dof index (inner vector), this
+                                     * bool is true if the level
+                                     * degree of freedom is on the
+                                     * refinement edge towards the
+                                     * lower level.
+                                     *
+                                     * @todo Clean up names
+                                     */
+    std::vector<std::vector<bool> > dofs_on_refinement_boundary; 
+                                    /**
+                                     * The constraints of the global
+                                     * system.
+                                     */
+    SmartPointer<const ConstraintMatrix> constraints;
 };
 
 
index 3fa2e06042be79991bb8cdecc5242eaa58a45d41..804421ac47a4ccb59ac7450d9a5c5ecc74c07e78 100644 (file)
@@ -15,7 +15,6 @@
 #define __deal2__mg_transfer_templates_h
 
 #include <lac/sparse_matrix.h>
-#include <lac/constraint_matrix.h>
 #include <grid/tria_iterator.h>
 #include <fe/fe.h>
 #include <multigrid/mg_base.h>
@@ -79,10 +78,15 @@ MGTransferPrebuilt<VECTOR>::copy_from_mg(
                                       // to the coarse level, but
                                       // have fine level basis
                                       // functions
+  dst = 0;
   for (unsigned int level=0;level<mg_dof_handler.get_tria().n_levels();++level)
+  {
     for (IT i= copy_indices[level].begin();
         i != copy_indices[level].end();++i)
       dst(i->first) = src[level](i->second);
+  }
+  if (constraints != 0)
+    constraints->condense(dst);
 }
 
 
@@ -91,7 +95,7 @@ template <class VECTOR>
 template <int dim, class OutVector, int spacedim>
 void
 MGTransferPrebuilt<VECTOR>::copy_from_mg_add (
-  const MGDoFHandler<dim,spacedim>              &mg_dof_handler,
+  const MGDoFHandler<dim,spacedim>mg_dof_handler,
   OutVector                            &dst,
   const MGLevelObject<VECTOR> &src) const
 {
@@ -108,6 +112,67 @@ MGTransferPrebuilt<VECTOR>::copy_from_mg_add (
       dst(i->first) += src[level](i->second);
 }
 
+template <class VECTOR>
+template <int dim, int spacedim>
+void MGTransferPrebuilt<VECTOR>::find_dofs_on_refinement_edges (
+    const MGDoFHandler<dim,spacedim>& mg_dof) 
+{
+  for (unsigned int level = 0; level<mg_dof.get_tria().n_levels(); ++level)
+  {
+    std::vector<bool> tmp (mg_dof.n_dofs(level));
+    dofs_on_refinement_edge.push_back (tmp);
+    dofs_on_refinement_boundary.push_back (tmp);
+  }
+
+  const unsigned int   dofs_per_cell   = mg_dof.get_fe().dofs_per_cell;
+  const unsigned int   dofs_per_face   = mg_dof.get_fe().dofs_per_face;
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  std::vector<unsigned int> face_dof_indices (dofs_per_face);
+
+  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof.begin(),
+           endc = mg_dof.end();
+
+  for (; cell!=endc; ++cell)
+  {
+    std::vector<bool> cell_dofs(dofs_per_cell);
+    std::vector<bool> boundary_cell_dofs(dofs_per_cell);
+    const unsigned int level = cell->level();
+    cell->get_mg_dof_indices (local_dof_indices);
+    for (unsigned int face_nr=0; 
+        face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+    {
+      typename DoFHandler<dim>::face_iterator face = cell->face(face_nr);
+      if(!cell->at_boundary(face_nr))
+      {
+        //interior face
+        typename MGDoFHandler<dim>::cell_iterator neighbor 
+          = cell->neighbor(face_nr);
+        // Do refinement face
+        // from the coarse side
+        if (neighbor->level() < cell->level())
+        {
+          for(unsigned int j=0; j<dofs_per_face; ++j)
+            cell_dofs[mg_dof.get_fe().face_to_cell_index(j,face_nr)] = true;
+        }
+      }
+      //boundary face
+      else
+        for(unsigned int j=0; j<dofs_per_face; ++j)
+          boundary_cell_dofs[mg_dof.get_fe().face_to_cell_index(j,face_nr)] = true;
+    }//faces
+    for(unsigned int i=0; i<dofs_per_cell; ++i)
+    {
+      if(cell_dofs[i])
+      {
+        dofs_on_refinement_edge[level][local_dof_indices[i]] = true;
+        dofs_on_refinement_boundary[level][local_dof_indices[i]] = true;
+      }
+      if(boundary_cell_dofs[i])
+        dofs_on_refinement_boundary[level][local_dof_indices[i]] = true;
+    }
+  }
+}
 
 template <class VECTOR>
 unsigned int
index 26ec75f3856e6a7741c33df3dc5e0ecb7f506aef..15050e2412704c7246b3a375c503649a1daebcba 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -103,6 +103,7 @@ class Multigrid : public Subscriptor
              const MGSmootherBase<VECTOR>& pre_smooth,
              const MGSmootherBase<VECTOR>& post_smooth,
              Cycle cycle = v_cycle);
+    
                                     /**
                                      * Experimental constructor for
                                      * cases in which no MGDoFHandler
@@ -214,6 +215,37 @@ class Multigrid : public Subscriptor
                                      */
     void Tvmult_add(VECTOR& dst, const VECTOR& src) const;
     
+                                    /**
+                                     * Set additional matrices to
+                                     * correct residual computation
+                                     * at refinement edges. Since we
+                                     * only smoothen in the interior
+                                     * of the refined part of the
+                                     * mesh, the coupling across the
+                                     * refinement edge is
+                                     * missing. This coupling is
+                                     * provided by these two
+                                     * matrices.
+                                     *
+                                     * @note While
+                                     * <tt>edge_out.vmult</tt> is
+                                     * used, for the second argument,
+                                     * we use
+                                     * <tt>edge_in.Tvmult</tt>. Thus,
+                                     * <tt>edge_in</tt> should be
+                                     * assembled in transposed
+                                     * form. This saves a second
+                                     * sparsity pattern for
+                                     * <tt>edge_in</tt>. In
+                                     * particular, for symmetric
+                                     * operators, both arguments can
+                                     * refer to the same matrix,
+                                     * saving assembling of one of
+                                     * them.
+                                     */
+    void set_edge_matrices (const MGMatrixBase<VECTOR>& edge_out,
+                           const MGMatrixBase<VECTOR>& edge_in);
+
                                     /**
                                      * Set additional matrices to
                                      * correct residual computation
@@ -224,9 +256,25 @@ class Multigrid : public Subscriptor
                                      * correspond tu the edge fluxes
                                      * at the refinement edge between
                                      * two levels.
-                                     */
-    void set_edge_matrices (const MGMatrixBase<VECTOR>& edge_down,
-                           const MGMatrixBase<VECTOR>& edge_up);
+                                     *
+                                     * @note While
+                                     * <tt>edge_down.vmult</tt> is
+                                     * used, for the second argument,
+                                     * we use
+                                     * <tt>edge_up.Tvmult</tt>. Thus,
+                                     * <tt>edge_up</tt> should be
+                                     * assembled in transposed
+                                     * form. This saves a second
+                                     * sparsity pattern for
+                                     * <tt>edge_up</tt>. In
+                                     * particular, for symmetric
+                                     * operators, both arguments can
+                                     * refer to the same matrix,
+                                     * saving assembling of one of
+                                     * them.
+                                     */
+    void set_edge_flux_matrices (const MGMatrixBase<VECTOR>& edge_down,
+                                const MGMatrixBase<VECTOR>& edge_up);
 
                                     /**
                                      * Return the finest level for
@@ -396,13 +444,39 @@ class Multigrid : public Subscriptor
                                      */
     SmartPointer<const MGSmootherBase<VECTOR>,Multigrid<VECTOR> > post_smooth;
     
+                                    /**
+                                     * Edge matrix from the interior
+                                     * of the refined part to the
+                                     * refinement edge.
+                                     *
+                                     * @note Only <tt>vmult</tt> is
+                                     * used for these matrices.
+                                     */
+    SmartPointer<const MGMatrixBase<VECTOR> > edge_out;
+
+                                    /**
+                                     * Transpose edge matrix from the
+                                     * refinement edge to the
+                                     * interior of the refined part.
+                                     *
+                                     * @note Only <tt>Tvmult</tt> is
+                                     * used for these matrices.
+                                     */
+    SmartPointer<const MGMatrixBase<VECTOR> > edge_in;
+
                                     /**
                                      * Edge matrix from fine to coarse.
+                                     *
+                                     * @note Only <tt>vmult</tt> is
+                                     * used for these matrices.
                                      */
     SmartPointer<const MGMatrixBase<VECTOR>,Multigrid<VECTOR> > edge_down;
 
                                     /**
                                      * Transpose edge matrix from coarse to fine.
+                                     *
+                                     * @note Only <tt>Tvmult</tt> is
+                                     * used for these matrices.
                                      */
     SmartPointer<const MGMatrixBase<VECTOR>,Multigrid<VECTOR> > edge_up;
 
index e53e8ccc0e26726e9af39a44afd4658f34b48e15..e0710af3889b733756db8cf7ae5b899551bdf2e4 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -42,6 +42,8 @@ Multigrid<VECTOR>::Multigrid (const unsigned int minlevel,
                transfer(&transfer, typeid(*this).name()),
                pre_smooth(&pre_smooth, typeid(*this).name()),
                post_smooth(&post_smooth, typeid(*this).name()),
+               edge_out(0, typeid(*this).name()),
+               edge_in(0, typeid(*this).name()),
                edge_down(0, typeid(*this).name()),
                edge_up(0, typeid(*this).name()),
                debug(0)
@@ -103,11 +105,19 @@ void
 Multigrid<VECTOR>::set_edge_matrices (const MGMatrixBase<VECTOR>& down,
                                      const MGMatrixBase<VECTOR>& up)
 {
-  edge_down = &down;
-  edge_up = &up;
+  edge_out = &down;
+  edge_in = &up;
 }
 
 
+template <class VECTOR>
+void
+Multigrid<VECTOR>::set_edge_flux_matrices (const MGMatrixBase<VECTOR>& down,
+                                          const MGMatrixBase<VECTOR>& up)
+{
+  edge_down = &down;
+  edge_up = &up;
+}
 
 
 template <class VECTOR>
@@ -150,10 +160,12 @@ Multigrid<VECTOR>::level_v_step(const unsigned int level)
   for (unsigned int l = level;l>minlevel;--l)
     {
       t[l-1] = 0.;
+      if (l==level && edge_out != 0)
+       edge_out->vmult_add(level, t[level], solution[level]);
+      
       if (l==level && edge_down != 0)
-       {
-         edge_down->vmult(level, t[level-1], solution[level]);
-       }
+       edge_down->vmult(level, t[level-1], solution[level]);
+      
       transfer->restrict_and_add (l, t[l-1], t[l]);
       defect[l-1] -= t[l-1];
     }
@@ -173,12 +185,18 @@ Multigrid<VECTOR>::level_v_step(const unsigned int level)
   
   solution[level] += t[level];
   
-  if (edge_up != 0)
+  if (edge_in != 0)
     {
-      edge_up->Tvmult(level, t[level], solution[level-1]);
+      edge_in->Tvmult(level, t[level], solution[level]);
       defect[level] -= t[level];
     }
   
+  if (edge_up != 0)
+    {
+     edge_up->Tvmult(level, t[level], solution[level-1]);
+     defect[level] -= t[level];
+   }
+  
   if (debug>2)
     deallog << "V-cycle  Defect norm   " << defect[level].l2_norm()
            << std::endl;
@@ -261,8 +279,12 @@ Multigrid<VECTOR>::level_step(const unsigned int level,
                                   // coarse-level defect already
                                   // contain the global defect, the
                                   // refined parts its restriction.
+  if (edge_out != 0)
+    edge_out->vmult_add(level, t[level], solution[level]);
+  
   if (edge_down != 0)
     edge_down->vmult_add(level, defect2[level-1], solution[level]);
+  
   transfer->restrict_and_add (level, defect2[level-1], t[level]);
                                   // do recursion
   solution[level-1] = 0.;
@@ -296,11 +318,12 @@ Multigrid<VECTOR>::level_step(const unsigned int level,
   solution[level] += t[level];
 
   
-  if (edge_up != 0)
-    {
-      edge_up->Tvmult(level, t[level], solution[level-1]);
-    }
+  if (edge_in != 0)
+    edge_in->Tvmult(level, t[level], solution[level]);
 
+  if (edge_up != 0)
+    edge_up->Tvmult(level, t[level], solution[level-1]);
+  
   t[level].sadd(-1.,-1.,defect2[level],1.,defect[level]);
 
   if (debug>2)
index 763e06e2beb76bbf8854b054c74ae1686635f5ad..22c40f93edfebadb5c2419486b481e17f24dacec 100644 (file)
@@ -907,7 +907,6 @@ MGTools::make_flux_sparsity_pattern_edge (
 }
 
 
-
 template <int dim, int spacedim>
 void
 MGTools::count_dofs_per_component (
diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc
deleted file mode 100644 (file)
index 76a5a9a..0000000
+++ /dev/null
@@ -1,139 +0,0 @@
-//---------------------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//---------------------------------------------------------------------------
-
-
-#include <lac/vector.h>
-#include <lac/block_vector.h>
-#include <lac/sparse_matrix.h>
-#include <lac/block_sparse_matrix.h>
-#include <lac/constraint_matrix.h>
-#include <fe/fe.h>
-#include <grid/tria.h>
-#include <grid/tria_iterator.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
-#include <multigrid/mg_smoother.h>
-#include <multigrid/mg_smoother.templates.h>
-
-#include <algorithm>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-
-#if deal_II_dimension == 1
-
-template <class VECTOR>
-MGSmootherContinuous<VECTOR>::MGSmootherContinuous (
-  const MGDoFHandler<1> &/*mg_dof*/,
-  const unsigned int steps)
-               :
-               steps(steps)
-{
-  Assert (false, ExcNotImplemented());
-}
-
-#endif
-
-
-#if deal_II_dimension > 1
-
-template <class VECTOR>
-template <int dim, int spacedim>
-MGSmootherContinuous<VECTOR>::MGSmootherContinuous (
-  const MGDoFHandler<dim,spacedim> &mg_dof,
-  const unsigned int steps)
-               :
-               steps(steps)
-{
-  const unsigned int n_levels = mg_dof.get_tria().n_levels();
-  
-                                  // allocate the right number of
-                                  // elements
-  interior_boundary_dofs.resize (n_levels-1);
-
-                                  // use a temporary to store the
-                                  // indices. this allows to not allocate
-                                  // to much memory by later copying the
-                                  // content of this vector to its final
-                                  // destination
-  std::vector<unsigned int> boundary_dofs;
-
-                                  // temporary to hold the dof indices
-                                  // on a face between to levels
-  std::vector<unsigned int> dofs_on_face (mg_dof.get_fe().dofs_per_face);
-
-  for (unsigned int level=1; level<n_levels; ++level)
-    {
-      boundary_dofs.clear ();
-
-                                      // for each cell on this level:
-                                      // find out whether a face is
-                                      // at the boundary of this level's
-                                      // cells and if so add the dofs
-                                      // to the interior boundary dofs
-      for (typename MGDoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
-          cell != mg_dof.end(level); ++cell)
-       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-         if ((cell->neighbor(face).state() == IteratorState::valid) &&
-             (static_cast<unsigned int>(cell->neighbor(face)->level())
-              == level-1))
-           {
-                                              // get indices of this face
-             cell->face(face)->get_mg_dof_indices (level, dofs_on_face);
-                                              // append them to the levelwise
-                                              // list
-             boundary_dofs.insert (boundary_dofs.end(),
-                                   dofs_on_face.begin(),
-                                   dofs_on_face.end());
-           };
-
-                                      // now sort the list of interior boundary
-                                      // dofs and eliminate duplicates
-      std::sort (boundary_dofs.begin(), boundary_dofs.end());
-      boundary_dofs.erase (std::unique (boundary_dofs.begin(),
-                                       boundary_dofs.end()),
-                          boundary_dofs.end());
-
-                                      // now finally copy the result
-                                      // for this level its destination
-      interior_boundary_dofs[level-1] = boundary_dofs;
-    };      
-}
-
-#endif
-
-
-// explicit instantiations
-
-// don't do the following instantiation in 1d, since there is a specialized
-// function there
-#if deal_II_dimension > 1
-template MGSmootherContinuous<Vector<float> >::MGSmootherContinuous (
-  const MGDoFHandler<deal_II_dimension>&,
-  const unsigned int);
-template MGSmootherContinuous<Vector<double> >::MGSmootherContinuous (
-  const MGDoFHandler<deal_II_dimension>&,
-  const unsigned int);
-template MGSmootherContinuous<BlockVector<float> >::MGSmootherContinuous (
-  const MGDoFHandler<deal_II_dimension>&,
-  const unsigned int);
-template MGSmootherContinuous<BlockVector<double> >::MGSmootherContinuous (
-  const MGDoFHandler<deal_II_dimension>&,
-  const unsigned int);
-#endif
-
-DEAL_II_NAMESPACE_CLOSE
index 0ff5d6ca11df5a91fcd7a0825e260403b729cb12..e85e42f8fde3cd53039de03d1b769d5a1b010396 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -32,7 +32,8 @@ void
 MGTools::apply_boundary_values (
   const std::set<unsigned int> &boundary_dofs,
   SparseMatrix<number>& matrix,
-  const bool preserve_symmetry)
+  const bool preserve_symmetry,
+  const bool ignore_zeros)
 {
                                   // if no boundary values are to be applied
                                   // simply return
@@ -95,10 +96,9 @@ MGTools::apply_boundary_values (
                                       //
                                       // store the new rhs entry to make
                                       // the gauss step more efficient
-      matrix.set (dof_number, dof_number,
-                 first_nonzero_diagonal_entry);
-
-
+      if(!ignore_zeros)
+       matrix.set (dof_number, dof_number,
+                   first_nonzero_diagonal_entry);
                                       // if the user wants to have
                                       // the symmetry of the matrix
                                       // preserved, and if the
@@ -431,10 +431,10 @@ MGTools::apply_boundary_values (
 
 template void MGTools::apply_boundary_values (
   const std::set<unsigned int>&,
-  SparseMatrix<float>&, const bool);
+  SparseMatrix<float>&, const bool, const bool);
 template void MGTools::apply_boundary_values (
   const std::set<unsigned int>&,
-  SparseMatrix<double>&, const bool);
+  SparseMatrix<double>&, const bool, const bool);
 template void MGTools::apply_boundary_values (
   const std::set<unsigned int>&,
   BlockSparseMatrix<float>&, const bool);
index a2e7137008349d344a85098954c018adb1fbec83..04921170bd0c9475c2f41149876de5ae11cea7ea 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -12,6 +12,7 @@
 //---------------------------------------------------------------------------
 
 #include <base/logstream.h>
+#include <base/function.h>
 
 #include <lac/vector.h>
 #include <lac/block_vector.h>
@@ -23,6 +24,7 @@
 #include <fe/fe.h>
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_tools.h>
 #include <multigrid/mg_transfer.h>
 #include <multigrid/mg_transfer.templates.h>
 
@@ -35,6 +37,7 @@ template <int dim, int spacedim>
 void MGTransferPrebuilt<number>::build_matrices (
   const MGDoFHandler<dim,spacedim> &mg_dof)
 {
+  //start building the matrices here
   const unsigned int n_levels      = mg_dof.get_tria().n_levels();
   const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
   
@@ -156,9 +159,42 @@ void MGTransferPrebuilt<number>::build_matrices (
              }
          }
     }
-                                  // For non-DG elements, add a
-                                  // condense operator here.
+  // impose boundary conditions 
+  // but only in the column of 
+  // the prolongation matrix 
+  //TODO: this way is not very efficient
+  std::vector<std::set<unsigned int> > boundary_indices (mg_dof.get_tria().n_levels());
+  typename FunctionMap<dim>::type boundary;
+  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
+  boundary[0] = &homogeneous_dirichlet_bc;
+  MGTools::make_boundary_list(mg_dof, boundary, boundary_indices);
+
+  for (unsigned int level=0; level<n_levels-1; ++level)
+  {
+    if (boundary_indices[level].size() == 0)
+      continue;
+
+    const unsigned int n_dofs = prolongation_matrices[level]->m();
+    for (unsigned int i=0; i<n_dofs; ++i)
+    {
+      SparseMatrix<double>::iterator anfang = prolongation_matrices[level]->begin(i),
+        ende = prolongation_matrices[level]->end(i); 
+      for(; anfang != ende; ++anfang)
+      {
+        std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
+          endd = boundary_indices[level].end();
+        for (; dof != endd; ++dof)
+        {
+          const unsigned int column_number = anfang->column();
+          const unsigned int dof_number = *dof;
+          if(dof_number == column_number)
+            prolongation_matrices[level]->set(i,dof_number,0);
+        }
+      }
+    }
+  }
   
+  find_dofs_on_refinement_edges (mg_dof);
                                   // Prepare indices to copy from a
                                   // global vector to the parts of an
                                   // MGVector corresponding to active
@@ -188,6 +224,7 @@ void MGTransferPrebuilt<number>::build_matrices (
          
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
+              if(!dofs_on_refinement_edge[level][level_dof_indices[i]])
              copy_indices[level].insert(
                std::make_pair(global_dof_indices[i], level_dof_indices[i]));
            }
index 92966843a22c076ee655f85b8f0d768308a30a0b..b7508931f823bbc9f44d44ea9dfd01dd255305a3 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -101,6 +101,16 @@ Multigrid<VECTOR>::print_vector (const unsigned int level,
 #endif
 
 
+template<class VECTOR>
+MGTransferPrebuilt<VECTOR>::MGTransferPrebuilt ()
+{} 
+
+
+template<class VECTOR>
+MGTransferPrebuilt<VECTOR>::MGTransferPrebuilt (const ConstraintMatrix &c)
+ :
+   constraints(&c)
+{} 
 
 
 template <class VECTOR>
@@ -133,10 +143,6 @@ void MGTransferPrebuilt<VECTOR>::restrict_and_add (
   prolongation_matrices[from_level-1]->Tvmult_add (dst, src);
 }
 
-
-
-
-
 template <typename number>
 MGTransferBlock<number>::MGTransferBlock ()
                :
@@ -312,21 +318,6 @@ void MGTransferSelect<number>::restrict_and_add (
 }
 
 
-template <class VECTOR>
-void
-MGSmootherContinuous<VECTOR>::set_zero_interior_boundary (
-  const unsigned int level,
-  VECTOR&            u) const
-{
-  if (level==0)
-    return;
-  else
-    for (std::vector<unsigned int>::const_iterator p=interior_boundary_dofs[level-1].begin();
-        p!=interior_boundary_dofs[level-1].end(); ++p)
-      u(*p) = 0;
-}
-
-
 //----------------------------------------------------------------------//
 
 template <typename number>
@@ -383,17 +374,5 @@ template class MGTransferSelect<double>;
 template class MGTransferBlockSelect<float>;
 template class MGTransferBlockSelect<double>;
 
-template
-void MGSmootherContinuous<Vector<double> >::set_zero_interior_boundary (
-  const unsigned int, Vector<double>&) const;
-template
-void MGSmootherContinuous<Vector<float> >::set_zero_interior_boundary (
-  const unsigned int, Vector<float>&) const;
-template
-void MGSmootherContinuous<BlockVector<double> >::set_zero_interior_boundary (
-  const unsigned int, BlockVector<double>&) const;
-template
-void MGSmootherContinuous<BlockVector<float> >::set_zero_interior_boundary (
-  const unsigned int, BlockVector<float>&) const;
 
 DEAL_II_NAMESPACE_CLOSE
index 2961469b4627c9c06f03833f6a9b18344923d67b..ca908a2997bf40726c4dc25a02dd7465033d487c 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors */
+/*    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -533,7 +533,7 @@ Step39<dim>::solve()
   Multigrid<Vector<double> > mg(mg_dof_handler, mgmatrix,
                                mg_coarse, mg_transfer,
                                mg_smoother, mg_smoother);
-  mg.set_edge_matrices(mgdown, mgup);
+  mg.set_edge_flux_matrices(mgdown, mgup);
   mg.set_debug(0);
   mg_smoother.set_debug(0);
   

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.