]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
change function calls and data structures in mesh worker to conform to work stream
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Feb 2010 03:06:44 +0000 (03:06 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Feb 2010 03:06:44 +0000 (03:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@20536 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/include/numerics/mesh_worker_info.h
deal.II/deal.II/include/numerics/mesh_worker_info.templates.h
deal.II/deal.II/include/numerics/mesh_worker_loop.h
deal.II/deal.II/source/numerics/mesh_worker_info.cc
deal.II/examples/step-38/step-38.cc
deal.II/examples/step-39/step-39.cc

index f330e125c7a98646acc5b6e0b5b5551306694cc4..8eecafe1003e50c1b2fec5f93f1f42651f6a07f3 100644 (file)
@@ -995,7 +995,7 @@ namespace MeshWorker
     Functional<number>::assemble(const DoFInfo<dim>& info)
     {
       for (unsigned int i=0;i<results.size();++i)
-       results[i] += info.J[i];
+       results[i] += info.value(i);
     }
     
 
@@ -1007,8 +1007,8 @@ namespace MeshWorker
     {
       for (unsigned int i=0;i<results.size();++i)
        {
-         results[i] += info1.J[i];
-         results[i] += info2.J[i];
+         results[i] += info1.value(i);
+         results[i] += info2.value(i);
        }
     }
 
@@ -1053,13 +1053,13 @@ namespace MeshWorker
     inline void
     CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
     {
-      for (unsigned int i=0;i<info.J.size();++i)
+      for (unsigned int i=0;i<info.n_values();++i)
        {
          if (separate_faces &&
              info.face_number != deal_II_numbers::invalid_unsigned_int)
-           results(1)->block(i)(info.face->user_index()) += info.J[i];
+           results(1)->block(i)(info.face->user_index()) += info.value(i);
          else
-           results(0)->block(i)(info.cell->user_index()) += info.J[i];
+           results(0)->block(i)(info.cell->user_index()) += info.value(i);
        }
     }
     
@@ -1070,19 +1070,19 @@ namespace MeshWorker
     CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
                                    const DoFInfo<dim>& info2)
     {
-      for (unsigned int i=0;i<info1.J.size();++i)
+      for (unsigned int i=0;i<info1.n_values();++i)
        {
          if (separate_faces)
            {
-             const double J = info1.J[i] + info2.J[i];
+             const double J = info1.value(i) + info2.value(i);
              results(1)->block(i)(info1.face->user_index()) += J;
              if (info2.face != info1.face)
                results(1)->block(i)(info2.face->user_index()) += J;
            }
          else
            {
-             results(0)->block(i)(info1.cell->user_index()) += .5*info1.J[i];
-             results(0)->block(i)(info2.cell->user_index()) += .5*info2.J[i];
+             results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i);
+             results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
            }
        }
     }
@@ -1114,8 +1114,8 @@ namespace MeshWorker
     ResidualSimple<VECTOR>::assemble(const DoFInfo<dim>& info)
     {
       for (unsigned int k=0;k<residuals.size();++k)
-       for (unsigned int i=0;i<info.R[k].block(0).size();++i)
-         (*residuals(k))(info.indices[i]) += info.R[k].block(0)(i);
+       for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
+         (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
     }
 
     
@@ -1127,10 +1127,10 @@ namespace MeshWorker
     {
       for (unsigned int k=0;k<residuals.size();++k)
        {
-         for (unsigned int i=0;i<info1.R[k].block(0).size();++i)
-           (*residuals(k))(info1.indices[i]) += info1.R[k].block(0)(i);
-         for (unsigned int i=0;i<info2.R[k].block(0).size();++i)
-           (*residuals(k))(info2.indices[i]) += info2.R[k].block(0)(i);
+         for (unsigned int i=0;i<info1.vector(k).block(0).size();++i)
+           (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
+         for (unsigned int i=0;i<info2.vector(k).block(0).size();++i)
+           (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
        }
     }
 
@@ -1187,7 +1187,7 @@ namespace MeshWorker
       const DoFInfo<dim>& info)
     {
       for (unsigned int i=0;i<residuals.size();++i)
-       assemble(*residuals(i), info.R[i], info.indices);
+       assemble(*residuals(i), info.vector(i), info.indices);
     }
 
     
@@ -1200,8 +1200,8 @@ namespace MeshWorker
     {
       for (unsigned int i=0;i<residuals.size();++i)
        {
-         assemble(*residuals(i), info1.R[i], info1.indices);
-         assemble(*residuals(i), info2.R[i], info2.indices);
+         assemble(*residuals(i), info1.vector(i), info1.indices);
+         assemble(*residuals(i), info2.vector(i), info2.indices);
        }
     }
 
@@ -1256,7 +1256,7 @@ namespace MeshWorker
     inline void
     MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info)
     {
-      assemble(info.M1[0].matrix, info.indices, info.indices);
+      assemble(info.matrix(0,false).matrix, info.indices, info.indices);
     }
     
 
@@ -1266,10 +1266,10 @@ namespace MeshWorker
     MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info1,
                                   const DoFInfo<dim>& info2)
     {
-      assemble(info1.M1[0].matrix, info1.indices, info1.indices);
-      assemble(info1.M2[0].matrix, info1.indices, info2.indices);
-      assemble(info2.M1[0].matrix, info2.indices, info2.indices);
-      assemble(info2.M2[0].matrix, info2.indices, info1.indices);
+      assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices);
+      assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices);
+      assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices);
+      assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices);
     }
     
     
@@ -1354,7 +1354,7 @@ namespace MeshWorker
     MGMatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info)
     {
       const unsigned int level = info.cell->level();
-      assemble((*matrix)[level], info.M1[0].matrix, info.indices, info.indices);
+      assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices);
     }
     
 
@@ -1369,10 +1369,10 @@ namespace MeshWorker
       
       if (level1 == level2)
        {
-         assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices);
-         assemble((*matrix)[level1], info1.M2[0].matrix, info1.indices, info2.indices);
-         assemble((*matrix)[level1], info2.M1[0].matrix, info2.indices, info2.indices);
-         assemble((*matrix)[level1], info2.M2[0].matrix, info2.indices, info1.indices);
+         assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+         assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
+         assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
+         assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
        }
       else
        {
@@ -1380,9 +1380,9 @@ namespace MeshWorker
                                           // Do not add info2.M1,
                                           // which is done by
                                           // the coarser cell
-         assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices);
-         assemble_transpose((*flux_up)[level1],info1.M2[0].matrix, info2.indices, info1.indices);
-         assemble((*flux_down)[level1], info2.M2[0].matrix, info2.indices, info1.indices);
+         assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+         assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices);
+         assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
        }
     }
     
@@ -1474,7 +1474,7 @@ namespace MeshWorker
          const unsigned int row = matrices[i]->row;
          const unsigned int col = matrices[i]->column;
 
-         assemble(matrices[i]->matrix, info.M1[i].matrix, row, col, info.indices, info.indices);
+         assemble(matrices[i]->matrix, info.matrix(i,false).matrix, row, col, info.indices, info.indices);
        }
     }
 
@@ -1493,10 +1493,10 @@ namespace MeshWorker
          const unsigned int row = matrices[i]->row;
          const unsigned int col = matrices[i]->column;
 
-         assemble(matrices[i]->matrix, info1.M1[i].matrix, row, col, info1.indices, info1.indices);
-         assemble(matrices[i]->matrix, info1.M2[i].matrix, row, col, info1.indices, info2.indices);
-         assemble(matrices[i]->matrix, info2.M1[i].matrix, row, col, info2.indices, info2.indices);
-         assemble(matrices[i]->matrix, info2.M2[i].matrix, row, col, info2.indices, info1.indices);
+         assemble(matrices[i]->matrix, info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices);
+         assemble(matrices[i]->matrix, info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices);
+         assemble(matrices[i]->matrix, info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices);
+         assemble(matrices[i]->matrix, info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices);
        }
     }
 
@@ -1607,7 +1607,7 @@ namespace MeshWorker
          const unsigned int row = matrices[i]->row;
          const unsigned int col = matrices[i]->column;
 
-         assemble(matrices[i]->matrix[level], info.M1[i].matrix, row, col,
+         assemble(matrices[i]->matrix[level], info.matrix(i,false).matrix, row, col,
                   info.indices, info.indices, level, level);
        }
     }
@@ -1632,10 +1632,10 @@ namespace MeshWorker
 
          if (level1 == level2)
            {
-             assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col, info1.indices, info1.indices, level1, level1);
-             assemble(matrices[i]->matrix[level1], info1.M2[i].matrix, row, col, info1.indices, info2.indices, level1, level2);
-             assemble(matrices[i]->matrix[level1], info2.M1[i].matrix, row, col, info2.indices, info2.indices, level2, level2);
-             assemble(matrices[i]->matrix[level1], info2.M2[i].matrix, row, col, info2.indices, info1.indices, level2, level1);
+             assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices, level1, level1);
+             assemble(matrices[i]->matrix[level1], info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices, level1, level2);
+             assemble(matrices[i]->matrix[level1], info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices, level2, level2);
+             assemble(matrices[i]->matrix[level1], info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices, level2, level1);
            }
          else
            {
@@ -1645,11 +1645,11 @@ namespace MeshWorker
                                                   // Do not add M22,
                                                   // which is done by
                                                   // the coarser cell
-                 assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col,
+                 assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col,
                           info1.indices, info1.indices, level1, level1);             
-                 assemble(flux_up[i]->matrix[level1], info1.M2[i].matrix, row, col,
+                 assemble(flux_up[i]->matrix[level1], info1.matrix(i,true).matrix, row, col,
                           info1.indices, info2.indices, level1, level2, true);
-                 assemble(flux_down[i]->matrix[level1], info2.M2[i].matrix, row, col,
+                 assemble(flux_down[i]->matrix[level1], info2.matrix(i,true).matrix, row, col,
                           info2.indices, info1.indices, level2, level1);
                }
            }
index 8028b762f477757a9adb0a72904ad6a389fc3e6e..96219c8ae798a10fdd6d188b2832ae00dddd7337 100644 (file)
@@ -96,6 +96,66 @@ namespace MeshWorker
                                        */
       void reinit(const BlockIndices& local_sizes);
 
+                                      /**
+                                       * The number of scalar values.
+                                       */
+      unsigned int n_values () const;
+      
+                                      /**
+                                       * The number of vectors.
+                                       */
+      unsigned int n_vectors () const;
+      
+                                      /**
+                                       * The number of matrices.
+                                       */
+      unsigned int n_matrices () const;
+      
+                                      /**
+                                       * Access scalar value at index
+                                       * @p i.
+                                       */
+      number& value(unsigned int i);
+
+                                      /**
+                                       * Read scalar value at index
+                                       * @p i.
+                                       */
+      number value(unsigned int i) const;
+
+                                      /**
+                                       * Access vector at index @p i.
+                                       */
+      BlockVector<number>& vector(unsigned int i);
+      
+                                      /**
+                                       * Read vector at index @p i.
+                                       */
+      const BlockVector<number>& vector(unsigned int i) const;
+      
+                                      /**
+                                       * Access matrix at index @p
+                                       * i. For results on internal
+                                       * faces, a true value for @p
+                                       * external refers to the flux
+                                       * between cells, while false
+                                       * refers to entries coupling
+                                       * inside the cell.
+                                       */
+      MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false);
+      
+                                      /**
+                                       * Read matrix at index @p
+                                       * i. For results on internal
+                                       * faces, a true value for @p
+                                       * external refers to the flux
+                                       * between cells, while false
+                                       * refers to entries coupling
+                                       * inside the cell.
+                                       */
+      const MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false) const;
+      
+    private:
                                       /**
                                        * The local numbers,
                                        * computed on a cell or on a
@@ -135,8 +195,8 @@ namespace MeshWorker
 
 
 /**
- * Basic info class only containing information on geometry and
- * degrees of freedom of the mesh object.
+ * A class containing information on geometry and degrees of freedom
+ * of a mesh object.
  *
  * The information in these objects is usually used by one of the
  * Assembler classes. It is also the kind of information which is
@@ -145,15 +205,9 @@ namespace MeshWorker
  *
  * In addition to the information on degrees of freedom stored in this
  * class, it also provides the local computation space for the worker
- * object operating on it. This space is provided by the base class
- * template DATATYPE. This base class will automatically
+ * object operating on it in LocalResults. This base class will automatically
  * reinitialized on each cell, but initial setup is up to the user and
- * should be done when initialize() for this class is called. The
- * currently available base classes are
- * <ul>
- * <li> LocalVectors
- * <li> LocalMatrices
- * </ul>
+ * should be done when initialize() for this class is called.
  *
  * This class operates in two different modes, corresponding to the
  * data models discussed in the Assembler namespace documentation.
@@ -163,7 +217,9 @@ namespace MeshWorker
  * BlockInfo::initialize_local(). If this function has been used, or
  * the vector has been changed from zero-length, then local dof
  * indices stored in this object will automatically be renumbered to
- * reflect local block structure.
+ * reflect local block structure. This means, the first entries in
+ * #indices will refer to the first block of the system, then comes
+ * the second block and so on.
  *
  * The BlockInfo object is stored as a pointer. Therefore, if the
  * block structure changes, for instance because of mesh refinement,
@@ -258,7 +314,7 @@ namespace MeshWorker
       SmartPointer<const BlockInfo,DoFInfo<dim,spacedim> > block_info;
 
     private:
-                                      /// Fill index vector
+                                      /// Fill index vector with active indices
       void get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator& c);
 
                                       /// Fill index vector with level indices
@@ -329,24 +385,16 @@ namespace MeshWorker
  * @author Guido Kanschat, 2009
  */
   template<int dim, class FEVALUESBASE, int spacedim = dim>
-  class IntegrationInfo : public DoFInfo<dim, spacedim>
+  class IntegrationInfo
   {
     private:
                                       /// vector of FEValues objects
       std::vector<boost::shared_ptr<FEVALUESBASE> > fevalv;
     public:
                                       /**
-                                       * Constructor forwarding
-                                       * information to DoFInfo.
+                                       * Constructor.
                                        */
-      IntegrationInfo(const BlockInfo& block_info);
-
-                                      /**
-                                       * Constructor forwarding
-                                       * information to DoFInfo.
-                                       */
-      template <class DH>
-      IntegrationInfo(const DH& dof_handler);
+      IntegrationInfo();
 
                                       /**
                                        * Build all internal
@@ -376,7 +424,8 @@ namespace MeshWorker
       void initialize(const FiniteElement<dim,spacedim>& el,
                      const Mapping<dim,spacedim>& mapping,
                      const Quadrature<FEVALUES::integral_dimension>& quadrature,
-                     const UpdateFlags flags);
+                     const UpdateFlags flags,
+                     const BlockInfo* local_block_info = 0);
 
                                       /**
                                        * Initialize the data
@@ -467,41 +516,16 @@ namespace MeshWorker
                                        * Reinitialize internal data
                                        * structures for use on a cell.
                                        */
-      template <class DHCellIterator>
-      void reinit(const DHCellIterator& c);
-
+      void reinit(const DoFInfo<dim, spacedim>& i);
+      
                                       /**
-                                       * Reinitialize internal data
-                                       * structures for use on a face.
-                                       */
-      template <class DHCellIterator, class DHFaceIterator>
-      void reinit(const DHCellIterator& c,
-                 const DHFaceIterator& f,
-                 const unsigned int fn);
-
-                                      /**
-                                       * Reinitialize internal data
-                                       * structures for use on a subface.
-                                       */
-      template <class DHCellIterator, class DHFaceIterator>
-      void reinit(const DHCellIterator& c,
-                 const DHFaceIterator& f,
-                 const unsigned int fn,
-                 const unsigned int sn);
-
-
-                                      /**
-                                       * @deprecated This is the
-                                       * old version not using
-                                       * VectorSelector.
-                                       *
                                        * Use the finite element
                                        * functions in #global_data
                                        * and fill the vectors
                                        * #values, #gradients and
                                        * #hessians.
                                        */
-      void fill_local_data(const bool split_fevalues);
+      void fill_local_data(const DoFInfo<dim, spacedim>& info, const bool split_fevalues);
 
                                       /**
                                        * The global data vector
@@ -554,15 +578,13 @@ namespace MeshWorker
       template <typename T>
       IntegrationInfoBox(const T&);
 
-      template <class WORKER, class ASSEMBLER>
+      template <class WORKER>
       void initialize(const WORKER&,
-                     ASSEMBLER &assembler,
                      const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping);
 
-      template <class WORKER, class ASSEMBLER, typename VECTOR>
+      template <class WORKER, typename VECTOR>
       void initialize(const WORKER&,
-                     ASSEMBLER &assembler,
                      const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping,
                      const NamedData<VECTOR*>& data);
@@ -573,6 +595,7 @@ namespace MeshWorker
       boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > boundary_data;
       boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > face_data;
 
+      DoFInfo<dim, spacedim> dof_info;
       CellInfo cell_info;
       FaceInfo boundary_info;
       FaceInfo face_info;
@@ -663,6 +686,101 @@ namespace MeshWorker
   }
 
 
+  template <typename number>
+  inline
+  unsigned int
+  LocalResults<number>::n_values() const
+  {
+    return J.size();
+  }
+  
+
+  template <typename number>
+  inline
+  unsigned int
+  LocalResults<number>::n_vectors() const
+  {
+    return R.size();
+  }
+  
+
+  template <typename number>
+  inline
+  unsigned int
+  LocalResults<number>::n_matrices() const
+  {
+    return M1.size();
+  }
+  
+
+  template <typename number>
+  inline
+  number&
+  LocalResults<number>::value(unsigned int i)
+  {
+    AssertIndexRange(i,J.size());
+    return J[i];
+  }
+  
+
+  template <typename number>
+  inline
+  BlockVector<number>&
+  LocalResults<number>::vector(unsigned int i)
+  {
+    AssertIndexRange(i,R.size());
+    return R[i];
+  }
+  
+  
+  template <typename number>
+  inline
+  MatrixBlock<FullMatrix<number> >&
+  LocalResults<number>::matrix(unsigned int i, bool external)
+  {
+    if (external)
+      {
+       AssertIndexRange(i,M2.size());
+       return M2[i];
+      }
+    AssertIndexRange(i,M1.size());
+    return M1[i];
+  }
+  
+  template <typename number>
+  inline
+  number
+  LocalResults<number>::value(unsigned int i) const
+  {
+    AssertIndexRange(i,J.size());
+    return J[i];
+  }
+  
+
+  template <typename number>
+  inline
+  const BlockVector<number>&
+  LocalResults<number>::vector(unsigned int i) const
+  {
+    AssertIndexRange(i,R.size());
+    return R[i];
+  }
+  
+  
+  template <typename number>
+  inline
+  const MatrixBlock<FullMatrix<number> >&
+  LocalResults<number>::matrix(unsigned int i, bool external) const
+  {
+    if (external)
+      {
+       AssertIndexRange(i,M2.size());
+       return M2[i];
+      }
+    AssertIndexRange(i,M1.size());
+    return M1[i];
+  }
+  
 //----------------------------------------------------------------------//
 
   template <int dim, int spacedim>
@@ -746,17 +864,6 @@ namespace MeshWorker
 
 //----------------------------------------------------------------------//
 
-  template <int dim, class FVB, int spacedim>
-  template <class DH>
-  IntegrationInfo<dim,FVB,spacedim>::IntegrationInfo(const DH& dof_handler)
-                 :
-                 DoFInfo<dim, spacedim>(dof_handler),
-                 fevalv(0),
-                 multigrid(false),
-                 global_data(boost::shared_ptr<VectorDataBase<dim, spacedim> >(new VectorDataBase<dim, spacedim>))
-  {}
-
-
   template <int dim, class FVB, int spacedim>
   inline const FVB&
   IntegrationInfo<dim,FVB,spacedim>::fe_values() const
@@ -776,94 +883,55 @@ namespace MeshWorker
 
 
   template <int dim, class FVB, int spacedim>
-  template <class DHCellIterator>
-  inline void
-  IntegrationInfo<dim,FVB,spacedim>::reinit(const DHCellIterator& c)
-  {
-    DoFInfo<dim,spacedim>::reinit(c);
-    for (unsigned int i=0;i<fevalv.size();++i)
-      {
-       FVB& febase = *fevalv[i];
-       FEValues<dim>& fe = dynamic_cast<FEValues<dim>&> (febase);
-       fe.reinit(this->cell);
-      }
-
-    const bool split_fevalues = this->block_info != 0;
-    fill_local_data(split_fevalues);
-  }
-
-
-  template <int dim, class FVB, int spacedim>
-  template <class DHCellIterator, class DHFaceIterator>
   inline void
-  IntegrationInfo<dim,FVB,spacedim>::reinit(
-    const DHCellIterator& c,
-    const DHFaceIterator& f,
-    const unsigned int fn)
+  IntegrationInfo<dim,FVB,spacedim>::reinit(const DoFInfo<dim, spacedim>& info)
   {
-    DoFInfo<dim,spacedim>::reinit(c, f, fn);
     for (unsigned int i=0;i<fevalv.size();++i)
       {
        FVB& febase = *fevalv[i];
-       FEFaceValues<dim>& fe = dynamic_cast<FEFaceValues<dim>&> (febase);
-       fe.reinit(this->cell, fn);
+       if (info.sub_number != deal_II_numbers::invalid_unsigned_int)
+         {
+                                            // This is a subface
+           FESubfaceValues<dim>& fe = dynamic_cast<FESubfaceValues<dim>&> (febase);
+           fe.reinit(info.cell, info.face_number, info.sub_number);
+         }
+       else if (info.face_number != deal_II_numbers::invalid_unsigned_int)
+         {
+                                            // This is a face
+           FEFaceValues<dim>& fe = dynamic_cast<FEFaceValues<dim>&> (febase);
+           fe.reinit(info.cell, info.face_number);    
+         }
+       else
+         {
+                                            // This is a cell
+           FEValues<dim>& fe = dynamic_cast<FEValues<dim>&> (febase);
+           fe.reinit(info.cell);
+         }
       }
 
-    const bool split_fevalues = this->block_info != 0;
-    fill_local_data(split_fevalues);
+    const bool split_fevalues = info.block_info != 0;
+    fill_local_data(info, split_fevalues);
   }
 
 
-  template <int dim, class FVB, int spacedim>
-  template <class DHCellIterator, class DHFaceIterator>
-  inline void
-  IntegrationInfo<dim,FVB,spacedim>::reinit(
-    const DHCellIterator& c,
-    const DHFaceIterator& f,
-    const unsigned int fn,
-    const unsigned int sn)
-  {
-    DoFInfo<dim,spacedim>::reinit(c, f, fn, sn);
-    for (unsigned int i=0;i<fevalv.size();++i)
-      {
-       FVB& febase = *fevalv[i];
-       FESubfaceValues<dim>& fe = dynamic_cast<FESubfaceValues<dim>&> (febase);
-       fe.reinit(this->cell, fn, sn);
-      }
-
-    const bool split_fevalues = this->block_info != 0;
-    fill_local_data(split_fevalues);
-  }
-
 //----------------------------------------------------------------------//
 
   template <int dim, int sdim>
   template <typename T>
   IntegrationInfoBox<dim,sdim>::IntegrationInfoBox(const T& t)
                  :
-                 cell_info(t),
-                 boundary_info(t),
-                 face_info(t),
-                 subface_info(t),
-                 neighbor_info(t)
+                 dof_info(t)
   {}
 
 
   template <int dim, int sdim>
-  template <class WORKER, class ASSEMBLER>
+  template <class WORKER>
   void
   IntegrationInfoBox<dim,sdim>::
   initialize(const WORKER& integrator,
-            ASSEMBLER &assembler,
             const FiniteElement<dim,sdim>& el,
             const Mapping<dim,sdim>& mapping)
   {
-    assembler.initialize_info(cell_info, false);
-    assembler.initialize_info(boundary_info, false);
-    assembler.initialize_info(face_info, true);
-    assembler.initialize_info(subface_info, true);
-    assembler.initialize_info(neighbor_info, true);
-
     cell_info.initialize<FEValues<dim,sdim> >(el, mapping, integrator.cell_quadrature,
                         integrator.cell_flags);
     boundary_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.boundary_quadrature,
@@ -878,16 +946,15 @@ namespace MeshWorker
 
 
   template <int dim, int sdim>
-  template <class WORKER, class ASSEMBLER, typename VECTOR>
+  template <class WORKER, typename VECTOR>
   void
   IntegrationInfoBox<dim,sdim>::initialize(
     const WORKER& integrator,
-    ASSEMBLER &assembler,
     const FiniteElement<dim,sdim>& el,
     const Mapping<dim,sdim>& mapping,
     const NamedData<VECTOR*>& data)
   {
-    initialize(integrator, assembler, el, mapping);
+    initialize(integrator, el, mapping);
     boost::shared_ptr<VectorData<VECTOR, dim, sdim> > p;
 
     p = boost::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (integrator.cell_selector));
index 3061cff4615ee1542e1f24dedb34d214371fa479..82baecbaa82c65ac45ab5f3832d002856ed529c9 100644 (file)
@@ -63,9 +63,8 @@ namespace MeshWorker
 //----------------------------------------------------------------------//
 
   template<int dim, class FVB, int sdim>
-  IntegrationInfo<dim,FVB,sdim>::IntegrationInfo(const BlockInfo& block_info)
+  IntegrationInfo<dim,FVB,sdim>::IntegrationInfo()
                  :
-                 DoFInfo<dim,sdim>(block_info),
                  fevalv(0),
                  multigrid(false),
                  global_data(boost::shared_ptr<VectorDataBase<dim, sdim> >(new VectorDataBase<dim, sdim>))
@@ -79,9 +78,10 @@ namespace MeshWorker
     const FiniteElement<dim,sdim>& el,
     const Mapping<dim,sdim>& mapping,
     const Quadrature<FEVALUES::integral_dimension>& quadrature,
-    const UpdateFlags flags)
+    const UpdateFlags flags,
+    const BlockInfo* block_info)
   {
-    if (this->block_info == 0 || this->block_info->local().size() == 0)
+    if (block_info == 0 || block_info->local().size() == 0)
       {
        fevalv.resize(1);             
        fevalv[0] = boost::shared_ptr<FVB> (
@@ -165,21 +165,21 @@ namespace MeshWorker
 
   template<int dim, class FVB, int sdim>
   void
-  IntegrationInfo<dim,FVB,sdim>::fill_local_data(bool split_fevalues)
+  IntegrationInfo<dim,FVB,sdim>::fill_local_data(const DoFInfo<dim, sdim>& info, bool split_fevalues)
   {
      if (split_fevalues)
        {
        unsigned int comp = 0;
                                         // Loop over all blocks
-       for (unsigned int b=0;b<this->block_info->local().size();++b)
+       for (unsigned int b=0;b<info.block_info->local().size();++b)
          {
-           const unsigned int fe_no = this->block_info->base_element(b);
+           const unsigned int fe_no = info.block_info->base_element(b);
            const FEValuesBase<dim>& fe = this->fe_values(fe_no);
            const unsigned int n_comp = fe.get_fe().n_components();
-           const unsigned int block_start = this->block_info->local().block_start(b);
-           const unsigned int block_size = this->block_info->local().block_size(b);
+           const unsigned int block_start = info.block_info->local().block_start(b);
+           const unsigned int block_size = info.block_info->local().block_size(b);
            
-           this->global_data->fill(values, gradients, hessians, fe, this->indices,
+           this->global_data->fill(values, gradients, hessians, fe, info.indices,
                                    comp, n_comp, block_start, block_size);
            comp += n_comp;
          }
@@ -188,8 +188,8 @@ namespace MeshWorker
        {
         const FEValuesBase<dim>& fe = this->fe_values(0);
         const unsigned int n_comp = fe.get_fe().n_components();
-        this->global_data->fill(values, gradients, hessians, fe, this->indices,
-                                0, n_comp, 0, this->indices.size());
+        this->global_data->fill(values, gradients, hessians, fe, info.indices,
+                                0, n_comp, 0, info.indices.size());
        }
   }
 }
index 140ac08a11f26b4177940849c0596a7cfad5c211..b48f92cc374b069244a057ac12479f8c1cb404ce 100644 (file)
@@ -66,14 +66,15 @@ namespace MeshWorker
  * @ingroup MeshWorker
  * @author Guido Kanschat, 2009
  */
-  template<class CELLINFO, class FACEINFO, class ITERATOR, typename ASSEMBLER>
+  template<class DOFINFO, class CELLINFO, class FACEINFO, class ITERATOR, typename ASSEMBLER>
   void loop(ITERATOR begin,
            typename identity<ITERATOR>::type end,
+           DOFINFO cell_dof_info,
            CELLINFO& cell_info, FACEINFO& boundary_info,
            FACEINFO& face_info, FACEINFO& subface_info, FACEINFO& neighbor_info,
-           const std_cxx1x::function<void (CELLINFO &)> &cell_worker,
-           const std_cxx1x::function<void (FACEINFO &)> &boundary_worker,
-           const std_cxx1x::function<void (FACEINFO &, FACEINFO &)> &face_worker,
+           const std_cxx1x::function<void (DOFINFO&, CELLINFO &)> &cell_worker,
+           const std_cxx1x::function<void (DOFINFO&, FACEINFO &)> &boundary_worker,
+           const std_cxx1x::function<void (DOFINFO&, DOFINFO&, FACEINFO &, FACEINFO &)> &face_worker,
            ASSEMBLER &assembler,
            bool cells_first = true)
   {
@@ -81,6 +82,13 @@ namespace MeshWorker
     const bool integrate_boundary      = (boundary_worker != 0);
     const bool integrate_interior_face = (face_worker != 0);
 
+    ;
+    DOFINFO face_dof_info = cell_dof_info;
+    DOFINFO neighbor_dof_info = cell_dof_info;
+    assembler.initialize_info(cell_dof_info, false);
+    assembler.initialize_info(face_dof_info, true);
+    assembler.initialize_info(neighbor_dof_info, true);
+    
                                     // Loop over all cells
     for (ITERATOR cell = begin; cell != end; ++cell)
       {
@@ -89,9 +97,10 @@ namespace MeshWorker
                                         // before faces
        if (integrate_cell && cells_first)
          {
-           cell_info.reinit(cell);
-           cell_worker(cell_info);
-           assembler.assemble (cell_info);
+           cell_dof_info.reinit(cell);
+           cell_info.reinit(cell_dof_info);
+           cell_worker(cell_dof_info, cell_info);
+           assembler.assemble(cell_dof_info);
          }
 
        if (integrate_interior_face || integrate_boundary)
@@ -102,9 +111,10 @@ namespace MeshWorker
                {
                  if (integrate_boundary)
                    {
-                     boundary_info.reinit(cell, face, face_no);
-                     boundary_worker(boundary_info);
-                     assembler.assemble (boundary_info);
+                     cell_dof_info.reinit(cell, face, face_no);
+                     boundary_info.reinit(cell_dof_info);
+                     boundary_worker(cell_dof_info, boundary_info);
+                     assembler.assemble(cell_dof_info);
                    }
                }
              else if (integrate_interior_face)
@@ -132,14 +142,18 @@ namespace MeshWorker
                        = cell->neighbor_of_coarser_neighbor(face_no);
                      typename ITERATOR::AccessorType::Container::face_iterator nface
                        = neighbor->face(neighbor_face_no.first);
-                     face_info.reinit(cell, face, face_no);
-                     subface_info.reinit(neighbor, nface, neighbor_face_no.first, neighbor_face_no.second);
+                     
+                     face_dof_info.reinit(cell, face, face_no);
+                     face_info.reinit(face_dof_info);
+                     neighbor_dof_info.reinit(neighbor, nface,
+                                              neighbor_face_no.first, neighbor_face_no.second);
+                     subface_info.reinit(neighbor_dof_info);
                                                       // Neighbor
                                                       // first to
                                                       // conform to
                                                       // old version
-                     face_worker(face_info, subface_info);
-                     assembler.assemble (face_info, subface_info);
+                     face_worker(face_dof_info, neighbor_dof_info, face_info, subface_info);
+                     assembler.assemble (face_dof_info, neighbor_dof_info);
                    }
                  else
                    {
@@ -162,12 +176,13 @@ namespace MeshWorker
                      unsigned int neighbor_face_no = cell->neighbor_of_neighbor(face_no);
                      Assert (neighbor->face(neighbor_face_no) == face, ExcInternalError());
                                                       // Regular interior face
-                     face_info.reinit(cell, face, face_no);
-                     neighbor_info.reinit(neighbor, neighbor->face(neighbor_face_no),
-                                          neighbor_face_no);
-
-                     face_worker(face_info, neighbor_info);
-                     assembler.assemble (face_info, neighbor_info);
+                     face_dof_info.reinit(cell, face, face_no);
+                     face_info.reinit(face_dof_info);
+                     neighbor_dof_info.reinit(neighbor, neighbor->face(neighbor_face_no),
+                                              neighbor_face_no);
+                     neighbor_info.reinit(neighbor_dof_info);
+                     face_worker(face_dof_info, neighbor_dof_info, face_info, neighbor_info);
+                     assembler.assemble(face_dof_info, neighbor_dof_info);
                    }
                }
            } // faces
@@ -175,9 +190,10 @@ namespace MeshWorker
                                         // have to be handled first
        if (integrate_cell && !cells_first)
          {
-           cell_info.reinit(cell);
-           cell_worker(cell_info);
-           assembler.assemble (cell_info);
+           cell_dof_info.reinit(cell);
+           cell_info.reinit(cell_dof_info);
+           cell_worker(cell_dof_info, cell_info);
+           assembler.assemble (cell_dof_info);
          }
       }
   }
@@ -192,14 +208,15 @@ namespace MeshWorker
   void integration_loop(ITERATOR begin,
                        typename identity<ITERATOR>::type end,
                        IntegrationInfoBox<dim, dim>& box,
-                       const std_cxx1x::function<void (CELLINFO &)> &cell_worker,
-                       const std_cxx1x::function<void (FACEINFO &)> &boundary_worker,
-                       const std_cxx1x::function<void (FACEINFO &, FACEINFO &)> &face_worker,
+                       const std_cxx1x::function<void (DoFInfo<dim>&, CELLINFO &)> &cell_worker,
+                       const std_cxx1x::function<void (DoFInfo<dim>&, FACEINFO &)> &boundary_worker,
+                       const std_cxx1x::function<void (DoFInfo<dim>&, DoFInfo<dim>&, FACEINFO &, FACEINFO &)> &face_worker,
                        ASSEMBLER &assembler,
                        bool cells_first = true)
   {
-    loop<CELLINFO,FACEINFO>
+    loop<DoFInfo<dim>,CELLINFO,FACEINFO>
       (begin, end,
+       box.dof_info,
        box.cell_info, box.boundary_info,
        box.face_info,
        box.subface_info, box.neighbor_info,
index cdd5bb741c32fc29cf6602c16c58eea2eef676dd..9d418fd23bc3cea5c7bf2f4888ea8f2ebd068b22 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2009 by the deal.II authors
+//    Copyright (C) 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
@@ -34,15 +34,15 @@ namespace MeshWorker
   template void IntegrationInfo<deal_II_dimension, FEValuesBase<deal_II_dimension> >
   ::initialize<FEValues<deal_II_dimension> >(
     const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
-    const Quadrature<FEValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+    const Quadrature<FEValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
   template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
   ::initialize<FEFaceValues<deal_II_dimension> >(
     const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
-    const Quadrature<FEFaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+    const Quadrature<FEFaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
   template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
   ::initialize<FESubfaceValues<deal_II_dimension> >(
     const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
-    const Quadrature<FESubfaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+    const Quadrature<FESubfaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
 }
 
 #endif
index ab3f260722b9d4dff7e3d52e7c239b131c458c8a..9d7a1ad814c6c18461f35cb1953ec80712491332 100644 (file)
@@ -39,6 +39,7 @@
                                 // Here come the new include files
                                 // for using the MeshWorker framework:
 #include <numerics/mesh_worker.h>
+#include <numerics/mesh_worker_info.h>
 #include <numerics/mesh_worker_loop.h>
 
 #include <iostream>
@@ -200,9 +201,10 @@ class DGMethod
                                     // types of arguments, but have
                                     // in fact other arguments
                                     // already bound.
-    static void integrate_cell_term (CellInfo& info);
-    static void integrate_boundary_term (FaceInfo& info);
-    static void integrate_face_term (FaceInfo& info1,
+    static void integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info);
+    static void integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info);
+    static void integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
+                                    MeshWorker::DoFInfo<dim>& dinfo2, FaceInfo& info1,
                                     FaceInfo& info2);
 };
 
@@ -378,19 +380,20 @@ void DGMethod<dim>::assemble_system ()
                                   // receives all the stuff we
                                   // created so far.
   MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
-  info_box.initialize(integration_worker, assembler, fe, mapping);
+  info_box.initialize(integration_worker, fe, mapping);
 
                                   // Finally, the integration loop
                                   // over all active cells
                                   // (determined by the first
-                                  // argument, which is an active iterator).
-  MeshWorker::integration_loop<CellInfo,FaceInfo>
+                                  // argument, which is an active
+                                  // iterator).
+  MeshWorker::integration_loop<CellInfo, FaceInfo, dim>
     (dof_handler.begin_active(), dof_handler.end(),
      info_box,
      &DGMethod<dim>::integrate_cell_term,
      &DGMethod<dim>::integrate_boundary_term,
      &DGMethod<dim>::integrate_face_term,
-     assembler);
+     assembler, true);
 }
 
 
@@ -412,7 +415,7 @@ void DGMethod<dim>::assemble_system ()
                                 // added soon).
 
 template <int dim>
-void DGMethod<dim>::integrate_cell_term (CellInfo& info)
+void DGMethod<dim>::integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info)
 {
                                   // First, let us retrieve some of
                                   // the objects used here from
@@ -422,7 +425,7 @@ void DGMethod<dim>::integrate_cell_term (CellInfo& info)
                                   // looks more complicated than
                                   // might seem necessary.
   const FEValuesBase<dim>& fe_v = info.fe_values();
-  FullMatrix<double>& local_matrix = info.M1[0].matrix;
+  FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
 
                                   // With these objects, we continue
@@ -458,11 +461,11 @@ void DGMethod<dim>::integrate_cell_term (CellInfo& info)
                                 // FESubfaceValues, in order to get access to
                                 // normal vectors.
 template <int dim>
-void DGMethod<dim>::integrate_boundary_term (FaceInfo& info)
+void DGMethod<dim>::integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info)
 {
   const FEFaceValuesBase<dim>& fe_v = info.fe_values();
-  FullMatrix<double>& local_matrix = info.M1[0].matrix;
-  Vector<double>& local_vector = info.R[0].block(0);
+  FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
+  Vector<double>& local_vector = dinfo.vector(0).block(0);
 
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
@@ -504,7 +507,9 @@ void DGMethod<dim>::integrate_boundary_term (FaceInfo& info)
                                 // for each cell and two for coupling
                                 // back and forth.
 template <int dim>
-void DGMethod<dim>::integrate_face_term (FaceInfo& info1,
+void DGMethod<dim>::integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
+                                        MeshWorker::DoFInfo<dim>& dinfo2,
+                                        FaceInfo& info1,
                                         FaceInfo& info2)
 {
                                   // For quadrature points, weights,
@@ -531,10 +536,10 @@ void DGMethod<dim>::integrate_face_term (FaceInfo& info1,
                                   // interior couplings of that cell,
                                   // while the second contains the
                                   // couplings between cells.
-  FullMatrix<double>& u1_v1_matrix = info1.M1[0].matrix;
-  FullMatrix<double>& u2_v1_matrix = info1.M2[0].matrix;
-  FullMatrix<double>& u1_v2_matrix = info2.M2[0].matrix;
-  FullMatrix<double>& u2_v2_matrix = info2.M1[0].matrix;
+  FullMatrix<double>& u1_v1_matrix = dinfo1.matrix(0,false).matrix;
+  FullMatrix<double>& u2_v1_matrix = dinfo1.matrix(0,true).matrix;
+  FullMatrix<double>& u1_v2_matrix = dinfo2.matrix(0,true).matrix;
+  FullMatrix<double>& u2_v2_matrix = dinfo2.matrix(0,false).matrix;
 
                                   // Here, following the previous
                                   // functions, we would have the
index c39e0a8b74bd1d0a75e58ed38d5582f2e20003fe..7de4081c58b63122490349979cd1a34808012b5f 100644 (file)
@@ -97,9 +97,13 @@ template <int dim>
 class MatrixIntegrator : public Subscriptor
 {
   public:
-    static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
-    static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
-    static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+    static void cell(MeshWorker::DoFInfo<dim>& dinfo,
+                    typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+    static void bdry(MeshWorker::DoFInfo<dim>& dinfo,
+                    typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+    static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+                    MeshWorker::DoFInfo<dim>& dinfo2,
+                    typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
                     typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
 };
 
@@ -112,10 +116,12 @@ class MatrixIntegrator : public Subscriptor
                                 // degrees of freedom associated
                                 // with the shape functions.
 template <int dim>
-void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
+void MatrixIntegrator<dim>::cell(
+  MeshWorker::DoFInfo<dim>& dinfo,
+  typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
 {
   const FEValuesBase<dim>& fe = info.fe_values();
-  FullMatrix<double>& local_matrix = info.M1[0].matrix;
+  FullMatrix<double>& local_matrix = dinfo.matrix(0,false).matrix;
   
   for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
     for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
@@ -127,13 +133,15 @@ void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::Ce
                                 // On boundary faces, we use the
                                 // Nitsche boundary condition
 template <int dim>
-void MatrixIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void MatrixIntegrator<dim>::bdry(
+  MeshWorker::DoFInfo<dim>& dinfo,
+  typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
 {
   const FEFaceValuesBase<dim>& fe = info.fe_values();
-  FullMatrix<double>& local_matrix = info.M1[0].matrix;
+  FullMatrix<double>& local_matrix = dinfo.matrix(0,false).matrix;
   
   const unsigned int deg = fe.get_fe().tensor_degree();
-  const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+  const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
   
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
     for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
@@ -146,24 +154,27 @@ void MatrixIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::Fa
 
 
 template <int dim>
-void MatrixIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
-                                typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
+void MatrixIntegrator<dim>::face(
+  MeshWorker::DoFInfo<dim>& dinfo1,
+  MeshWorker::DoFInfo<dim>& dinfo2,
+  typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+  typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
 {
   const FEFaceValuesBase<dim>& fe1 = info1.fe_values();
   const FEFaceValuesBase<dim>& fe2 = info2.fe_values();
-  FullMatrix<double>& matrix_v1u1 = info1.M1[0].matrix;
-  FullMatrix<double>& matrix_v1u2 = info1.M2[0].matrix;
-  FullMatrix<double>& matrix_v2u1 = info2.M2[0].matrix;
-  FullMatrix<double>& matrix_v2u2 = info2.M1[0].matrix;
+  FullMatrix<double>& matrix_v1u1 = dinfo1.matrix(0,false).matrix;
+  FullMatrix<double>& matrix_v1u2 = dinfo1.matrix(0,true).matrix;
+  FullMatrix<double>& matrix_v2u1 = dinfo2.matrix(0,true).matrix;
+  FullMatrix<double>& matrix_v2u2 = dinfo2.matrix(0,false).matrix;
   
   const unsigned int deg = fe1.get_fe().tensor_degree();
-  double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
-  double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
-  if (info1.cell->has_children() ^ info2.cell->has_children())
+  double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+  double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
+  if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
     {
-      Assert (info1.face == info2.face, ExcInternalError());
-      Assert (info1.face->has_children(), ExcInternalError());
-//      Assert (info1.cell->has_children(), ExcInternalError());
+      Assert (dinfo1.face == dinfo2.face, ExcInternalError());
+      Assert (dinfo1.face->has_children(), ExcInternalError());
+//      Assert (dinfo1.cell->has_children(), ExcInternalError());
       penalty1 *= 2;
     }
 const double penalty = penalty1 + penalty2;
@@ -196,29 +207,31 @@ template <int dim>
 class RHSIntegrator : public Subscriptor
 {
   public:
-    static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
-    static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
-    static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+    static void cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+    static void bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+    static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+                    MeshWorker::DoFInfo<dim>& dinfo2,
+                    typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
                     typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
 };
 
 
 template <int dim>
-void RHSIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo&)
+void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>&, typename MeshWorker::IntegrationWorker<dim>::CellInfo&)
 {}
 
 
 template <int dim>
-void RHSIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void RHSIntegrator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
 {
   const FEFaceValuesBase<dim>& fe = info.fe_values();
-  Vector<double>& local_vector = info.R[0].block(0);
+  Vector<double>& local_vector = dinfo.vector(0).block(0);
   
   std::vector<double> boundary_values(fe.n_quadrature_points);
   exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
   
   const unsigned int deg = fe.get_fe().tensor_degree();
-  const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+  const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
   
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
     for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
@@ -229,7 +242,9 @@ void RHSIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceI
 
 
 template <int dim>
-void RHSIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo&,
+void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim>&,
+                             MeshWorker::DoFInfo<dim>&,
+                             typename MeshWorker::IntegrationWorker<dim>::FaceInfo&,
                              typename MeshWorker::IntegrationWorker<dim>::FaceInfo&)
 {}
 
@@ -238,29 +253,31 @@ template <int dim>
 class Estimator : public Subscriptor
 {
   public:
-    static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
-    static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
-    static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+    static void cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+    static void bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+    static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+                    MeshWorker::DoFInfo<dim>& dinfo2,
+                    typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
                     typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
 };
 
 
 template <int dim>
-void Estimator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
+void Estimator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
 {
   const FEValuesBase<dim>& fe = info.fe_values();
   
   const std::vector<Tensor<2,dim> >& DDuh = info.hessians[0][0];
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
     {
-      const double t = info.cell->diameter() * trace(DDuh[k]);
-      info.J[0] +=  t*t * fe.JxW(k);
+      const double t = dinfo.cell->diameter() * trace(DDuh[k]);
+      dinfo.value(0) +=  t*t * fe.JxW(k);
     }
 }
 
 
 template <int dim>
-void Estimator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void Estimator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
 {
   const FEFaceValuesBase<dim>& fe = info.fe_values();
   
@@ -270,16 +287,18 @@ void Estimator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo&
   const std::vector<double>& uh = info.values[0][0];
   
   const unsigned int deg = fe.get_fe().tensor_degree();
-  const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+  const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
   
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
-    info.J[0] += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
+    dinfo.value(0) += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
                 * fe.JxW(k);
 }
 
 
 template <int dim>
-void Estimator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+void Estimator<dim>::face(MeshWorker::DoFInfo<dim>& dinfo1,
+                         MeshWorker::DoFInfo<dim>& dinfo2,
+                         typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
                          typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
 {
   const FEFaceValuesBase<dim>& fe = info1.fe_values();
@@ -289,19 +308,19 @@ void Estimator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo&
   const std::vector<Tensor<1,dim> >& Duh2 = info2.gradients[0][0];
   
   const unsigned int deg = fe.get_fe().tensor_degree();
-  const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
-  const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
+  const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+  const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
   const double penalty = penalty1 + penalty2;
-  const double h = info1.face->measure();
+  const double h = dinfo1.face->measure();
   
   for (unsigned k=0;k<fe.n_quadrature_points;++k)
     {
       double diff1 = uh1[k] - uh2[k];
       double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
-      info1.J[0] += (penalty * diff1*diff1 + h * diff2*diff2)
+      dinfo1.value(0) += (penalty * diff1*diff1 + h * diff2*diff2)
                    * fe.JxW(k);
     }
-  info2.J[0] = info1.J[0];
+  dinfo2.value(0) = dinfo1.value(0);
 }
 
 
@@ -425,8 +444,8 @@ Step39<dim>::assemble_matrix()
 
   assembler.initialize(matrix);
   MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
-  info_box.initialize(integration_worker, assembler, fe, mapping);
-  MeshWorker::integration_loop<CellInfo, FaceInfo>(
+  info_box.initialize(integration_worker, fe, mapping);
+  MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
     dof_handler.begin_active(), dof_handler.end(),
     info_box,
     &MatrixIntegrator<dim>::cell,
@@ -451,8 +470,8 @@ Step39<dim>::assemble_mg_matrix()
   assembler.initialize(mg_matrix);
   assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
   MeshWorker::IntegrationInfoBox<dim> info_box(mg_dof_handler);
-  info_box.initialize(integration_worker, assembler, fe, mapping);
-  MeshWorker::integration_loop<CellInfo, FaceInfo>(
+  info_box.initialize(integration_worker, fe, mapping);
+  MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
     mg_dof_handler.begin(), mg_dof_handler.end(),
     info_box,
     &MatrixIntegrator<dim>::cell,
@@ -479,9 +498,9 @@ Step39<dim>::assemble_right_hand_side()
   data.add(rhs, "RHS");
   assembler.initialize(data);
   MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
-  info_box.initialize(integration_worker, assembler, fe, mapping);
+  info_box.initialize(integration_worker, fe, mapping);
   
-  MeshWorker::integration_loop<CellInfo, FaceInfo>(
+  MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
     dof_handler.begin_active(), dof_handler.end(),
     info_box,
     &RHSIntegrator<dim>::cell,
@@ -600,8 +619,8 @@ Step39<dim>::estimate()
   out_data.add(est, "cells");
   assembler.initialize(out_data, false);
   MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
-  info_box.initialize(integration_worker, assembler, fe, mapping, solution_data);
-  MeshWorker::integration_loop<CellInfo, FaceInfo> (
+  info_box.initialize(integration_worker, fe, mapping, solution_data);
+  MeshWorker::integration_loop<CellInfo, FaceInfo, dim> (
     dof_handler.begin_active(), dof_handler.end(),
     info_box,
     &Estimator<dim>::cell,

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.