]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
start splitting assembler.h
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Sep 2011 16:15:40 +0000 (16:15 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Sep 2011 16:15:40 +0000 (16:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@24273 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/meshworker/assembler.h
deal.II/include/deal.II/meshworker/functional.h [new file with mode: 0644]
deal.II/include/deal.II/meshworker/simple.h [new file with mode: 0644]

index 2f968c3f7e003833d9ac7df28ab52059c25978e7..54acf9779f56b027072373723205716705327927 100644 (file)
@@ -18,6 +18,8 @@
 #include <deal.II/base/mg_level_object.h>
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/functional.h>
+#include <deal.II/meshworker/simple.h>
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 
 
@@ -78,229 +80,6 @@ namespace MeshWorker
  */
   namespace Assembler
   {
-/**
- * The class assembling local contributions to a functional into the
- * global functionals.
- *
- * 
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <typename number = double>
-    class Functional
-    {
-      public:
-                                        /**
-                                         * Initialize local data to
-                                         * store functionals. The
-                                         * number <tt>n</tt> is the
-                                         * number of functionals to
-                                         * be computed.
-                                         */
-       void initialize(unsigned int n);
-                                        /**
-                                         * Initialize the local data
-                                         * in the DoFInfo object used
-                                         * later for assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face);
-       
-                                        /**
-                                         * Assemble the local values
-                                         * into the global vectors.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local values
-                                         * into the global vectors.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-
-                                        /**
-                                         * The value of the ith entry
-                                         * in #results.
-                                         */
-       number operator() (unsigned int i) const;
-      private:
-                                        /**
-                                         * The values into which the
-                                         * results are added.
-                                         */
-       std::vector<double> results;
-    };
-
-/**
- * Compute cell and face contributions of one or several functionals,
- * typically for error estimates.
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <typename number = double>
-    class CellsAndFaces
-    {
-      public:
-                                        /**
-                                         * The data type for
-                                         * communicating the cell and
-                                         * face vectors.
-                                         */
-       typedef NamedData<BlockVector<number>*> DataVectors;
-       
-                                        /**
-                                         * The initialization
-                                         * function, specifying the
-                                         * #results vectors and
-                                         * whether face data should
-                                         * be collected separately.
-                                         *
-                                         * #results should contain
-                                         * two block vectors named
-                                         * "cells" and "faces" (the
-                                         * latter only if
-                                         * #separate_faces is
-                                         * true). In each of the two,
-                                         * each block should have
-                                         * equal size and be large
-                                         * enough to accomodate all
-                                         * user indices set in the
-                                         * cells and faces covered by
-                                         * the loop it is used
-                                         * in. Typically, for
-                                         * estimators, this is
-                                         * Triangulation::n_active_cells()
-                                         * and
-                                         * Triangulation::n_faces(),
-                                         * respectively.
-                                         *
-                                         * The use of BlockVector may
-                                         * seem cumbersome, but it
-                                         * allows us to assemble
-                                         * several functionals at the
-                                         * same time, one in each
-                                         * block. The typical
-                                         * situation for error
-                                         * estimate is just having a
-                                         * single block in each vector.
-                                         */
-       void initialize(DataVectors& results,
-                       bool separate_faces = true);
-                                        /**
-                                         * Initialize the local data
-                                         * in the
-                                         * DoFInfo
-                                         * object used later for
-                                         * assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-                                        /**
-                                         * Assemble the local values
-                                         * into the global vectors.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local values
-                                         * into the global vectors.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-
-                                        /**
-                                         * The value of the ith entry
-                                         * in #results.
-                                         */
-       number operator() (unsigned int i) const;
-      private:
-       DataVectors results;
-       bool separate_faces;
-    };
-    
-
-/**
- * Assemble residuals without block structure.
- *
- * The data structure for this Assembler class is a simple vector on
- * each cell with entries from zero to
- * FiniteElementData::dofs_per_cell and a simple global vector with
- * entries numbered from zero to DoFHandler::n_dofs(). No BlockInfo is
- * required and the global vector may be any type of vector having
- * element access through <tt>operator() (unsigned int)</tt>
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <class VECTOR>
-    class ResidualSimple
-    {
-      public:
-       void initialize(NamedData<VECTOR*>& results);
-                                        /**
-                                         * Initialize the constraints. 
-                                         */
-        void initialize(const ConstraintMatrix& constraints);
-                                        /**
-                                         * Initialize the local data
-                                         * in the DoFInfo object used
-                                         * later for assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-                                        /**
-                                         * Assemble the local residuals
-                                         * into the global residuals.
-                                         */
-       template <class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local residuals
-                                         * into the global residuals.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-      private:
-                                        /**
-                                         * The global residal vectors
-                                         * filled by assemble().
-                                         */
-       NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
-      /**
-       * A pointer to the object containing constraints.
-       */
-      SmartPointer<const ConstraintMatrix,ResidualSimple<VECTOR> > constraints;
-    };
-    
 /**
  * Assemble local residuals into global residuals.
  *
@@ -343,373 +122,61 @@ namespace MeshWorker
                                          * DoFInfo
                                          * object used later for
                                          * assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-
-                                        /**
-                                         * Assemble the local residuals
-                                         * into the global residuals.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local residuals
-                                         * into the global residuals.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-      private:
-                                        /**
-                                         * Assemble a single local
-                                         * residual into the global.
-                                         */
-       void assemble(VECTOR& global,
-                     const BlockVector<double>& local,
-                     const std::vector<unsigned int>& dof);
-       
-                                        /**
-                                         * The global matrices,
-                                         * stored as a vector of
-                                         * pointers.
-                                         */
-       NamedData<SmartPointer<VECTOR,
-          ResidualLocalBlocksToGlobalBlocks<VECTOR> > > residuals;
-       
-      /**
-       * A pointer to the object containing the block structure.
-       */
-      SmartPointer<const BlockInfo,
-        ResidualLocalBlocksToGlobalBlocks<VECTOR> > block_info;
-      /**
-       * A pointer to the object containing constraints.
-       */
-      SmartPointer<const ConstraintMatrix,
-        ResidualLocalBlocksToGlobalBlocks<VECTOR> > constraints;
-   };
-
-
-/**
- * Assemble local matrices into a single global matrix without using
- * block structure.
- *
- * After being initialized with a SparseMatrix object (or another
- * matrix offering the same functionality as SparseMatrix::add()),
- * this class can be used in a MeshWorker::loop() to assemble the cell
- * and face matrices into the global matrix.
- *
- * @todo On locally refined meshes, a ConstraintMatrix should be used
- * to automatically eliminate hanging nodes.
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <class MATRIX>
-    class MatrixSimple
-    {
-      public:
-                                        /**
-                                         * Constructor, initializing
-                                         * the #threshold, which
-                                         * limits how small numbers
-                                         * may be to be entered into
-                                         * the matrix.
-                                         */
-       MatrixSimple(double threshold = 1.e-12);
-       
-                                        /**
-                                         * Store the result matrix
-                                         * for later assembling.
-                                         */
-       void initialize(MATRIX& m);
-                                        /**
-                                         * Initialize the constraints. 
-                                         */
-        void initialize(const ConstraintMatrix& constraints);
-                                        /**
-                                         * Initialize the local data
-                                         * in the
-                                         * DoFInfo
-                                         * object used later for
-                                         * assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-                                        /**
-                                         * Assemble the matrix
-                                         * DoFInfo::M1[0]
-                                         * into the global matrix.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local
-                                         * matrices in the info
-                                         * objects into the global
-                                         * matrix.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-      private:
-                                        /**
-                                         * Assemble a single matrix
-                                         * into #matrix.
-                                         */
-       void assemble(const FullMatrix<double>& M,
-                     const std::vector<unsigned int>& i1,
-                     const std::vector<unsigned int>& i2);
-       
-                                        /**
-                                         * The global matrix being
-                                         * assembled.
-                                         */
-       SmartPointer<MATRIX,MatrixSimple<MATRIX> > matrix;
-      /**
-       * A pointer to the object containing constraints.
-       */
-        SmartPointer<const ConstraintMatrix,MatrixSimple<MATRIX> > constraints;
-       
-                                        /**
-                                         * The smallest positive
-                                         * number that will be
-                                         * entered into the global
-                                         * matrix. All smaller
-                                         * absolute values will be
-                                         * treated as zero and will
-                                         * not be assembled.
-                                         */
-       const double threshold;
-       
-    };
-    
-    
-/**
- * Assemble local matrices into level matrices without using
- * block structure.
- *
- * @todo The matrix structures needed for assembling level matrices
- * with local refinement and continuous elements are missing.
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <class MATRIX>
-    class MGMatrixSimple
-    {
-      public:
-                                        /**
-                                         * Constructor, initializing
-                                         * the #threshold, which
-                                         * limits how small numbers
-                                         * may be to be entered into
-                                         * the matrix.
-                                         */
-       MGMatrixSimple(double threshold = 1.e-12);
-       
-                                        /**
-                                         * Store the result matrix
-                                         * for later assembling.
-                                         */
-       void initialize(MGLevelObject<MATRIX>& m);
-
-                                        /**
-                                         * Initialize the multilevel
-                                          * constraints. 
-                                         */
-        void initialize(const MGConstrainedDoFs& mg_constrained_dofs);
-
-                                        /**
-                                         * Initialize the matrices
-                                         * #flux_up and #flux_down
-                                         * used for local refinement
-                                         * with discontinuous
-                                         * Galerkin methods.
-                                         */
-       void initialize_fluxes(MGLevelObject<MATRIX>& flux_up,
-                              MGLevelObject<MATRIX>& flux_down);
-       
-                                        /**
-                                         * Initialize the matrices
-                                         * #interface_in and #interface_out
-                                         * used for local refinement
-                                         * with continuous
-                                         * Galerkin methods.
-                                         */
-
-       void initialize_interfaces(MGLevelObject<MATRIX>& interface_in,
-                              MGLevelObject<MATRIX>& interface_out);
-                                        /**
-                                         * Initialize the local data
-                                         * in the
-                                         * DoFInfo
-                                         * object used later for
-                                         * assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-                                        /**
-                                         * Assemble the matrix
-                                         * DoFInfo::M1[0]
-                                         * into the global matrix.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local
-                                         * matrices in the info
-                                         * objects into the global
-                                         * matrices.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-      private:
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-    void assemble(MATRIX& G,
-                 const FullMatrix<double>& M,
-                 const std::vector<unsigned int>& i1,
-                 const std::vector<unsigned int>& i2);
-
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-    void assemble(MATRIX& G,
-                 const FullMatrix<double>& M,
-                 const std::vector<unsigned int>& i1,
-                 const std::vector<unsigned int>& i2,
-                  const unsigned int level);
-       
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-
-    void assemble_up(MATRIX& G,
-                           const FullMatrix<double>& M,
-                           const std::vector<unsigned int>& i1,
-                           const std::vector<unsigned int>& i2,
-                            const unsigned int level = -1);
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-
-    void assemble_down(MATRIX& G,
-                           const FullMatrix<double>& M,
-                           const std::vector<unsigned int>& i1,
-                           const std::vector<unsigned int>& i2,
-                            const unsigned int level = -1);
-
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-
-    void assemble_in(MATRIX& G,
-                 const FullMatrix<double>& M,
-                 const std::vector<unsigned int>& i1,
-                 const std::vector<unsigned int>& i2,
-                  const unsigned int level = -1);
-
-                                        /**
-                                         * Assemble a single matrix
-                                         * into a global matrix.
-                                         */
-
-    void assemble_out(MATRIX& G,
-                 const FullMatrix<double>& M,
-                 const std::vector<unsigned int>& i1,
-                 const std::vector<unsigned int>& i2,
-                  const unsigned int level = -1);
-       
-                                        /**
-                                         * The global matrix being
-                                         * assembled.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > matrix;
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
        
+
                                         /**
-                                         * The matrix used for face
-                                         * flux terms across the
-                                         * refinement edge, coupling
-                                         * coarse to fine.
+                                         * Assemble the local residuals
+                                         * into the global residuals.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_up;
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
        
                                         /**
-                                         * The matrix used for face
-                                         * flux terms across the
-                                         * refinement edge, coupling
-                                         * fine to coarse.
+                                         * Assemble both local residuals
+                                         * into the global residuals.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_down;
-
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+      private:
                                         /**
-                                         * The matrix used for face
-                                         * contributions for continuous
-                                          * elements across the
-                                         * refinement edge, coupling
-                                         * coarse to fine.
+                                         * Assemble a single local
+                                         * residual into the global.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_in;
+       void assemble(VECTOR& global,
+                     const BlockVector<double>& local,
+                     const std::vector<unsigned int>& dof);
        
                                         /**
-                                         * The matrix used for face
-                                         * contributions for continuous
-                                          * elements across the
-                                         * refinement edge, coupling
-                                         * fine to coarse.
+                                         * The global matrices,
+                                         * stored as a vector of
+                                         * pointers.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_out;
+       NamedData<SmartPointer<VECTOR,
+          ResidualLocalBlocksToGlobalBlocks<VECTOR> > > residuals;
+       
+      /**
+       * A pointer to the object containing the block structure.
+       */
+      SmartPointer<const BlockInfo,
+        ResidualLocalBlocksToGlobalBlocks<VECTOR> > block_info;
       /**
        * A pointer to the object containing constraints.
        */
-        SmartPointer<const MGConstrainedDoFs,MGMatrixSimple<MATRIX> > mg_constrained_dofs;
-       
-                                        /**
-                                         * The smallest positive
-                                         * number that will be
-                                         * entered into the global
-                                         * matrix. All smaller
-                                         * absolute values will be
-                                         * treated as zero and will
-                                         * not be assembled.
-                                         */
-       const double threshold;
-       
-    };
-    
-    
+      SmartPointer<const ConstraintMatrix,
+        ResidualLocalBlocksToGlobalBlocks<VECTOR> > constraints;
+   };
+
+
 /**
  * A helper class assembling local matrices into global matrices.
  *
@@ -947,318 +414,62 @@ namespace MeshWorker
          const std::vector<unsigned int>& dof1,
          const std::vector<unsigned int>& dof2,
          unsigned int level1,
-         unsigned int level2,
-         bool transpose = false);
-       
-                                        /**
-                                         * The level matrices,
-                                         * stored as a vector of
-                                         * pointers.
-                                         */
-       MatrixPtrVectorPtr matrices;
-       
-                                        /**
-                                         * The flux matrix between
-                                         * the fine and the coarse
-                                         * level at refinement edges.
-                                         */
-       MatrixPtrVectorPtr flux_down;
-       
-                                        /**
-                                         * The flux matrix between
-                                         * the coarse and the fine
-                                         * level at refinement edges.
-                                         */
-       MatrixPtrVectorPtr flux_up;     
-
-                                        /**
-                                         * The interface matrix between
-                                         * the fine and the coarse
-                                         * level at refinement edges.
-                                         */
-       MatrixPtrVectorPtr interface_out;
-       
-                                        /**
-                                         * The interface matrix between
-                                         * the coarse and the fine
-                                         * level at refinement edges.
-                                         */
-       MatrixPtrVectorPtr interface_in;        
-      
-      /**
-       * A pointer to the object containing the block structure.
-       */
-      SmartPointer<const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
-      
-                                        /**
-                                         * The smallest positive
-                                         * number that will be
-                                         * entered into the global
-                                         * matrix. All smaller
-                                         * absolute values will be
-                                         * treated as zero and will
-                                         * not be assembled.
-                                         */
-       const double threshold;
-       
-    };
-
-/**
- * Assemble a simple matrix and a simple right hand side at once. We
- * use a combination of MatrixSimple and ResidualSimple to achieve
- * this. Cell and face operators should fill the matrix and vector
- * objects in LocalResults and this class will assemble
- * them into matrix and vector objects.
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
-    template <class MATRIX, class VECTOR>
-    class SystemSimple :
-       private MatrixSimple<MATRIX>,
-       private ResidualSimple<VECTOR>
-    {
-      public:
-                                        /**
-                                         * Constructor setting the
-                                         * threshold value in
-                                         * MatrixSimple.
-                                         */
-       SystemSimple(double threshold = 1.e-12);
-
-                                        /**
-                                         * Store the two objects data
-                                         * is assembled into.
-                                         */
-       void initialize(MATRIX& m, VECTOR& rhs);
-       
-                                        /**
-                                         * Initialize the local data
-                                         * in the
-                                         * DoFInfo
-                                         * object used later for
-                                         * assembling.
-                                         *
-                                         * The info object refers to
-                                         * a cell if
-                                         * <code>!face</code>, or
-                                         * else to an interior or
-                                         * boundary face.
-                                         */
-       template <class DOFINFO>
-       void initialize_info(DOFINFO& info, bool face) const;
-       
-                                        /**
-                                         * Assemble the matrix
-                                         * DoFInfo::M1[0]
-                                         * into the global matrix.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info);
-       
-                                        /**
-                                         * Assemble both local
-                                         * matrices in the info
-                                         * objects into the global
-                                         * matrix.
-                                         */
-       template<class DOFINFO>
-       void assemble(const DOFINFO& info1,
-                     const DOFINFO& info2);
-    };
-
-    
-//----------------------------------------------------------------------//
-    
-    template <typename number>
-    inline void
-    Functional<number>::initialize(unsigned int n)
-    {
-      results.resize(n);
-      std::fill(results.begin(), results.end(), 0.);
-    }
-    
-    
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    Functional<number>::initialize_info(DOFINFO& info, bool)
-    {
-      info.initialize_numbers(results.size());
-    }
-    
-    
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    Functional<number>::assemble(const DOFINFO& info)
-    {
-      for (unsigned int i=0;i<results.size();++i)
-       results[i] += info.value(i);
-    }
-    
-
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    Functional<number>::assemble(const DOFINFO& info1,
-                                const DOFINFO& info2)
-    {
-      for (unsigned int i=0;i<results.size();++i)
-       {
-         results[i] += info1.value(i);
-         results[i] += info2.value(i);
-       }
-    }
-
-    
-    template <typename number>
-    inline number
-    Functional<number>::operator() (unsigned int i) const
-    {
-      AssertIndexRange(i, results.size());
-      return results[i];
-    }
-    
-//----------------------------------------------------------------------//
-
-    template <typename number>
-    inline void
-    CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
-    {
-      Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
-      if (sep)
-       {
-         Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
-         AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
-       }
-      
-      results = r;
-      separate_faces = sep;
-    }
-
-    
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    CellsAndFaces<number>::initialize_info(DOFINFO& info, bool) const
-    {
-      info.initialize_numbers(results(0)->n_blocks());
-    }
-
-
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    CellsAndFaces<number>::assemble(const DOFINFO& info)
-    {
-      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.value(i);
-         else
-           results(0)->block(i)(info.cell->user_index()) += info.value(i);
-       }
-    }
-    
-
-    template <typename number>
-    template <class DOFINFO>
-    inline void
-    CellsAndFaces<number>::assemble(const DOFINFO& info1,
-                                   const DOFINFO& info2)
-    {
-      for (unsigned int i=0;i<info1.n_values();++i)
-       {
-         if (separate_faces)
-           {
-             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.value(i);
-             results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
-           }
-       }
-    }
-    
-
-
-//----------------------------------------------------------------------//
-
-    template <class VECTOR>
-    inline void
-    ResidualSimple<VECTOR>::initialize(NamedData<VECTOR*>& results)
-    {
-      residuals = results;
-    }
-
-    template <class VECTOR>
-    inline void
-    ResidualSimple<VECTOR>::initialize(const ConstraintMatrix& c)
-    {
-      constraints = &c;
-    }
-
-    
-    template <class VECTOR>
-    template <class DOFINFO>
-    inline void
-    ResidualSimple<VECTOR>::initialize_info(DOFINFO& info, bool) const
-    {
-      info.initialize_vectors(residuals.size());
-    }
-
-    
-    template <class VECTOR>
-    template <class DOFINFO>
-    inline void
-    ResidualSimple<VECTOR>::assemble(const DOFINFO& info)
-    {
-      for (unsigned int k=0;k<residuals.size();++k)
-      {
-        if(constraints == 0)
-        {
-          for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
-            (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
-        }
-        else
-          constraints->distribute_local_to_global(
-              info.vector(k).block(0), info.indices, (*residuals(k)));
-      }
-    }
+         unsigned int level2,
+         bool transpose = false);
+       
+                                        /**
+                                         * The level matrices,
+                                         * stored as a vector of
+                                         * pointers.
+                                         */
+       MatrixPtrVectorPtr matrices;
+       
+                                        /**
+                                         * The flux matrix between
+                                         * the fine and the coarse
+                                         * level at refinement edges.
+                                         */
+       MatrixPtrVectorPtr flux_down;
+       
+                                        /**
+                                         * The flux matrix between
+                                         * the coarse and the fine
+                                         * level at refinement edges.
+                                         */
+       MatrixPtrVectorPtr flux_up;     
 
-    
-    template <class VECTOR>
-    template <class DOFINFO>
-    inline void
-    ResidualSimple<VECTOR>::assemble(const DOFINFO& info1,
-                                       const DOFINFO& info2)
-    {
-      for (unsigned int k=0;k<residuals.size();++k)
-       {
-          if(constraints == 0)
-          {
-            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);
-          }
-          else
-          {
-            constraints->distribute_local_to_global(
-              info1.vector(k).block(0), info1.indices, (*residuals(k)));
-            constraints->distribute_local_to_global(
-              info2.vector(k).block(0), info2.indices, (*residuals(k)));
-          }
-       }
-    }
+                                        /**
+                                         * The interface matrix between
+                                         * the fine and the coarse
+                                         * level at refinement edges.
+                                         */
+       MatrixPtrVectorPtr interface_out;
+       
+                                        /**
+                                         * The interface matrix between
+                                         * the coarse and the fine
+                                         * level at refinement edges.
+                                         */
+       MatrixPtrVectorPtr interface_in;        
+      
+      /**
+       * A pointer to the object containing the block structure.
+       */
+      SmartPointer<const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
+      
+                                        /**
+                                         * The smallest positive
+                                         * number that will be
+                                         * entered into the global
+                                         * matrix. All smaller
+                                         * absolute values will be
+                                         * treated as zero and will
+                                         * not be assembled.
+                                         */
+       const double threshold;
+       
+    };
 
-    
 //----------------------------------------------------------------------//
 
     template <class VECTOR>
@@ -1343,421 +554,6 @@ namespace MeshWorker
     }
 
 
-//----------------------------------------------------------------------//
-
-    template <class MATRIX>
-    inline
-    MatrixSimple<MATRIX>::MatrixSimple(double threshold)
-                   :
-                   threshold(threshold)
-    {}
-    
-    
-    template <class MATRIX>
-    inline void
-    MatrixSimple<MATRIX>::initialize(MATRIX& m)
-    {
-      matrix = &m;
-    }
-
-
-    template <class MATRIX>
-    inline void
-    MatrixSimple<MATRIX>::initialize(const ConstraintMatrix& c)
-    {
-      constraints = &c;
-    }
-    
-
-    template <class MATRIX >
-    template <class DOFINFO>
-    inline void
-    MatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
-    {
-      info.initialize_matrices(1, face);
-    }
-
-
-
-    template <class MATRIX>
-    inline void
-    MatrixSimple<MATRIX>::assemble(const FullMatrix<double>& M,
-                                  const std::vector<unsigned int>& i1,
-                                  const std::vector<unsigned int>& i2)
-    {
-      AssertDimension(M.m(), i1.size());
-      AssertDimension(M.n(), i2.size());
-     
-     if(constraints == 0)
-     { 
-      for (unsigned int j=0; j<i1.size(); ++j)
-       for (unsigned int k=0; k<i2.size(); ++k)
-         if (std::fabs(M(j,k)) >= threshold)
-           matrix->add(i1[j], i2[k], M(j,k));
-     }
-     else
-       constraints->distribute_local_to_global(
-           M, i1, i2, *matrix);
-    }
-    
-    
-    template <class MATRIX>
-    template <class DOFINFO>
-    inline void
-    MatrixSimple<MATRIX>::assemble(const DOFINFO& info)
-    {
-      assemble(info.matrix(0,false).matrix, info.indices, info.indices);
-    }
-    
-
-    template <class MATRIX>
-    template <class DOFINFO>
-    inline void
-    MatrixSimple<MATRIX>::assemble(const DOFINFO& info1,
-                                  const DOFINFO& info2)
-    {
-      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);
-    }
-    
-    
-//----------------------------------------------------------------------//
-
-    template <class MATRIX>
-    inline
-    MGMatrixSimple<MATRIX>::MGMatrixSimple(double threshold)
-                   :
-                   threshold(threshold)
-    {}
-    
-    
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::initialize(MGLevelObject<MATRIX>& m)
-    {
-      matrix = &m;
-    }
-    
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::initialize(const MGConstrainedDoFs& c)
-    {
-      mg_constrained_dofs = &c;
-    }
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::initialize_fluxes(
-      MGLevelObject<MATRIX>& up, MGLevelObject<MATRIX>& down)
-    {
-      flux_up = &up;
-      flux_down = &down;
-    }
-    
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::initialize_interfaces(
-      MGLevelObject<MATRIX>& in, MGLevelObject<MATRIX>& out)
-    {
-      interface_in = &in;
-      interface_out = &out;
-    }
-
-
-    template <class MATRIX >
-    template <class DOFINFO>
-    inline void
-    MGMatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
-    {
-      info.initialize_matrices(1, face);
-    }
-
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2)
-    {
-      AssertDimension(M.m(), i1.size());
-      AssertDimension(M.n(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              G.add(i1[j], i2[k], M(j,k));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-            {
-                if(!mg_constrained_dofs->continuity_across_refinement_edges())
-                  G.add(i1[j], i2[k], M(j,k));
-            }
-      }
-    }
-
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2,
-      const unsigned int level)
-    {
-      AssertDimension(M.m(), i1.size());
-      AssertDimension(M.n(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              G.add(i1[j], i2[k], M(j,k));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-              {
-                if (mg_constrained_dofs->set_boundary_values())
-                {
-                  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                        !mg_constrained_dofs->is_boundary_index(level, i2[k]))
-                      ||
-                      (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                       mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
-                       i1[j] == i2[k]))
-                    G.add(i1[j], i2[k], M(j,k));
-                }
-                else
-                    G.add(i1[j], i2[k], M(j,k));
-            }
-      }
-    }
-    
-    
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble_up(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2,
-      const unsigned int level)
-    {
-      AssertDimension(M.n(), i1.size());
-      AssertDimension(M.m(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(k,j)) >= threshold)
-              G.add(i1[j], i2[k], M(k,j));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(k,j)) >= threshold)
-              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                G.add(i1[j], i2[k], M(k,j));
-      }
-    }
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble_down(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2,
-      const unsigned int level)
-    {
-      AssertDimension(M.m(), i1.size());
-      AssertDimension(M.n(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              G.add(i1[j], i2[k], M(j,k));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                  G.add(i1[j], i2[k], M(j,k));
-      }
-    }
-
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble_in(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2,
-      const unsigned int level)
-    {
-      AssertDimension(M.m(), i1.size());
-      AssertDimension(M.n(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              G.add(i1[j], i2[k], M(j,k));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-              {
-                if (mg_constrained_dofs->set_boundary_values())
-                {
-                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
-                      ||
-                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                       mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                       i1[j] == i2[k]))
-                    G.add(i1[j], i2[k], M(j,k));
-                }
-                else
-                  G.add(i1[j], i2[k], M(j,k));
-              }
-      }
-    }
-    
-    template <class MATRIX>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble_out(
-      MATRIX& G,
-      const FullMatrix<double>& M,
-      const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2,
-      const unsigned int level)
-    {
-      AssertDimension(M.n(), i1.size());
-      AssertDimension(M.m(), i2.size());
-      
-      if(mg_constrained_dofs == 0)
-      {
-      for (unsigned int j=0; j<i1.size(); ++j)
-       for (unsigned int k=0; k<i2.size(); ++k)
-         if (std::fabs(M(k,j)) >= threshold)
-           G.add(i1[j], i2[k], M(k,j));
-      }
-      else
-      {
-        for (unsigned int j=0; j<i1.size(); ++j)
-          for (unsigned int k=0; k<i2.size(); ++k)
-            if (std::fabs(M(j,k)) >= threshold)
-              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-              {
-                if (mg_constrained_dofs->set_boundary_values())
-                {
-                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) 
-                      ||
-                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                       mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                       i1[j] == i2[k]))
-                    G.add(i1[j], i2[k], M(k,j));
-                }
-                else
-                    G.add(i1[j], i2[k], M(k,j));
-              }
-      }
-    }
-    
-    
-    template <class MATRIX>
-    template <class DOFINFO>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info)
-    {
-      const unsigned int level = info.cell->level();
-      assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
-
-      if(mg_constrained_dofs != 0)
-      {
-        assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
-        assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
-      }
-    }
-    
-
-    template <class MATRIX>
-    template <class DOFINFO>
-    inline void
-    MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info1,
-                                    const DOFINFO& info2)
-    {
-      const unsigned int level1 = info1.cell->level();
-      const unsigned int level2 = info2.cell->level();
-      
-      if (level1 == level2)
-       {
-          if(mg_constrained_dofs == 0)
-          {
-            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
-          {
-            assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
-            assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
-            assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
-            assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
-          }
-       }
-      else
-       {
-         Assert(level1 > level2, ExcInternalError());
-                                          // Do not add info2.M1,
-                                          // which is done by
-                                          // the coarser cell
-          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);
-          if(level1>0)
-          {
-          assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
-          assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
-          }
-       }
-    }
-    
-    
 //----------------------------------------------------------------------//
     
     template <class MATRIX, typename number>
@@ -2052,61 +848,6 @@ namespace MeshWorker
        }
     }
 
-//----------------------------------------------------------------------//
-
-    template <class MATRIX, class VECTOR>
-    SystemSimple<MATRIX,VECTOR>::SystemSimple(double t)
-                   :
-                   MatrixSimple<MATRIX>(t)
-    {}
-
-
-    template <class MATRIX, class VECTOR>
-    inline void
-    SystemSimple<MATRIX,VECTOR>::initialize(MATRIX& m, VECTOR& rhs)
-    {
-      NamedData<VECTOR*> data;
-      VECTOR* p = &rhs;
-      data.add(p, "right hand side");
-      
-      MatrixSimple<MATRIX>::initialize(m);
-      ResidualSimple<VECTOR>::initialize(data);
-    }
-
-    
-    template <class MATRIX, class VECTOR>
-    template <class DOFINFO>
-    inline void
-    SystemSimple<MATRIX,VECTOR>::initialize_info(DOFINFO& info,
-                                                bool face) const
-    {
-      MatrixSimple<MATRIX>::initialize_info(info, face);
-      ResidualSimple<VECTOR>::initialize_info(info, face);
-    }
-
-
-    template <class MATRIX, class VECTOR>
-    template <class DOFINFO>
-    inline void
-    SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info)
-    {
-      MatrixSimple<MATRIX>::assemble(info);
-      ResidualSimple<VECTOR>::assemble(info);
-    }
-
-
-    template <class MATRIX, class VECTOR>
-    template <class DOFINFO>
-    inline void
-    SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info1,
-                                         const DOFINFO& info2)
-    {
-      MatrixSimple<MATRIX>::assemble(info1, info2);
-      ResidualSimple<VECTOR>::assemble(info1, info2);
-    }    
-  }
-}
-
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/deal.II/include/deal.II/meshworker/functional.h b/deal.II/include/deal.II/meshworker/functional.h
new file mode 100644 (file)
index 0000000..3a74174
--- /dev/null
@@ -0,0 +1,311 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__mesh_worker_functional_h
+#define __deal2__mesh_worker_functional_h
+
+#include <deal.II/base/named_data.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/mg_level_object.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+  namespace Assembler
+  {
+/**
+ * The class assembling local contributions to a functional into
+ * global functionals.
+ *
+ * 
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <typename number = double>
+    class Functional
+    {
+      public:
+                                        /**
+                                         * Initialize local data to
+                                         * store functionals. The
+                                         * number <tt>n</tt> is the
+                                         * number of functionals to
+                                         * be computed.
+                                         */
+       void initialize(unsigned int n);
+                                        /**
+                                         * Initialize the local data
+                                         * in the DoFInfo object used
+                                         * later for assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face);
+       
+                                        /**
+                                         * Assemble the local values
+                                         * into the global vectors.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local values
+                                         * into the global vectors.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+
+                                        /**
+                                         * The value of the ith entry
+                                         * in #results.
+                                         */
+       number operator() (unsigned int i) const;
+      private:
+                                        /**
+                                         * The values into which the
+                                         * results are added.
+                                         */
+       std::vector<double> results;
+    };
+
+/**
+ * Compute cell and face contributions of one or several functionals,
+ * typically for error estimates.
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <typename number = double>
+    class CellsAndFaces
+    {
+      public:
+                                        /**
+                                         * The data type for
+                                         * communicating the cell and
+                                         * face vectors.
+                                         */
+       typedef NamedData<BlockVector<number>*> DataVectors;
+       
+                                        /**
+                                         * The initialization
+                                         * function, specifying the
+                                         * #results vectors and
+                                         * whether face data should
+                                         * be collected separately.
+                                         *
+                                         * #results should contain
+                                         * two block vectors named
+                                         * "cells" and "faces" (the
+                                         * latter only if
+                                         * #separate_faces is
+                                         * true). In each of the two,
+                                         * each block should have
+                                         * equal size and be large
+                                         * enough to accomodate all
+                                         * user indices set in the
+                                         * cells and faces covered by
+                                         * the loop it is used
+                                         * in. Typically, for
+                                         * estimators, this is
+                                         * Triangulation::n_active_cells()
+                                         * and
+                                         * Triangulation::n_faces(),
+                                         * respectively.
+                                         *
+                                         * The use of BlockVector may
+                                         * seem cumbersome, but it
+                                         * allows us to assemble
+                                         * several functionals at the
+                                         * same time, one in each
+                                         * block. The typical
+                                         * situation for error
+                                         * estimate is just having a
+                                         * single block in each vector.
+                                         */
+       void initialize(DataVectors& results,
+                       bool separate_faces = true);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
+       
+                                        /**
+                                         * Assemble the local values
+                                         * into the global vectors.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local values
+                                         * into the global vectors.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+
+                                        /**
+                                         * The value of the ith entry
+                                         * in #results.
+                                         */
+       number operator() (unsigned int i) const;
+      private:
+       DataVectors results;
+       bool separate_faces;
+    };
+//----------------------------------------------------------------------//
+    
+    template <typename number>
+    inline void
+    Functional<number>::initialize(unsigned int n)
+    {
+      results.resize(n);
+      std::fill(results.begin(), results.end(), 0.);
+    }
+    
+    
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    Functional<number>::initialize_info(DOFINFO& info, bool)
+    {
+      info.initialize_numbers(results.size());
+    }
+    
+    
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    Functional<number>::assemble(const DOFINFO& info)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       results[i] += info.value(i);
+    }
+    
+
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    Functional<number>::assemble(const DOFINFO& info1,
+                                const DOFINFO& info2)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       {
+         results[i] += info1.value(i);
+         results[i] += info2.value(i);
+       }
+    }
+
+    
+    template <typename number>
+    inline number
+    Functional<number>::operator() (unsigned int i) const
+    {
+      AssertIndexRange(i, results.size());
+      return results[i];
+    }
+    
+//----------------------------------------------------------------------//
+
+    template <typename number>
+    inline void
+    CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
+    {
+      Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
+      if (sep)
+       {
+         Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
+         AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
+       }
+      
+      results = r;
+      separate_faces = sep;
+    }
+
+    
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    CellsAndFaces<number>::initialize_info(DOFINFO& info, bool) const
+    {
+      info.initialize_numbers(results(0)->n_blocks());
+    }
+
+
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    CellsAndFaces<number>::assemble(const DOFINFO& info)
+    {
+      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.value(i);
+         else
+           results(0)->block(i)(info.cell->user_index()) += info.value(i);
+       }
+    }
+    
+
+    template <typename number>
+    template <class DOFINFO>
+    inline void
+    CellsAndFaces<number>::assemble(const DOFINFO& info1,
+                                   const DOFINFO& info2)
+    {
+      for (unsigned int i=0;i<info1.n_values();++i)
+       {
+         if (separate_faces)
+           {
+             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.value(i);
+             results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
+           }
+       }
+    }
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h
new file mode 100644 (file)
index 0000000..81970cc
--- /dev/null
@@ -0,0 +1,1021 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__mesh_worker_assembler_h
+#define __deal2__mesh_worker_assembler_h
+
+#include <deal.II/base/named_data.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/mg_level_object.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/functional.h>
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+/**
+ * @file
+ * @brief MeshWorker::Assember objects for systems without block structure
+ *
+ * The header containing the classes
+ * MeshWorker::Assember::MatrixSimple,
+ * MeshWorker::Assember::MGMatrixSimple,
+ * MeshWorker::Assember::ResidualSimple, and
+ * MeshWorker::Assember::SystemSimple.
+ */
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+  namespace Assembler
+  {
+/**
+ * Assemble residuals without block structure.
+ *
+ * The data structure for this Assembler class is a simple vector on
+ * each cell with entries from zero to
+ * FiniteElementData::dofs_per_cell and a simple global vector with
+ * entries numbered from zero to DoFHandler::n_dofs(). No BlockInfo is
+ * required and the global vector may be any type of vector having
+ * element access through <tt>operator() (unsigned int)</tt>
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <class VECTOR>
+    class ResidualSimple
+    {
+      public:
+       void initialize(NamedData<VECTOR*>& results);
+                                        /**
+                                         * Initialize the constraints. 
+                                         */
+        void initialize(const ConstraintMatrix& constraints);
+                                        /**
+                                         * Initialize the local data
+                                         * in the DoFInfo object used
+                                         * later for assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
+       
+                                        /**
+                                         * Assemble the local residuals
+                                         * into the global residuals.
+                                         */
+       template <class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local residuals
+                                         * into the global residuals.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+      private:
+                                        /**
+                                         * The global residal vectors
+                                         * filled by assemble().
+                                         */
+       NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
+      /**
+       * A pointer to the object containing constraints.
+       */
+      SmartPointer<const ConstraintMatrix,ResidualSimple<VECTOR> > constraints;
+    };
+    
+/**
+ * Assemble local matrices into a single global matrix without using
+ * block structure.
+ *
+ * After being initialized with a SparseMatrix object (or another
+ * matrix offering the same functionality as SparseMatrix::add()),
+ * this class can be used in a MeshWorker::loop() to assemble the cell
+ * and face matrices into the global matrix.
+ *
+ * @todo On locally refined meshes, a ConstraintMatrix should be used
+ * to automatically eliminate hanging nodes.
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <class MATRIX>
+    class MatrixSimple
+    {
+      public:
+                                        /**
+                                         * Constructor, initializing
+                                         * the #threshold, which
+                                         * limits how small numbers
+                                         * may be to be entered into
+                                         * the matrix.
+                                         */
+       MatrixSimple(double threshold = 1.e-12);
+       
+                                        /**
+                                         * Store the result matrix
+                                         * for later assembling.
+                                         */
+       void initialize(MATRIX& m);
+                                        /**
+                                         * Initialize the constraints. 
+                                         */
+        void initialize(const ConstraintMatrix& constraints);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
+       
+                                        /**
+                                         * Assemble the matrix
+                                         * DoFInfo::M1[0]
+                                         * into the global matrix.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local
+                                         * matrices in the info
+                                         * objects into the global
+                                         * matrix.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+      private:
+                                        /**
+                                         * Assemble a single matrix
+                                         * into #matrix.
+                                         */
+       void assemble(const FullMatrix<double>& M,
+                     const std::vector<unsigned int>& i1,
+                     const std::vector<unsigned int>& i2);
+       
+                                        /**
+                                         * The global matrix being
+                                         * assembled.
+                                         */
+       SmartPointer<MATRIX,MatrixSimple<MATRIX> > matrix;
+      /**
+       * A pointer to the object containing constraints.
+       */
+        SmartPointer<const ConstraintMatrix,MatrixSimple<MATRIX> > constraints;
+       
+                                        /**
+                                         * The smallest positive
+                                         * number that will be
+                                         * entered into the global
+                                         * matrix. All smaller
+                                         * absolute values will be
+                                         * treated as zero and will
+                                         * not be assembled.
+                                         */
+       const double threshold;
+       
+    };
+    
+    
+/**
+ * Assemble local matrices into level matrices without using
+ * block structure.
+ *
+ * @todo The matrix structures needed for assembling level matrices
+ * with local refinement and continuous elements are missing.
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <class MATRIX>
+    class MGMatrixSimple
+    {
+      public:
+                                        /**
+                                         * Constructor, initializing
+                                         * the #threshold, which
+                                         * limits how small numbers
+                                         * may be to be entered into
+                                         * the matrix.
+                                         */
+       MGMatrixSimple(double threshold = 1.e-12);
+       
+                                        /**
+                                         * Store the result matrix
+                                         * for later assembling.
+                                         */
+       void initialize(MGLevelObject<MATRIX>& m);
+
+                                        /**
+                                         * Initialize the multilevel
+                                          * constraints. 
+                                         */
+        void initialize(const MGConstrainedDoFs& mg_constrained_dofs);
+
+                                        /**
+                                         * Initialize the matrices
+                                         * #flux_up and #flux_down
+                                         * used for local refinement
+                                         * with discontinuous
+                                         * Galerkin methods.
+                                         */
+       void initialize_fluxes(MGLevelObject<MATRIX>& flux_up,
+                              MGLevelObject<MATRIX>& flux_down);
+       
+                                        /**
+                                         * Initialize the matrices
+                                         * #interface_in and #interface_out
+                                         * used for local refinement
+                                         * with continuous
+                                         * Galerkin methods.
+                                         */
+
+       void initialize_interfaces(MGLevelObject<MATRIX>& interface_in,
+                              MGLevelObject<MATRIX>& interface_out);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
+       
+                                        /**
+                                         * Assemble the matrix
+                                         * DoFInfo::M1[0]
+                                         * into the global matrix.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local
+                                         * matrices in the info
+                                         * objects into the global
+                                         * matrices.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+      private:
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+    void assemble(MATRIX& G,
+                 const FullMatrix<double>& M,
+                 const std::vector<unsigned int>& i1,
+                 const std::vector<unsigned int>& i2);
+
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+    void assemble(MATRIX& G,
+                 const FullMatrix<double>& M,
+                 const std::vector<unsigned int>& i1,
+                 const std::vector<unsigned int>& i2,
+                  const unsigned int level);
+       
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+
+    void assemble_up(MATRIX& G,
+                           const FullMatrix<double>& M,
+                           const std::vector<unsigned int>& i1,
+                           const std::vector<unsigned int>& i2,
+                            const unsigned int level = -1);
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+
+    void assemble_down(MATRIX& G,
+                           const FullMatrix<double>& M,
+                           const std::vector<unsigned int>& i1,
+                           const std::vector<unsigned int>& i2,
+                            const unsigned int level = -1);
+
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+
+    void assemble_in(MATRIX& G,
+                 const FullMatrix<double>& M,
+                 const std::vector<unsigned int>& i1,
+                 const std::vector<unsigned int>& i2,
+                  const unsigned int level = -1);
+
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+
+    void assemble_out(MATRIX& G,
+                 const FullMatrix<double>& M,
+                 const std::vector<unsigned int>& i1,
+                 const std::vector<unsigned int>& i2,
+                  const unsigned int level = -1);
+       
+                                        /**
+                                         * The global matrix being
+                                         * assembled.
+                                         */
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > matrix;
+       
+                                        /**
+                                         * The matrix used for face
+                                         * flux terms across the
+                                         * refinement edge, coupling
+                                         * coarse to fine.
+                                         */
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_up;
+       
+                                        /**
+                                         * The matrix used for face
+                                         * flux terms across the
+                                         * refinement edge, coupling
+                                         * fine to coarse.
+                                         */
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_down;
+
+                                        /**
+                                         * The matrix used for face
+                                         * contributions for continuous
+                                          * elements across the
+                                         * refinement edge, coupling
+                                         * coarse to fine.
+                                         */
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_in;
+       
+                                        /**
+                                         * The matrix used for face
+                                         * contributions for continuous
+                                          * elements across the
+                                         * refinement edge, coupling
+                                         * fine to coarse.
+                                         */
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_out;
+      /**
+       * A pointer to the object containing constraints.
+       */
+        SmartPointer<const MGConstrainedDoFs,MGMatrixSimple<MATRIX> > mg_constrained_dofs;
+       
+                                        /**
+                                         * The smallest positive
+                                         * number that will be
+                                         * entered into the global
+                                         * matrix. All smaller
+                                         * absolute values will be
+                                         * treated as zero and will
+                                         * not be assembled.
+                                         */
+       const double threshold;
+       
+    };
+    
+    
+/**
+ * Assemble a simple matrix and a simple right hand side at once. We
+ * use a combination of MatrixSimple and ResidualSimple to achieve
+ * this. Cell and face operators should fill the matrix and vector
+ * objects in LocalResults and this class will assemble
+ * them into matrix and vector objects.
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+    template <class MATRIX, class VECTOR>
+    class SystemSimple :
+       private MatrixSimple<MATRIX>,
+       private ResidualSimple<VECTOR>
+    {
+      public:
+                                        /**
+                                         * Constructor setting the
+                                         * threshold value in
+                                         * MatrixSimple.
+                                         */
+       SystemSimple(double threshold = 1.e-12);
+
+                                        /**
+                                         * Store the two objects data
+                                         * is assembled into.
+                                         */
+       void initialize(MATRIX& m, VECTOR& rhs);
+       
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The info object refers to
+                                         * a cell if
+                                         * <code>!face</code>, or
+                                         * else to an interior or
+                                         * boundary face.
+                                         */
+       template <class DOFINFO>
+       void initialize_info(DOFINFO& info, bool face) const;
+       
+                                        /**
+                                         * Assemble the matrix
+                                         * DoFInfo::M1[0]
+                                         * into the global matrix.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info);
+       
+                                        /**
+                                         * Assemble both local
+                                         * matrices in the info
+                                         * objects into the global
+                                         * matrix.
+                                         */
+       template<class DOFINFO>
+       void assemble(const DOFINFO& info1,
+                     const DOFINFO& info2);
+    };
+
+    
+//----------------------------------------------------------------------//
+
+    template <class VECTOR>
+    inline void
+    ResidualSimple<VECTOR>::initialize(NamedData<VECTOR*>& results)
+    {
+      residuals = results;
+    }
+
+    template <class VECTOR>
+    inline void
+    ResidualSimple<VECTOR>::initialize(const ConstraintMatrix& c)
+    {
+      constraints = &c;
+    }
+
+    
+    template <class VECTOR>
+    template <class DOFINFO>
+    inline void
+    ResidualSimple<VECTOR>::initialize_info(DOFINFO& info, bool) const
+    {
+      info.initialize_vectors(residuals.size());
+    }
+
+    
+    template <class VECTOR>
+    template <class DOFINFO>
+    inline void
+    ResidualSimple<VECTOR>::assemble(const DOFINFO& info)
+    {
+      for (unsigned int k=0;k<residuals.size();++k)
+      {
+        if(constraints == 0)
+        {
+          for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
+            (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
+        }
+        else
+          constraints->distribute_local_to_global(
+              info.vector(k).block(0), info.indices, (*residuals(k)));
+      }
+    }
+
+    
+    template <class VECTOR>
+    template <class DOFINFO>
+    inline void
+    ResidualSimple<VECTOR>::assemble(const DOFINFO& info1,
+                                       const DOFINFO& info2)
+    {
+      for (unsigned int k=0;k<residuals.size();++k)
+       {
+          if(constraints == 0)
+          {
+            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);
+          }
+          else
+          {
+            constraints->distribute_local_to_global(
+              info1.vector(k).block(0), info1.indices, (*residuals(k)));
+            constraints->distribute_local_to_global(
+              info2.vector(k).block(0), info2.indices, (*residuals(k)));
+          }
+       }
+    }
+
+    
+//----------------------------------------------------------------------//
+
+    template <class MATRIX>
+    inline
+    MatrixSimple<MATRIX>::MatrixSimple(double threshold)
+                   :
+                   threshold(threshold)
+    {}
+    
+    
+    template <class MATRIX>
+    inline void
+    MatrixSimple<MATRIX>::initialize(MATRIX& m)
+    {
+      matrix = &m;
+    }
+
+
+    template <class MATRIX>
+    inline void
+    MatrixSimple<MATRIX>::initialize(const ConstraintMatrix& c)
+    {
+      constraints = &c;
+    }
+    
+
+    template <class MATRIX >
+    template <class DOFINFO>
+    inline void
+    MatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
+    {
+      info.initialize_matrices(1, face);
+    }
+
+
+
+    template <class MATRIX>
+    inline void
+    MatrixSimple<MATRIX>::assemble(const FullMatrix<double>& M,
+                                  const std::vector<unsigned int>& i1,
+                                  const std::vector<unsigned int>& i2)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+     
+     if(constraints == 0)
+     { 
+      for (unsigned int j=0; j<i1.size(); ++j)
+       for (unsigned int k=0; k<i2.size(); ++k)
+         if (std::fabs(M(j,k)) >= threshold)
+           matrix->add(i1[j], i2[k], M(j,k));
+     }
+     else
+       constraints->distribute_local_to_global(
+           M, i1, i2, *matrix);
+    }
+    
+    
+    template <class MATRIX>
+    template <class DOFINFO>
+    inline void
+    MatrixSimple<MATRIX>::assemble(const DOFINFO& info)
+    {
+      assemble(info.matrix(0,false).matrix, info.indices, info.indices);
+    }
+    
+
+    template <class MATRIX>
+    template <class DOFINFO>
+    inline void
+    MatrixSimple<MATRIX>::assemble(const DOFINFO& info1,
+                                  const DOFINFO& info2)
+    {
+      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);
+    }
+    
+    
+//----------------------------------------------------------------------//
+
+    template <class MATRIX>
+    inline
+    MGMatrixSimple<MATRIX>::MGMatrixSimple(double threshold)
+                   :
+                   threshold(threshold)
+    {}
+    
+    
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize(MGLevelObject<MATRIX>& m)
+    {
+      matrix = &m;
+    }
+    
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize(const MGConstrainedDoFs& c)
+    {
+      mg_constrained_dofs = &c;
+    }
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize_fluxes(
+      MGLevelObject<MATRIX>& up, MGLevelObject<MATRIX>& down)
+    {
+      flux_up = &up;
+      flux_down = &down;
+    }
+    
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize_interfaces(
+      MGLevelObject<MATRIX>& in, MGLevelObject<MATRIX>& out)
+    {
+      interface_in = &in;
+      interface_out = &out;
+    }
+
+
+    template <class MATRIX >
+    template <class DOFINFO>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
+    {
+      info.initialize_matrices(1, face);
+    }
+
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              G.add(i1[j], i2[k], M(j,k));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+            {
+                if(!mg_constrained_dofs->continuity_across_refinement_edges())
+                  G.add(i1[j], i2[k], M(j,k));
+            }
+      }
+    }
+
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              G.add(i1[j], i2[k], M(j,k));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+              {
+                if (mg_constrained_dofs->set_boundary_values())
+                {
+                  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                        !mg_constrained_dofs->is_boundary_index(level, i2[k]))
+                      ||
+                      (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                       mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
+                       i1[j] == i2[k]))
+                    G.add(i1[j], i2[k], M(j,k));
+                }
+                else
+                    G.add(i1[j], i2[k], M(j,k));
+            }
+      }
+    }
+    
+    
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_up(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.n(), i1.size());
+      AssertDimension(M.m(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(k,j)) >= threshold)
+              G.add(i1[j], i2[k], M(k,j));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(k,j)) >= threshold)
+              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+                G.add(i1[j], i2[k], M(k,j));
+      }
+    }
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_down(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              G.add(i1[j], i2[k], M(j,k));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+                  G.add(i1[j], i2[k], M(j,k));
+      }
+    }
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_in(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              G.add(i1[j], i2[k], M(j,k));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+              {
+                if (mg_constrained_dofs->set_boundary_values())
+                {
+                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                      ||
+                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                       mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                       i1[j] == i2[k]))
+                    G.add(i1[j], i2[k], M(j,k));
+                }
+                else
+                  G.add(i1[j], i2[k], M(j,k));
+              }
+      }
+    }
+    
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_out(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.n(), i1.size());
+      AssertDimension(M.m(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+      for (unsigned int j=0; j<i1.size(); ++j)
+       for (unsigned int k=0; k<i2.size(); ++k)
+         if (std::fabs(M(k,j)) >= threshold)
+           G.add(i1[j], i2[k], M(k,j));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(j,k)) >= threshold)
+              if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+                  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+              {
+                if (mg_constrained_dofs->set_boundary_values())
+                {
+                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) 
+                      ||
+                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                       mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                       i1[j] == i2[k]))
+                    G.add(i1[j], i2[k], M(k,j));
+                }
+                else
+                    G.add(i1[j], i2[k], M(k,j));
+              }
+      }
+    }
+    
+    
+    template <class MATRIX>
+    template <class DOFINFO>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info)
+    {
+      const unsigned int level = info.cell->level();
+      assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
+
+      if(mg_constrained_dofs != 0)
+      {
+        assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
+        assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
+      }
+    }
+    
+
+    template <class MATRIX>
+    template <class DOFINFO>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info1,
+                                    const DOFINFO& info2)
+    {
+      const unsigned int level1 = info1.cell->level();
+      const unsigned int level2 = info2.cell->level();
+      
+      if (level1 == level2)
+       {
+          if(mg_constrained_dofs == 0)
+          {
+            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
+          {
+            assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
+            assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
+            assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
+            assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          }
+       }
+      else
+       {
+         Assert(level1 > level2, ExcInternalError());
+                                          // Do not add info2.M1,
+                                          // which is done by
+                                          // the coarser cell
+          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);
+          if(level1>0)
+          {
+          assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          }
+       }
+    }
+//----------------------------------------------------------------------//
+
+    template <class MATRIX, class VECTOR>
+    SystemSimple<MATRIX,VECTOR>::SystemSimple(double t)
+                   :
+                   MatrixSimple<MATRIX>(t)
+    {}
+
+
+    template <class MATRIX, class VECTOR>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::initialize(MATRIX& m, VECTOR& rhs)
+    {
+      NamedData<VECTOR*> data;
+      VECTOR* p = &rhs;
+      data.add(p, "right hand side");
+      
+      MatrixSimple<MATRIX>::initialize(m);
+      ResidualSimple<VECTOR>::initialize(data);
+    }
+
+    
+    template <class MATRIX, class VECTOR>
+    template <class DOFINFO>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::initialize_info(DOFINFO& info,
+                                                bool face) const
+    {
+      MatrixSimple<MATRIX>::initialize_info(info, face);
+      ResidualSimple<VECTOR>::initialize_info(info, face);
+    }
+
+
+    template <class MATRIX, class VECTOR>
+    template <class DOFINFO>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info)
+    {
+      MatrixSimple<MATRIX>::assemble(info);
+      ResidualSimple<VECTOR>::assemble(info);
+    }
+
+
+    template <class MATRIX, class VECTOR>
+    template <class DOFINFO>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info1,
+                                         const DOFINFO& info2)
+    {
+      MatrixSimple<MATRIX>::assemble(info1, info2);
+      ResidualSimple<VECTOR>::assemble(info1, info2);
+    }    
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.