]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
automatically create quadrature if none is provided
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 8 Oct 2010 20:48:45 +0000 (20:48 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 8 Oct 2010 20:48:45 +0000 (20:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@22295 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/mesh_worker_info.h
deal.II/examples/step-39/step-39.cc
deal.II/lac/include/lac/matrix_block.h

index ced843f5f9d9fb63b555dff99e1ed99c7a7ed472..f70e2449972365e0b80a177ee07751be16e9ab2e 100644 (file)
@@ -526,6 +526,10 @@ namespace MeshWorker
  * <li> It stores information on finite element vectors and whether
  * their data should be used to compute values or derivatives of
  * functions at quadrature points.
+ *
+ * <li> It makes educated guesses on quadrature rules and update
+ * flags, so that minimal code has to be written when default
+ * parameters are sufficient.
  * </ol>
  *
  * In order to allow for sufficient generality, a few steps have to be
@@ -541,10 +545,16 @@ namespace MeshWorker
  * the selectors filled before and adds all the flags needed to get
  * the selection. Additional flags can be set with add_update_flags().
  *
- * Finally, we need to choose quadrature formulas. If you choose to
- * use Gauss formulas only, use initialize_gauss_quadrature() with
- * appropriate values. Otherwise, you can fill the variables
- * #cell_quadrature, #boundary_quadrature and #face_quadrature directly.
+ * Finally, we need to choose quadrature formulas. In the simplest
+ * case, you might be happy with the default settings, which are
+ * <i>n</i>-point Gauss formulas. If only derivatives of the shape
+ * functions are used (#update_values is not set) <i>n</i> equals the
+ * highest polynomial degree in the FiniteElement, if #update_values
+ * is set, <i>n</i> is one higher than this degree.  If you choose to
+ * use Gauss formulas of other size, use initialize_gauss_quadrature()
+ * with appropriate values. Otherwise, you can fill the variables
+ * #cell_quadrature, #boundary_quadrature and #face_quadrature
+ * directly.
  *
  * In order to save time, you can set the variables boundary_fluxes
  * and interior_fluxes of the base class to false, thus telling the
@@ -571,17 +581,56 @@ namespace MeshWorker
                                        * Default constructor.
                                        */
       IntegrationInfoBox ();
-      
+
+                                      /**
+                                       * Initialize the
+                                       * IntegrationInfo objects
+                                       * contained.
+                                       *
+                                       * Before doing so, add update
+                                       * flags necessary to produce
+                                       * the data needed and also
+                                       * set uninitialized quadrature
+                                       * rules to Gauss formulas,
+                                       * which integrate polynomial
+                                       * bilinear forms exactly.
+                                       */
       void initialize(const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping,
                      const BlockInfo* block_info = 0);
 
+                                      /**
+                                       * Initialize the
+                                       * IntegrationInfo objects
+                                       * contained.
+                                       *
+                                       * Before doing so, add update
+                                       * flags necessary to produce
+                                       * the data needed and also
+                                       * set uninitialized quadrature
+                                       * rules to Gauss formulas,
+                                       * which integrate polynomial
+                                       * bilinear forms exactly.
+                                       */
       template <typename VECTOR>
       void initialize(const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping,
                      const NamedData<VECTOR*>& data,
                      const BlockInfo* block_info = 0);
 
+                                      /**
+                                       * Initialize the
+                                       * IntegrationInfo objects
+                                       * contained.
+                                       *
+                                       * Before doing so, add update
+                                       * flags necessary to produce
+                                       * the data needed and also
+                                       * set uninitialized quadrature
+                                       * rules to Gauss formulas,
+                                       * which integrate polynomial
+                                       * bilinear forms exactly.
+                                       */
       template <typename VECTOR>
       void initialize(const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping,
@@ -637,17 +686,26 @@ namespace MeshWorker
                            const bool face = true,
                            const bool neighbor = true);
 
-                                      /** Assign n-point Gauss
+                                      /**
+                                       * Assign n-point Gauss
                                        * quadratures to each of the
                                        * quadrature rules. Here, a
                                        * size of zero points means
                                        * that no loop over these grid
                                        * entities should be
                                        * performed.
+                                       *
+                                       * If the parameter
+                                       * <tt>force</tt> is true, then
+                                       * all quadrature sets are
+                                       * filled with new quadrature
+                                       * ruels. If it is false, then
+                                       * only empty rules are changed.
                                        */
       void initialize_gauss_quadrature(unsigned int n_cell_points,
                                       unsigned int n_boundary_points,
-                                      unsigned int n_face_points);
+                                      unsigned int n_face_points,
+                                      const bool force = true);
 
                                       /**
                                        * The memory used by this object.
@@ -1303,42 +1361,105 @@ namespace MeshWorker
   }
 
 
+//----------------------------------------------------------------------//
+  
+  template <>
+  inline
+  void
+  IntegrationInfoBox<1,1>::initialize_gauss_quadrature(
+    const unsigned int cp,
+    const unsigned int,
+    const unsigned int,
+    bool force)
+  {
+    if (force || cell_quadrature.size() == 0)
+      cell_quadrature = QGauss<1>(cp);
+  }
+
+
+  template <>
+  inline
+  void
+  IntegrationInfoBox<1,2>::initialize_gauss_quadrature(
+    const unsigned int cp,
+    const unsigned int,
+    const unsigned int,
+    bool force)
+  {
+    if (force || cell_quadrature.size() == 0)
+      cell_quadrature = QGauss<1>(cp);
+  }
+
+  
   template <int dim, int sdim>
   inline
   void
-  IntegrationInfoBox<dim,sdim>::
-  initialize(const FiniteElement<dim,sdim>& el,
-            const Mapping<dim,sdim>& mapping,
-            const BlockInfo* block_info)
+  IntegrationInfoBox<dim,sdim>::initialize_gauss_quadrature(
+    unsigned int cp,
+    unsigned int bp,
+    unsigned int fp,
+    bool force)
   {
-    initialize_update_flags();
-    
-    cell.template initialize<FEValues<dim,sdim> >(el, mapping, cell_quadrature,
-                                                 cell_flags, block_info);
-    boundary.template initialize<FEFaceValues<dim,sdim> >(el, mapping, boundary_quadrature,
-                                                         boundary_flags, block_info);
-    face.template initialize<FEFaceValues<dim,sdim> >(el, mapping, face_quadrature,
-                                                     face_flags, block_info);
-    subface.template initialize<FESubfaceValues<dim,sdim> >(el, mapping, face_quadrature,
-                                                           face_flags, block_info);
-    neighbor.template initialize<FEFaceValues<dim,sdim> >(el, mapping, face_quadrature,
-                                                         neighbor_flags, block_info);
+    if (force || cell_quadrature.size() == 0)
+      cell_quadrature = QGauss<dim>(cp);
+    if (force || boundary_quadrature.size() == 0)
+      boundary_quadrature = QGauss<dim-1>(bp);
+    if (force || face_quadrature.size() == 0)
+      face_quadrature = QGauss<dim-1>(fp);
+  }
+
+
+  template <int dim, int sdim>
+  inline
+  void
+  IntegrationInfoBox<dim,sdim>::add_update_flags_all (const UpdateFlags flags)
+  {
+    add_update_flags(flags, true, true, true, true);
+  }
+
+  
+  template <int dim, int sdim>
+  inline
+  void
+  IntegrationInfoBox<dim,sdim>::add_update_flags_cell (const UpdateFlags flags)
+  {
+    add_update_flags(flags, true, false, false, false);
   }
 
+  
+  template <int dim, int sdim>
+  inline
+  void
+  IntegrationInfoBox<dim,sdim>::add_update_flags_boundary (const UpdateFlags flags)
+  {
+    add_update_flags(flags, false, true, false, false);
+  }
 
+  
+  template <int dim, int sdim>
+  inline
+  void
+  IntegrationInfoBox<dim,sdim>::add_update_flags_face (const UpdateFlags flags)
+  {
+    add_update_flags(flags, false, false, true, true);
+  }
+  
+  
   template <>
   inline
   void
-  IntegrationInfoBox<1,1>::
-  initialize(const FiniteElement<1,1>& el,
-            const Mapping<1,1>& mapping,
-            const BlockInfo* block_info)
+  IntegrationInfoBox<1,1>::initialize(
+    const FiniteElement<1,1>& el,
+    const Mapping<1,1>& mapping,
+    const BlockInfo* block_info)
   {
     initialize_update_flags();
+    initialize_gauss_quadrature(
+      (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), 1, 1, false);
     
     const int dim = 1;
     const int sdim = 1;
-
+    
     cell.initialize<FEValues<dim,sdim> >(el, mapping, cell_quadrature,
                                         cell_flags, block_info);
   }
@@ -1348,12 +1469,14 @@ namespace MeshWorker
   template <>
   inline
   void
-  IntegrationInfoBox<1,2>::
-  initialize(const FiniteElement<1,2>& el,
-            const Mapping<1,2>& mapping,
-            const BlockInfo* block_info)
+  IntegrationInfoBox<1,2>::initialize(
+    const FiniteElement<1,2>& el,
+    const Mapping<1,2>& mapping,
+    const BlockInfo* block_info)
   {
     initialize_update_flags();
+    initialize_gauss_quadrature(
+      (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), 1, 1, false);
     
     const int dim = 1;
     const int sdim = 2;
@@ -1363,7 +1486,32 @@ namespace MeshWorker
   }
 
 
-//----------------------------------------------------------------------//
+  template <int dim, int sdim>
+  inline
+  void
+  IntegrationInfoBox<dim,sdim>::initialize(
+    const FiniteElement<dim,sdim>& el,
+    const Mapping<dim,sdim>& mapping,
+    const BlockInfo* block_info)
+  {
+    initialize_update_flags();
+    initialize_gauss_quadrature(
+      (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(),
+      (boundary_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(),
+      (face_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), false);
+    
+    cell.template initialize<FEValues<dim,sdim> >(el, mapping, cell_quadrature,
+                                                 cell_flags, block_info);
+    boundary.template initialize<FEFaceValues<dim,sdim> >(el, mapping, boundary_quadrature,
+                                                         boundary_flags, block_info);
+    face.template initialize<FEFaceValues<dim,sdim> >(el, mapping, face_quadrature,
+                                                     face_flags, block_info);
+    subface.template initialize<FESubfaceValues<dim,sdim> >(el, mapping, face_quadrature,
+                                                           face_flags, block_info);
+    neighbor.template initialize<FEFaceValues<dim,sdim> >(el, mapping, face_quadrature,
+                                                         neighbor_flags, block_info);
+  }
+
 
   template <int dim, int sdim>
   template <typename VECTOR>
@@ -1441,79 +1589,6 @@ namespace MeshWorker
   {}
 
 
-  template <int dim, int sdim>
-  inline
-  void
-  IntegrationInfoBox<dim,sdim>::initialize_gauss_quadrature(
-    unsigned int cp,
-    unsigned int bp,
-    unsigned int fp)
-  {
-    cell_quadrature = QGauss<dim>(cp);
-    boundary_quadrature = QGauss<dim-1>(bp);
-    face_quadrature = QGauss<dim-1>(fp);
-  }
-
-
-  template <>
-  inline
-  void
-  IntegrationInfoBox<1,1>::
-  initialize_gauss_quadrature(const unsigned int cp,
-                             const unsigned int,
-                             const unsigned int)
-  {
-    cell_quadrature = QGauss<1>(cp);
-  }
-
-
-  template <>
-  inline
-  void
-  IntegrationInfoBox<1,2>::
-  initialize_gauss_quadrature(const unsigned int cp,
-                             const unsigned int,
-                             const unsigned int)
-  {
-    cell_quadrature = QGauss<1>(cp);
-  }
-
-  
-  template <int dim, int sdim>
-  inline
-  void
-  IntegrationInfoBox<dim,sdim>::add_update_flags_all (const UpdateFlags flags)
-  {
-    add_update_flags(flags, true, true, true, true);
-  }
-
-  
-  template <int dim, int sdim>
-  inline
-  void
-  IntegrationInfoBox<dim,sdim>::add_update_flags_cell (const UpdateFlags flags)
-  {
-    add_update_flags(flags, true, false, false, false);
-  }
-
-  
-  template <int dim, int sdim>
-  inline
-  void
-  IntegrationInfoBox<dim,sdim>::add_update_flags_boundary (const UpdateFlags flags)
-  {
-    add_update_flags(flags, false, true, false, false);
-  }
-
-  
-  template <int dim, int sdim>
-  inline
-  void
-  IntegrationInfoBox<dim,sdim>::add_update_flags_face (const UpdateFlags flags)
-  {
-    add_update_flags(flags, false, false, true, true);
-  }
-  
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 335bec6d7810c5b78c9751d20d1ebe567837bf2b..e47b9deff3cb1306b6891d2ccf6b0fd511f014a3 100644 (file)
@@ -570,27 +570,33 @@ Step39<dim>::assemble_matrix()
                                   // such that they always point to
                                   // the current cell. To this end,
                                   // we need to tell it first, where
-                                  // and what to compute,
+                                  // and what to compute. Since we
+                                  // are not doing anything fancy, we
+                                  // can rely on their standard
+                                  // choice for quadrature rules.
+                                  //
+                                  // Since their default update flags
+                                  // are minimal, we add what we need
+                                  // additionally, namely the values
+                                  // and gradients of shape functions
+                                  // on all objects (cells, boundary
+                                  // and interior faces). Afterwards,
+                                  // we are ready to initialize the
+                                  // container, which will create all
+                                  // necessary FEValuesBase objects
+                                  // for integration.
   MeshWorker::IntegrationInfoBox<dim> info_box;
-                                  // namely, which quadrature
-                                  // formulas to use and
-  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
-  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
-                                  // which values to update in these
-                                  // points. Update flags are
-                                  // initialized to some default
-                                  // values to be able to
-                                  // integrate. Here, we add what we
-                                  // need additionally, namely the
-                                  // values and gradients of shape
-                                  // functions on all objects (cells,
-                                  // boundary and interior faces).
   UpdateFlags update_flags = update_values | update_gradients;
   info_box.add_update_flags_all(update_flags);
   info_box.initialize(fe, mapping);
 
                                   // This is the object into which we
-                                  // integrate local data.
+                                  // integrate local data. It is
+                                  // filled by the local integration
+                                  // routines in MatrixIntegrator and
+                                  // then used by the assembler to
+                                  // distribute the information into
+                                  // the global matrix.
   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
 
                                   // Finally, we need an object that
@@ -628,8 +634,6 @@ void
 Step39<dim>::assemble_mg_matrix()
 {
   MeshWorker::IntegrationInfoBox<dim> info_box;
-  const unsigned int n_gauss_points = mg_dof_handler.get_fe().tensor_degree()+1;
-  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
   UpdateFlags update_flags = update_values | update_gradients;
   info_box.add_update_flags_all(update_flags);
   info_box.initialize(fe, mapping);
@@ -673,8 +677,6 @@ void
 Step39<dim>::assemble_right_hand_side()
 {
   MeshWorker::IntegrationInfoBox<dim> info_box;
-  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
-  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
   UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
   info_box.add_update_flags_all(update_flags);
   info_box.initialize(fe, mapping);
index 3560b7d993cbef9cead670ac04828df8f4f4b238..1dba79fd2a72dc3a884f7a678432103ec90047e5 100644 (file)
@@ -479,13 +479,35 @@ class MatrixBlockVector
  */
 template <class MATRIX>
 class MGMatrixBlockVector
-  : private NamedData<std_cxx1x::shared_ptr<MGLevelObject<MatrixBlock<MATRIX> > > >
+  : public Subscriptor
 {
   public:
                                     /**
                                      * The type of object stored.
                                      */
     typedef MGLevelObject<MatrixBlock<MATRIX> > value_type;
+                                    /**
+                                     * Constructor, determining which
+                                     * matrices should be stored.
+                                     *
+                                     * If <tt>edge_matrices</tt> is
+                                     * true, then objects for edge
+                                     * matrices for discretizations
+                                     * with degrees of freedom on
+                                     * faces are allocated.
+                                     *
+                                     * If <tt>edge_flux_matrices</tt>
+                                     * is true, then objects for DG
+                                     * fluxes on the refinement edge
+                                     * are allocated.
+                                     */
+    MGMatrixBlockVector(const bool edge_matrices = false,
+                       const bool edge_flux_matrices = false);
+
+                                    /**
+                                     * The number of blocks.
+                                     */
+    unsigned int size () const;
     
                                     /**
                                      * Add a new matrix block at the
@@ -504,8 +526,35 @@ class MGMatrixBlockVector
                                      * reinitializes each matrix in
                                      * the vector with the correct
                                      * pattern from the block system.
+                                     *
+                                     * This function reinitializes
+                                     * the level matrices.
                                      */
-    void reinit(const MGLevelObject<BlockSparsityPattern>& sparsity);
+    void reinit_matrix(const MGLevelObject<BlockSparsityPattern>& sparsity);
+                                    /**
+                                     * For matrices using a
+                                     * SparsityPattern, this function
+                                     * reinitializes each matrix in
+                                     * the vector with the correct
+                                     * pattern from the block system.
+                                     *
+                                     * This function reinitializes
+                                     * the matrices for degrees of
+                                     * freedom on the refinement edge.
+                                     */
+    void reinit_edge(const MGLevelObject<BlockSparsityPattern>& sparsity);
+                                    /**
+                                     * For matrices using a
+                                     * SparsityPattern, this function
+                                     * reinitializes each matrix in
+                                     * the vector with the correct
+                                     * pattern from the block system.
+                                     *
+                                     * This function reinitializes
+                                     * the flux matrices over the
+                                     * refinement edge.
+                                     */
+    void reinit_edge_flux(const MGLevelObject<BlockSparsityPattern>& sparsity);
     
                                     /**
                                      * Clears the object.
@@ -528,26 +577,88 @@ class MGMatrixBlockVector
     
                                     /**
                                      * Access a constant reference to
-                                     * the block at position <i>i</i>.
+                                     * the matrix block at position <i>i</i>.
                                      */
     const value_type& block(unsigned int i) const;
     
                                     /**
                                      * Access a reference to
-                                     * the block at position <i>i</i>.
+                                     * the matrix block at position <i>i</i>.
                                      */
     value_type& block(unsigned int i);
     
+                                    /**
+                                     * Access a constant reference to
+                                     * the edge matrix block at position <i>i</i>.
+                                     */
+    const value_type& block_in(unsigned int i) const;
+    
+                                    /**
+                                     * Access a reference to
+                                     * the edge matrix block at position <i>i</i>.
+                                     */
+    value_type& block_in(unsigned int i);
+    
+                                    /**
+                                     * Access a constant reference to
+                                     * the edge matrix block at position <i>i</i>.
+                                     */
+    const value_type& block_out(unsigned int i) const;
+    
+                                    /**
+                                     * Access a reference to
+                                     * the edge matrix block at position <i>i</i>.
+                                     */
+    value_type& block_out(unsigned int i);
+    
+                                    /**
+                                     * Access a constant reference to
+                                     * the  edge flux matrix block at position <i>i</i>.
+                                     */
+    const value_type& block_up(unsigned int i) const;
+    
+                                    /**
+                                     * Access a reference to
+                                     * the  edge flux matrix block at position <i>i</i>.
+                                     */
+    value_type& block_up(unsigned int i);
+    
+                                    /**
+                                     * Access a constant reference to
+                                     * the  edge flux matrix block at position <i>i</i>.
+                                     */
+    const value_type& block_down(unsigned int i) const;
+    
+                                    /**
+                                     * Access a reference to
+                                     * the edge flux matrix block at position <i>i</i>.
+                                     */
+    value_type& block_down(unsigned int i);
+    
                                     /**
                                      * The memory used by this object.
                                      */
     unsigned int memory_consumption () const;
-
-    NamedData<std_cxx1x::shared_ptr<value_type> >::subscribe;
-    NamedData<std_cxx1x::shared_ptr<value_type> >::unsubscribe;
-    NamedData<std_cxx1x::shared_ptr<value_type> >::size;
-
   private:
+                                    /// Clear one of the matrix objects
+    void clear_object(NamedData<MGLevelObject<MatrixBlock<MATRIX> > >&);
+    
+                                    /// Flag for storing #matrices_in and #matrices_out
+    const bool edge_matrices;
+
+                                    /// Flag for storing #flux_matrices_up and #flux_matrices_down
+    const bool edge_flux_matrices;
+    
+                                    /// The level matrices
+    NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices;
+              /// The matrix from the interior of a level to the refinement edge
+    NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices_in;
+              /// The matrix from the refinement edge to the interior of a level
+    NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices_out;
+              /// The DG flux from a level to the lower level
+    NamedData<MGLevelObject<MatrixBlock<MATRIX> > > flux_matrices_down;
+              /// The DG flux from the lower level to a level
+    NamedData<MGLevelObject<MatrixBlock<MATRIX> > > flux_matrices_up;
 };
 
 
@@ -871,18 +982,46 @@ MatrixBlockVector<MATRIX>::matrix(unsigned int i)
 
 //----------------------------------------------------------------------//
 
+template <class MATRIX>
+inline
+MGMatrixBlockVector<MATRIX>::MGMatrixBlockVector(
+  const bool e, const bool f)
+               :
+               edge_matrices(e),
+               edge_flux_matrices(f)
+{}
+
+
+template <class MATRIX>
+inline
+unsigned int
+MGMatrixBlockVector<MATRIX>::size () const
+{
+  return matrices.size();
+}
+
+
 template <class MATRIX>
 inline void
 MGMatrixBlockVector<MATRIX>::add(
   unsigned int row, unsigned int column,
   const std::string& name)
 {
-  std_cxx1x::shared_ptr<MGLevelObject<MatrixBlock<MATRIX> > >
-    p(new value_type(0, 1));
-  (*p)[0].row = row;
-  (*p)[0].column = column;
+  MGLevelObject<MatrixBlock<MATRIX> > p(0, 1);
+  p[0].row = row;
+  p[0].column = column;
   
-  NamedData<std_cxx1x::shared_ptr<value_type> >::add(p, name);
+  matrices.add(p, name);
+  if (edge_matrices)
+    {
+      matrices_in.add(p, name);
+      matrices_out.add(p, name);
+    }
+  if (edge_flux_matrices)
+    {
+      flux_matrices_up.add(p, name);
+      flux_matrices_down.add(p, name);
+    }
 }
 
 
@@ -890,7 +1029,7 @@ template <class MATRIX>
 inline const MGLevelObject<MatrixBlock<MATRIX> >&
 MGMatrixBlockVector<MATRIX>::block(unsigned int i) const
 {
-  return *this->read(i);
+  return matrices.read(i);
 }
 
 
@@ -898,13 +1037,77 @@ template <class MATRIX>
 inline MGLevelObject<MatrixBlock<MATRIX> >&
 MGMatrixBlockVector<MATRIX>::block(unsigned int i)
 {
-  return *(*this)(i);
+  return matrices(i);
+}
+
+
+template <class MATRIX>
+inline const MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_in(unsigned int i) const
+{
+  return matrices_in.read(i);
+}
+
+
+template <class MATRIX>
+inline MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_in(unsigned int i)
+{
+  return matrices_in(i);
+}
+
+
+template <class MATRIX>
+inline const MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_out(unsigned int i) const
+{
+  return matrices_out.read(i);
+}
+
+
+template <class MATRIX>
+inline MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_out(unsigned int i)
+{
+  return matrices_out(i);
+}
+
+
+template <class MATRIX>
+inline const MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_up(unsigned int i) const
+{
+  return flux_matrices_up.read(i);
+}
+
+
+template <class MATRIX>
+inline MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_up(unsigned int i)
+{
+  return flux_matrices_up(i);
+}
+
+
+template <class MATRIX>
+inline const MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_down(unsigned int i) const
+{
+  return flux_matrices_down.read(i);
+}
+
+
+template <class MATRIX>
+inline MGLevelObject<MatrixBlock<MATRIX> >&
+MGMatrixBlockVector<MATRIX>::block_down(unsigned int i)
+{
+  return flux_matrices_down(i);
 }
 
 
 template <class MATRIX>
 inline void
-MGMatrixBlockVector<MATRIX>::reinit(const MGLevelObject<BlockSparsityPattern>& sparsity)
+MGMatrixBlockVector<MATRIX>::reinit_matrix(const MGLevelObject<BlockSparsityPattern>& sparsity)
 {
   for (unsigned int i=0;i<this->size();++i)
     {
@@ -923,6 +1126,70 @@ MGMatrixBlockVector<MATRIX>::reinit(const MGLevelObject<BlockSparsityPattern>& s
 }
 
 
+template <class MATRIX>
+inline void
+MGMatrixBlockVector<MATRIX>::reinit_edge(const MGLevelObject<BlockSparsityPattern>& sparsity)
+{
+  for (unsigned int i=0;i<this->size();++i)
+    {
+      MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
+      const unsigned int row = o[o.min_level()].row;
+      const unsigned int col = o[o.min_level()].column;
+
+      block_in(i).resize(sparsity.min_level(), sparsity.max_level());
+      block_out(i).resize(sparsity.min_level(), sparsity.max_level());
+      for (unsigned int level = o.min_level();level <= o.max_level();++level)
+       {
+         block_in(i)[level].row = row;
+         block_in(i)[level].column = col;
+         internal::reinit(block_in(i)[level], sparsity[level]);
+         block_out(i)[level].row = row;
+         block_out(i)[level].column = col;
+         internal::reinit(block_out(i)[level], sparsity[level]);
+       }
+    }
+}
+
+
+template <class MATRIX>
+inline void
+MGMatrixBlockVector<MATRIX>::reinit_edge_flux(const MGLevelObject<BlockSparsityPattern>& sparsity)
+{
+  for (unsigned int i=0;i<this->size();++i)
+    {
+      MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
+      const unsigned int row = o[o.min_level()].row;
+      const unsigned int col = o[o.min_level()].column;
+
+      block_up(i).resize(sparsity.min_level(), sparsity.max_level());
+      block_down(i).resize(sparsity.min_level(), sparsity.max_level());
+      for (unsigned int level = o.min_level();level <= o.max_level();++level)
+       {
+         block_up(i)[level].row = row;
+         block_up(i)[level].column = col;
+         internal::reinit(block_up(i)[level], sparsity[level]);
+         block_down(i)[level].row = row;
+         block_down(i)[level].column = col;
+         internal::reinit(block_down(i)[level], sparsity[level]);
+       }
+      
+    }
+}
+
+
+template <class MATRIX>
+inline void
+MGMatrixBlockVector<MATRIX>::clear_object(NamedData<MGLevelObject<MatrixBlock<MATRIX> > >& mo)
+{
+  for (unsigned int i=0;i<mo.size();++i)
+    {
+      MGLevelObject<MatrixBlock<MATRIX> >& o = mo(i);
+      for (unsigned int level = o.min_level();level <= o.max_level();++level)
+       o[level].matrix.clear();
+    }
+}
+
+
 template <class MATRIX>
 inline void
 MGMatrixBlockVector<MATRIX>::clear(bool really_clean)
@@ -933,12 +1200,11 @@ MGMatrixBlockVector<MATRIX>::clear(bool really_clean)
     }
   else
     {
-      for (unsigned int i=0;i<this->size();++i)
-       {
-         MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
-         for (unsigned int level = o.min_level();level <= o.max_level();++level)
-           o[level].matrix.clear();
-       }
+      clear_object(matrices);
+      clear_object(matrices_in);
+      clear_object(matrices_out);
+      clear_object(flux_matrices_up);
+      clear_object(flux_matrices_down);
     }
 }
 

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.