]> https://gitweb.dealii.org/ - dealii.git/commitdiff
first example for the generic MeshWorker::loop() is up and running
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 2 Oct 2009 00:08:27 +0000 (00:08 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 2 Oct 2009 00:08:27 +0000 (00:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@19654 0785d39b-7218-0410-832d-ea1e28bc413d

16 files changed:
deal.II/deal.II/include/numerics/mesh_worker.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_assembler.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_info.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_info.templates.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_loop.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.inst.in [new file with mode: 0644]
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/fe_field_function.inst.in
deal.II/deal.II/source/numerics/mesh_worker.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker_info.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker_info.inst.in [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in [new file with mode: 0644]
deal.II/examples/step-38/step-38.cc

diff --git a/deal.II/deal.II/include/numerics/mesh_worker.h b/deal.II/deal.II/include/numerics/mesh_worker.h
new file mode 100644 (file)
index 0000000..df1e067
--- /dev/null
@@ -0,0 +1,411 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
+//
+//    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_h
+#define __deal2__mesh_worker_h
+
+#include <base/geometry_info.h>
+#include <base/named_data.h>
+#include <lac/block_indices.h>
+#include <multigrid/mg_level_object.h>
+#include <numerics/mesh_worker_vector_selector.h>
+#include <numerics/mesh_worker_info.h>
+#include <numerics/mesh_worker_assembler.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+template<int,int> class DoFHandler;
+template<int,int> class MGDoFHandler;
+
+/**
+ * A collection of functions and classes for the mesh loops that are
+ * an ubiquitous part of each finite element program.
+ *
+ * The workhorse of this namespace is the loop() function, which implements a
+ * completely generic loop over all mesh cells.
+ *
+ * The loop() depends on certain objects handed to it as
+ * arguments. These objects are of two types, info objects like
+ * DoFInfo and IntegrationInfo and worker objects like LocalWorker and
+ * IntegrationWorker.
+ *
+ * Worker objects usually do two different jobs: first, they compute
+ * the local contribution of a cell or face to the global
+ * operation. Second, they assemble this local contribution into the
+ * global result, whether a functional, a form or a bilinear
+ * form. While the first job is particular to the problem being
+ * solved, the second is generic and only depends on the data
+ * structures. Therefore, base classes for workers assimbling into
+ * global data are provided in the namespace Assembler.
+ *
+ * <h3>Simplified interfaces</h3>
+ *
+ * Since the loop() is fairly general, a specialization
+ * integration_loop() is available, which is a wrapper around loop()
+ * with a simplified interface.
+ *
+ * This integration_loop() is complemented by a class
+ * AssemblingIntegrator, which separates local integration from
+ * assembling. In order to use it, follow this recipe: fist, create a
+ * class responsible for the local integration:
+ *
+ * @begin{code}
+ * @end{code}
+ *
+ * @author Guido Kanschat, 2009
+ */
+namespace MeshWorker
+{
+/**
+ * Template for a class for the objects doing the actual work on cells
+ * and faces.
+ *
+ * This class can serve as a base class for the actual worker class,
+ * since, while we do not use virtual functions, we provide the
+ * necessary interface for the mesh loops and the DoFInfo
+ * class here. Thus, they do not have to be reprogrammed.
+ *
+ * In particular, the mesh loops will require data elements
+ * #interior_fluxes and #boundary_fluxes to determine whether the
+ * loop over faces will be started at all.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template <int dim>
+  class LocalWorker
+  {
+    public:
+                                      /**
+                                       * Constructor, setting
+                                       * #interior_fluxes and
+                                       * #boundary_fluxes to @p
+                                       * true.
+                                       */
+      LocalWorker ();
+       
+                                      /**
+                                       * Do the work on a cell.
+                                       */
+      void cell(DoFInfo<dim>& cell);
+      
+                                      /**
+                                       * Do the work on a boundary face.
+                                       */
+      void bdry(DoFInfo<dim>& face);
+      
+                                      /**
+                                       * Do the work on an interior face.
+                                       */
+      void face(DoFInfo<dim>& face1, DoFInfo<dim>& face2);
+
+                                      /**
+                                       * Computations on interior
+                                       * faces are necessary.
+                                       */
+      bool interior_fluxes;
+      
+                                      /**
+                                       * Computations on interior
+                                       * faces are necessary.
+                                       */
+      bool boundary_fluxes;      
+  };
+
+/**
+ * Worker object for integration of functionals, residuals or matrices.
+ *
+ */
+  template <int dim>
+  class IntegrationWorker : public LocalWorker<dim>
+  {
+    public:
+                                      /**
+                                       * The info type expected by a
+                                       * cell integrator.
+                                       */
+      typedef IntegrationInfo<dim, FEValuesBase<dim, dim>, dim> CellInfo;
+      
+                                      /**
+                                       * The info type expected by a
+                                       * face integrator.
+                                       */
+      typedef IntegrationInfo<dim, FEFaceValuesBase<dim, dim>, dim> FaceInfo;
+      
+                                      /**
+                                       * Initialize default values.
+                                       */
+      IntegrationWorker();
+                                      /**
+                                       * Initialize the
+                                       * VectorSelector objects
+                                       * #cell_selector,
+                                       * #bdry_selector and
+                                       * #face_selector in order to
+                                       * save computational
+                                       * eeffort. If no selectors
+                                       * are used, then values for
+                                       * all named vectors in
+                                       * #global_data will be
+                                       * computed in all quadrature
+                                       * points.
+                                       *
+                                       * This function will also
+                                       * add UpdateFlags to the
+                                       * flags stored in this class.
+                                       */
+      void initialize_selectors(const VectorSelector& cell_selector,
+                               const VectorSelector& bdry_selector,
+                               const VectorSelector& face_selector);
+
+                                      /**
+                                       * Add a vector to some or
+                                       * all selectors.
+                                       */
+      void add_selector(const std::string& name, bool values, bool gradients, bool hessians,
+                       bool cell, bool bdry, bool face);
+
+                                      /**
+                                       * Add additional values for update.
+                                       */
+      void add_update_flags(const UpdateFlags flags, bool cell = true,
+                           bool bdry = true, bool face = true,
+                           bool ngbr = true);
+       
+                                      /** Assign n-point Gauss
+                                       * quadratures to each of the
+                                       * quadrature rules. Here, a
+                                       * size of zero points means
+                                       * that no loop over these grid
+                                       * entities should be
+                                       * performed.
+                                       */
+      void initialize_gauss_quadrature(unsigned int n_cell_points,
+                                      unsigned int n_bdry_points,
+                                      unsigned int n_face_points);
+
+                                      /**
+                                       * Select the vectors from
+                                       * DoFInfo::global_data
+                                       * that should be computed in
+                                       * the quadrature points on cells.
+                                       */
+      MeshWorker::VectorSelector cell_selector;
+       
+                                      /**
+                                       * Select the vectors from
+                                       * DoFInfo::global_data
+                                       * that should be computed in
+                                       * the quadrature points on
+                                       * boundary faces.
+                                       */
+      MeshWorker::VectorSelector bdry_selector;
+       
+                                      /**
+                                       * Select the vectors from
+                                       * DoFInfo::global_data
+                                       * that should be computed in
+                                       * the quadrature points on
+                                       * interior faces.
+                                       */
+      MeshWorker::VectorSelector face_selector;
+
+                                      /**
+                                       * The set of update flags
+                                       * for boundary cell integration.
+                                       *
+                                       * Defaults to
+                                       * #update_JxW_values.
+                                       */
+      UpdateFlags cell_flags;
+                                      /**
+                                       * The set of update flags
+                                       * for boundary face integration.
+                                       *
+                                       * Defaults to
+                                       * #update_JxW_values and
+                                       * #update_normal_vectors.
+                                       */
+      UpdateFlags bdry_flags;
+       
+                                      /**
+                                       * The set of update flags
+                                       * for interior face integration.
+                                       *
+                                       * Defaults to
+                                       * #update_JxW_values and
+                                       * #update_normal_vectors.
+                                       */
+      UpdateFlags face_flags;
+
+                                      /**
+                                       * The set of update flags
+                                       * for interior face integration.
+                                       *
+                                       * Defaults to
+                                       * #update_default, since
+                                       * quadrature weights are
+                                       * taken from the other cell.
+                                       */
+      UpdateFlags ngbr_flags;
+       
+                                      /**
+                                       * The quadrature rule used
+                                       * on cells.
+                                       */
+      Quadrature<dim> cell_quadrature;
+      
+                                      /**
+                                       * The quadrature rule used
+                                       * on boundary faces.
+                                       */
+      Quadrature<dim-1> bdry_quadrature;
+      
+                                      /**
+                                       * The quadrature rule used
+                                       * on interior faces.
+                                       */
+      Quadrature<dim-1> face_quadrature;
+  };
+      
+  
+/**
+ * A worker class for integration_loop(), separating assembling and
+ * local integration. Objects of this class rely on the base class
+ * provided by the template argument ASSEMBLER for assembling local
+ * data into global. The integration of local data is delegated to the
+ * other template parameter class INTEGRATOR.
+ *
+ * In order to use this class, it will be necessary to create an
+ * INTEGRATOR class following this template:
+ *
+ * @begin{code}
+ * template <int dim>
+ * class MyLocalOperator : public Subscriptor
+ * {
+ *   public:
+ *     void cell(typename IntegrationWorker<dim>::CellInfo& info) const;
+ *     void bdry(typename IntegrationWorker<dim>::FaceInfo& info) const;
+ *     void face(typename IntegrationWorker<dim>::FaceInfo& info1,
+ *               typename IntegrationWorker<dim>::FaceInfo& info2) const;
+ * };
+ * @end{code}
+ *
+ * This class will do whatever your problem requires locally on each
+ * cell and/or face. Once this class is defined, you choose a suitable
+ * assembler for your problem from the Assembler namespace and set up
+ * objects:
+ *
+ * @begin{code}
+ * MyLocalOperator<dim> myop;
+ *
+ * AssemblingIntegrator<dim, Assembler::MyAssembler, MyLocalOperator<dim> >
+ *   integrator(myop);
+ * @end{code}
+ *
+ * You do the necessary initializations of this @p integrator and then
+ * you have a worker object suitable for integration_loop().
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template <int dim, class ASSEMBLER, class INTEGRATOR>
+  class AssemblingIntegrator :
+      public IntegrationWorker<dim>,
+      public ASSEMBLER
+  {
+    public:
+                                      /**
+                                       * Constructor, initializing
+                                       * the local data.
+                                       */
+      AssemblingIntegrator(const INTEGRATOR& local);
+                                      /**
+                                       * The cell operator called by
+                                       * integration_loop().
+                                       */
+      void cell(typename IntegrationWorker<dim>::CellInfo& info);
+      
+                                      /**
+                                       * The local boundary operator
+                                       * called by
+                                       * integration_loop().
+                                       */
+      void bdry(typename IntegrationWorker<dim>::FaceInfo& info);
+      
+                                      /**
+                                       * The interior face operator
+                                       * called by
+                                       * integration_loop().
+                                       */
+      void face(typename IntegrationWorker<dim>::FaceInfo& info1,
+               typename IntegrationWorker<dim>::FaceInfo& info2);
+      
+    private:
+                                      /**
+                                       * Pointer to the object doing
+                                       * local integration.
+                                       */
+      SmartPointer<const INTEGRATOR> local;
+  };
+  
+//----------------------------------------------------------------------//
+  template <int dim>
+  inline
+  LocalWorker<dim>::LocalWorker()
+                 :
+                 interior_fluxes(true), boundary_fluxes(true)
+  {}
+
+//----------------------------------------------------------------------//
+
+  template <int dim, class ASSEMBLER, class INTEGRATOR>
+  inline
+  AssemblingIntegrator<dim,ASSEMBLER,INTEGRATOR>::AssemblingIntegrator(const INTEGRATOR& local)
+                 :
+                 local(&local, typeid(*this).name())
+  {}
+
+  
+  template <int dim, class ASSEMBLER, class INTEGRATOR>
+  inline void
+  AssemblingIntegrator<dim,ASSEMBLER,INTEGRATOR>::cell(
+    typename IntegrationWorker<dim>::CellInfo& info)
+  {
+    local->cell(info);
+    ASSEMBLER::assemble(info);
+  }
+
+  
+  template <int dim, class ASSEMBLER, class INTEGRATOR>
+  inline void
+  AssemblingIntegrator<dim,ASSEMBLER,INTEGRATOR>::bdry(
+    typename IntegrationWorker<dim>::FaceInfo& info)
+  {
+    local->bdry(info);
+    ASSEMBLER::assemble(info);
+  }
+
+  
+  template <int dim, class ASSEMBLER, class INTEGRATOR>
+  inline void
+  AssemblingIntegrator<dim,ASSEMBLER,INTEGRATOR>::face(
+    typename IntegrationWorker<dim>::FaceInfo& info1,
+    typename IntegrationWorker<dim>::FaceInfo& info2)
+  {
+    local->face(info1, info2);
+    ASSEMBLER::assemble(info1, info2);
+  }  
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h
new file mode 100644 (file)
index 0000000..b71701d
--- /dev/null
@@ -0,0 +1,1308 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
+//
+//    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 <numerics/mesh_worker_info.h>
+#include <base/named_data.h>
+#include <base/smartpointer.h>
+#include <lac/block_vector.h>
+#include <multigrid/mg_level_object.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+/**
+ * The namespace containing objects that can be used to assemble data
+ * computed on cells and faces into global objects. This can reach
+ * from collecting the total error estimate from cell and face
+ * contributions to assembling matrices and multilevel matrices.
+ *
+ * <h3>Data models</h3>
+ *
+ * The class chosen from this namespace determines which data model is
+ * used. For the local as well as the global objects, we have the
+ * choice between two models:
+ *
+ * <h4>The comprehensive data model</h4>
+ *
+ * This is the structure set up by the FESystem class. Globally, this
+ * means, data is assembled into one residual vector and into one
+ * matrix. These objects may be block vectors and block matrices, but
+ * the process of assembling ignores this fact.
+ *
+ * Similarly, there is only a single cell vector and cell matrix,
+ * respectively, which is indexed by all degrees of freedom of the
+ * FESystem. When building the cell matrix, it is necessary to
+ * distinguish between the different components of the system and
+ * select the right operator for each pair.
+ *
+ * <h4>The blocked data model</h4>
+ *
+ * Here, all the blocks are treated separately (in spite of using
+ * FESystem for its convenience in other places). For instance, no
+ * block matrix is assembled, but a list of blocks, which can be
+ * combined later by BlockMatrixArray. Locally, this means, that each
+ * matrix block of a system is generated separately and assembled into
+ * the corresponding global block.
+ *
+ * This approach is advantageous, if the number of matrices for each
+ * block position in the global system is different. For instance,
+ * block preconditioners for the Oseen problem require 3 pressure
+ * matrices, but only one divergence and one advection-diffusion
+ * operator for velocities.
+ *
+ * Additionally, this approach enables the construction of a system of
+ * equations from building blocks for each equation and coupling
+ * operator.
+ *
+ * Nevertheless, since a separate FEValues object must be created for
+ * each base element, it is not quite clear a priori, which data model
+ * is more efficient.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  namespace Assembler
+  {
+/**
+ * The class assembling local contributions to a functional into the
+ * global functionals.
+ *
+ * 
+ *
+ * @author Guido Kanschat, 2009
+ */
+    template <typename number = double>
+    class Functional
+    {
+      public:
+                                        /**
+                                         * Initialize local data to
+                                         * store functionals.
+                                         */
+       void initialize(unsigned int n);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face);
+       
+                                        /**
+                                         * Assemble the local values
+                                         * into the global vectors.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local values
+                                         * into the global vectors.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& 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.
+ *
+ * @author Guido Kanschat, 2009
+ */
+    template <typename number = double>
+    class CellsAndFaces
+    {
+      public:
+                                        /**
+                                         * The data type for
+                                         * communicating the cell and
+                                         * face vectors.
+                                         */
+       typedef NamedData<SmartPointer<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 second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face);
+       
+                                        /**
+                                         * Assemble the local values
+                                         * into the global vectors.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local values
+                                         * into the global vectors.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& 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>
+ */
+    template <class VECTOR>
+    class ResidualSimple
+    {
+      public:
+       void initialize(NamedData<SmartPointer<VECTOR> >& results);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face) const;
+       
+                                        /**
+                                         * Assemble the local residuals
+                                         * into the global residuals.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local residuals
+                                         * into the global residuals.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& info2);
+      private:
+                                        /**
+                                         * The global residal vectors
+                                         * filled by assemble().
+                                         */
+       NamedData<SmartPointer<VECTOR> > residuals;
+    };
+    
+/**
+ * Assemble local residuals into global residuals.
+ *
+ * The global residuals are expected as an FEVectors object.
+ * The local residuals are block vectors.
+ *
+ * Depending on whether the BlockInfo object was initialize with
+ * BlockInfo::initialize_local(), the comprehensive or block data
+ * model is used locally.
+ *
+ * In the block model, each of the blocks of the local vectors
+ * corresponds to the restriction of a single block of the system to
+ * this cell (@ref GlossBlock). Thus, the size of this local block is
+ * the number of degrees of freedom of the corresponding base element
+ * of the FESystem.
+ *
+ * @todo Comprehensive model currently not implemented.
+ *
+ * @author Guido Kanschat, 2009
+ */
+    template <class VECTOR>
+    class ResidualLocalBlocksToGlobalBlocks
+    {
+      public:
+                                        /**
+                                         * Copy the BlockInfo and the
+                                         * matrix pointers into local
+                                         * variables.
+                                         */
+       void initialize(NamedData<SmartPointer<VECTOR> >& residuals);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face) const;
+       
+
+                                        /**
+                                         * Assemble the local residuals
+                                         * into the global residuals.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local residuals
+                                         * into the global residuals.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& info2);
+      private:
+                                        /**
+                                         * Assemble a single local
+                                         * residual into the global.
+                                         */
+       void assemble(VECTOR global,
+                     const BlockVector<double>& local,
+                     const BlockIndices& bi,
+                     const std::vector<unsigned int>& dof);
+       
+                                        /**
+                                         * The global matrices,
+                                         * stored as a vector of
+                                         * pointers.
+                                         */
+       NamedData<SmartPointer<VECTOR> > residuals;
+    };
+
+
+/**
+ * Assemble local matrices into a single global matrix without using
+ * block structure.
+ *
+ * @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 local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face) const;
+       
+                                        /**
+                                         * Assemble the matrix
+                                         * DoFInfo::M1[0]
+                                         * into the global matrix.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local
+                                         * matrices in the info
+                                         * objects into the global
+                                         * matrix.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& 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> matrix;
+       
+                                        /**
+                                         * 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;
+       
+    };
+    
+    
+/**
+ * A helper class assembling local matrices into global matrices.
+ *
+ * The global matrices are expected as a vector of MatrixBlock
+ * objects, each containing a matrix object with a function
+ * corresponding to SparseMatrix::add() and information on the block
+ * row and column this matrix represents in a block system.
+ *
+ * The local matrices are expected as a similar vector of MatrixBlock
+ * objects, but containing a FullMatrix.
+ *
+ * Like with ResidualLocalBlocksToGlobalBlocks, the initialization of
+ * the BlockInfo object decides whether the comprehensive data model
+ * or the block model is used.
+ *
+ * In the comprehensive model, each of the LocalMatrixBlocks has
+ * coordinates (0,0) and dimensions equal to the number of degrees of
+ * freedom of the FESystem.
+ *
+ * In the comprehensive model, each block has its own block
+ * coordinates and the size depends on the associated
+ * FESystem::base_element(). These blocks can be generated separately
+ * and will be assembled into the correct matrix block by this object.
+ *
+ * @author Guido Kanschat, 2009
+ */
+    template <class MATRIX, typename number = double>
+    class MatrixLocalBlocksToGlobalBlocks
+    {
+      public:
+                                        /// The object that is stored
+       typedef boost::shared_ptr<MatrixBlock<MATRIX> > MatrixPtr;
+       
+                                        /**
+                                         * Constructor, initializing
+                                         * the #threshold, which
+                                         * limits how small numbers
+                                         * may be to be entered into
+                                         * the matrix.
+                                         */
+       MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
+       
+                                        /**
+                                         * Copy the BlockInfo and the
+                                         * matrix pointers into local
+                                         * variables and initialize
+                                         * cell matrix vectors.
+                                         */
+       void initialize(std::vector<MatrixPtr>& matrices);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face) const;
+       
+
+                                        /**
+                                         * Assemble the local matrices
+                                         * into the global matrices.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble all local matrices
+                                         * into the global matrices.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& info2);
+       
+      private:
+                                        /**
+                                         * Assemble a single local
+                                         * matrix into a global one.
+                                         */
+       void assemble(
+         MATRIX& global,
+         const FullMatrix<number>& local,
+         unsigned int block_row,
+         unsigned int block_col,
+         const std::vector<unsigned int>& dof1,
+         const std::vector<unsigned int>& dof2);
+       
+                                        /**
+                                         * The global matrices,
+                                         * stored as a vector of
+                                         * pointers.
+                                         */
+       std::vector<MatrixPtr> matrices;
+
+                                        /**
+                                         * 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;
+       
+    };
+
+/**
+ * A helper class assembling local matrices into global multilevel
+ * matrices. This class is the multilevel equivalent of
+ * MatrixLocalBlocksToGlobalBlocks and documentation of that class
+ * applies here to a large extend.
+ *
+ * The global matrices are expected as a vector of pointers to MatrixBlock
+ * objects, each containing a MGLevelObject with matrices with a function
+ * corresponding to SparseMatrix::add() and information on the block
+ * row and column this matrix represents in a block system.
+ *
+ * The local matrices are a similar vector of MatrixBlock objects, but
+ * containing a FullMatrix.
+ *
+ * If local refinement occurs, the Multigrid method needs more
+ * matrices, two for continuous elements and another two if numerical
+ * fluxes are computed on interfaces. The second set can be added
+ * using initialize_edge_flux(). Once added, the contributions in all
+ * participating matrices will be assembled from the cell and face
+ * matrices automatically.
+ *
+ * @author Guido Kanschat, 2009
+ */
+    template <class MATRIX, typename number = double>
+    class MGMatrixLocalBlocksToGlobalBlocks
+    {
+      public:
+                                        /// The object that is stored
+       typedef boost::shared_ptr<MatrixBlock<MGLevelObject<MATRIX> > > MatrixPtr;
+
+                                        /**
+                                         * Constructor, initializing
+                                         * the #threshold, which
+                                         * limits how small numbers
+                                         * may be to be entered into
+                                         * the matrix.
+                                         */
+       MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
+       
+                                        /**
+                                         * Copy the BlockInfo and the
+                                         * matrix pointers into local
+                                         * variables and initialize
+                                         * cell matrix vectors.
+                                         */
+       void initialize(std::vector<MatrixPtr>& matrices);
+
+                                        /**
+                                         * Multigrid methods on
+                                         * locally refined meshes
+                                         * need additional
+                                         * matrices. For
+                                         * discontinuous Galerkin
+                                         * methods, these are two
+                                         * flux matrices across the
+                                         * refinement edge, which are
+                                         * set by this method.
+                                         */
+       void initialize_edge_flux(std::vector<MatrixPtr>& up, std::vector<MatrixPtr>& down);
+                                        /**
+                                         * Initialize the local data
+                                         * in the
+                                         * DoFInfo
+                                         * object used later for
+                                         * assembling.
+                                         *
+                                         * The second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info,
+                            bool interior_face) const;
+       
+       
+                                        /**
+                                         * Assemble the local matrices
+                                         * into the global matrices.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble all local matrices
+                                         * into the global matrices.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& info2);
+       
+      private:
+                                        /**
+                                         * Assemble a single local
+                                         * matrix into a global one.
+                                         */
+       void assemble(
+         MATRIX& global,
+         const FullMatrix<number>& local,
+         unsigned int block_row,
+         unsigned int block_col,
+         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.
+                                         */
+       std::vector<MatrixPtr> matrices;
+       
+                                        /**
+                                         * The flux matrix between
+                                         * the fine and the coarse
+                                         * level at refinement edges.
+                                         */
+       std::vector<MatrixPtr> flux_down;
+       
+                                        /**
+                                         * The flux matrix between
+                                         * the coarse and the fine
+                                         * level at refinement edges.
+                                         */
+       std::vector<MatrixPtr> flux_up;
+       
+                                        /**
+                                         * 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.
+ *
+ * @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 second parameter is
+                                         * used to distinguish
+                                         * between the data used on
+                                         * cells and boundary faces
+                                         * on the one hand and
+                                         * interior faces on the
+                                         * other. Interior faces may
+                                         * require additional data
+                                         * being initialized.
+                                         */
+       template <int dim>
+       void initialize_info(DoFInfo<dim>& info, bool) const;
+       
+                                        /**
+                                         * Assemble the matrix
+                                         * DoFInfo::M1[0]
+                                         * into the global matrix.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info);
+       
+                                        /**
+                                         * Assemble both local
+                                         * matrices in the info
+                                         * objects into the global
+                                         * matrix.
+                                         */
+       template<int dim>
+       void assemble(const DoFInfo<dim>& info1,
+                     const DoFInfo<dim>& info2);
+    };
+
+    
+//----------------------------------------------------------------------//
+
+    template <class VECTOR>
+    inline void
+    ResidualSimple<VECTOR>::initialize(NamedData<SmartPointer<VECTOR> >& results)
+    {
+      residuals = results;
+    }
+
+    
+    template <class VECTOR>
+    template<int dim>
+    inline void
+    ResidualSimple<VECTOR>::initialize_info(DoFInfo<dim>& info, bool) const
+    {
+      info.initialize_vectors(residuals.size());
+    }
+
+    
+    template <class VECTOR>
+    template<int dim>
+    inline void
+    ResidualSimple<VECTOR>::assemble(const DoFInfo<dim>& info)
+    {
+      for (unsigned int k=0;k<residuals.size();++k)
+       for (unsigned int i=0;i<info.R[k].block(0).size();++i)
+         (*residuals(k))(info.indices[i]) += info.R[k].block(0)(i);
+    }
+
+    
+    template <class VECTOR>
+    template<int dim>
+    inline void
+    ResidualSimple<VECTOR>::assemble(const DoFInfo<dim>& info1,
+                                       const DoFInfo<dim>& info2)
+    {
+      for (unsigned int k=0;k<residuals.size();++k)
+       {
+         for (unsigned int i=0;i<info1.R[k].block(0).size();++i)
+           (*residuals(k))(info1.indices[i]) += info1.R[k].block(0)(i);
+         for (unsigned int i=0;i<info2.R[k].block(0).size();++i)
+           (*residuals(k))(info2.indices[i]) += info2.R[k].block(0)(i);
+       }
+    }
+
+    
+//----------------------------------------------------------------------//
+
+    template <class VECTOR>
+    inline void
+    ResidualLocalBlocksToGlobalBlocks<VECTOR>::initialize(NamedData<SmartPointer<VECTOR> >& m)
+    {
+      residuals = m;      
+    }
+
+
+    template <class VECTOR>
+    template <int dim>
+    inline void
+    ResidualLocalBlocksToGlobalBlocks<VECTOR>::initialize_info(
+      DoFInfo<dim>& info, bool) const
+    {
+      info.initialize_vectors(residuals.size());
+    }
+
+    template <class VECTOR>
+    inline void
+    ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
+      VECTOR global,
+      const BlockVector<double>& local,
+      const BlockIndices& bi,
+      const std::vector<unsigned int>& dof)
+    {
+      for (unsigned int b=0;b<local.n_blocks();++b)
+       for (unsigned int j=0;j<local.block(b).size();++j)
+         {
+                                            // The coordinates of
+                                            // the current entry in
+                                            // DoFHandler
+                                            // numbering, which
+                                            // differs from the
+                                            // block-wise local
+                                            // numbering we use in
+                                            // our local vectors
+           const unsigned int jcell = this->bi.local_to_global(b, j);
+           (*global)(dof[jcell]) += local.block(b)(j);
+         }
+    }
+
+    
+    template <class VECTOR>
+    template <int dim>
+    inline void
+    ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
+      const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<residuals.n_vectors();++i)
+       assemble(residuals.vector(i), info.R[i], info.block_info.local(), info.indices);
+    }
+
+    
+    template <class VECTOR>
+    template <int dim>
+    inline void
+    ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
+      const DoFInfo<dim>& info1,
+      const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<residuals.n_vectors();++i)
+       {
+         assemble(residuals.vector(i), info1.R[i], info1.block_info.local(), info1.indices);
+         assemble(residuals.vector(i), info2.R[i], info2.block_info.local(), info2.indices);
+       }
+    }
+
+
+//----------------------------------------------------------------------//
+
+    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 >
+    template <int dim>
+    inline void
+    MatrixSimple<MATRIX>::initialize_info(DoFInfo<dim>& info,
+                                         bool interior_face) const
+    {
+      info.initialize_matrices(1, interior_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());
+      
+      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));
+    }
+    
+    
+    template <class MATRIX>
+    template <int dim>
+    inline void
+    MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info)
+    {
+      assemble(info.M1[0].matrix, info.indices, info.indices);
+    }
+    
+
+    template <class MATRIX>
+    template <int dim>
+    inline void
+    MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info1,
+                                  const DoFInfo<dim>& info2)
+    {
+      assemble(info1.M1[0].matrix, info1.indices, info1.indices);
+      assemble(info1.M2[0].matrix, info1.indices, info2.indices);
+      assemble(info2.M1[0].matrix, info2.indices, info2.indices);
+      assemble(info2.M2[0].matrix, info2.indices, info1.indices);
+    }
+    
+    
+//----------------------------------------------------------------------//
+    
+    template <class MATRIX, typename number>
+    inline
+    MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MatrixLocalBlocksToGlobalBlocks(
+      double threshold)
+                   :
+                   threshold(threshold)
+    {}
+    
+    
+    template <class MATRIX, typename number>
+    inline void
+    MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(std::vector<MatrixPtr>& m)
+    {
+      matrices = m;
+    }
+    
+
+    
+    template <class MATRIX, typename number>
+    inline void
+    MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      MATRIX& global,
+      const FullMatrix<number>& local,
+      unsigned int block_row,
+      unsigned int block_col,
+      const std::vector<unsigned int>& dof1,
+      const std::vector<unsigned int>& dof2)
+    {
+      for (unsigned int j=0;j<local.n_rows();++j)
+       for (unsigned int k=0;k<local.n_cols();++k)
+         if (std::fabs(local(j,k)) >= threshold)
+           {
+                                              // The coordinates of
+                                              // the current entry in
+                                              // DoFHandler
+                                              // numbering, which
+                                              // differs from the
+                                              // block-wise local
+                                              // numbering we use in
+                                              // our local matrices
+             const unsigned int jcell = this->block_info.local.local_to_global(block_row, j);
+             const unsigned int kcell = this->block_info.local.local_to_global(block_col, k);
+
+                                              // The global dof
+                                              // indices to assemble
+                                              // in. Since we may
+                                              // have face matrices
+                                              // coupling two
+                                              // different cells, we
+                                              // provide two sets of
+                                              // dof indices.
+             const unsigned int jglobal = this->block_info.global.global_to_local(dof1[jcell]).second;
+             const unsigned int kglobal = this->block_info.global.global_to_local(dof2[kcell]).second;
+             
+             global.add(jglobal, kglobal, local(j,k));
+           }
+    }
+
+    
+    template <class MATRIX, typename number>
+    template <int dim>
+    inline void
+    MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<matrices.size();++i)
+       {
+                                          // Row and column index of
+                                          // the block we are dealing with
+         const unsigned int row = matrices[i]->row;
+         const unsigned int col = matrices[i]->column;
+
+         assemble(matrices[i]->matrix, info->M1[i].matrix, row, col, info.indices, info.indices);
+       }
+    }
+
+    
+    template <class MATRIX, typename number>
+    template <int dim>
+    inline void
+    MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      const DoFInfo<dim>& info1,
+      const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<matrices.size();++i)
+       {
+                                          // Row and column index of
+                                          // the block we are dealing with
+         const unsigned int row = matrices[i]->row;
+         const unsigned int col = matrices[i]->column;
+
+         assemble(matrices[i]->matrix, this->M11[i].matrix, row, col, info1.indices, info1.indices);
+         assemble(matrices[i]->matrix, this->M12[i].matrix, row, col, info1.indices, info2.indices);
+         assemble(matrices[i]->matrix, this->M21[i].matrix, row, col, info2.indices, info1.indices);
+         assemble(matrices[i]->matrix, this->M22[i].matrix, row, col, info2.indices, info2.indices);
+
+         this->M11[i].matrix = 0.;
+         this->M12[i].matrix = 0.;
+         this->M21[i].matrix = 0.;
+         this->M22[i].matrix = 0.;
+       }
+    }
+
+
+// ----------------------------------------------------------------------//
+    
+    template <class MATRIX, typename number>
+    inline
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MGMatrixLocalBlocksToGlobalBlocks(
+      double threshold)
+                   :
+                   threshold(threshold)
+    {}
+    
+    
+    template <class MATRIX, typename number>
+    inline void
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(std::vector<MatrixPtr>& m)
+    {
+      matrices = m;
+    }
+    
+
+    template <class MATRIX, typename number>
+    inline void
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_edge_flux(
+      std::vector<MatrixPtr>& up,
+      std::vector<MatrixPtr>& down)
+    {
+      flux_up = up;
+      flux_down = down;
+    }
+    
+    
+    template <class MATRIX, typename number>
+    inline void
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      MATRIX& global,
+      const FullMatrix<number>& local,
+      unsigned int block_row,
+      unsigned int block_col,
+      const std::vector<unsigned int>& dof1,
+      const std::vector<unsigned int>& dof2,
+      unsigned int level1,
+      unsigned int level2,
+      bool transpose)
+    {
+      for (unsigned int j=0;j<local.n_rows();++j)
+       for (unsigned int k=0;k<local.n_cols();++k)
+         if (std::fabs(local(j,k)) >= threshold)
+           {
+                                              // The coordinates of
+                                              // the current entry in
+                                              // DoFHandler
+                                              // numbering, which
+                                              // differs from the
+                                              // block-wise local
+                                              // numbering we use in
+                                              // our local matrices
+             const unsigned int jcell = this->block_info.local.local_to_global(block_row, j);
+             const unsigned int kcell = this->block_info.local.local_to_global(block_col, k);
+
+                                              // The global dof
+                                              // indices to assemble
+                                              // in. Since we may
+                                              // have face matrices
+                                              // coupling two
+                                              // different cells, we
+                                              // provide two sets of
+                                              // dof indices.
+             const unsigned int jglobal = this->block_info.levels[level1].global_to_local(dof1[jcell]).second;
+             const unsigned int kglobal = this->block_info.levels[level2].global_to_local(dof2[kcell]).second;
+
+             if (transpose)
+               global.add(kglobal, jglobal, local(j,k));
+             else
+               global.add(jglobal, kglobal, local(j,k));
+           }
+    }
+
+    
+    template <class MATRIX, typename number>
+    template <int dim>
+    inline void
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      const DoFInfo<dim>& info)
+    {
+      const unsigned int level = info.cell->level();
+      
+      for (unsigned int i=0;i<matrices.size();++i)
+       {
+                                          // Row and column index of
+                                          // the block we are dealing with
+         const unsigned int row = matrices[i]->row;
+         const unsigned int col = matrices[i]->column;
+
+         assemble(matrices[i]->matrix[level], this->M11[i].matrix, row, col,
+                  info.indices, info.indices, level, level);
+         
+         this->M11[i].matrix = 0.;
+       }
+    }
+
+    
+    template <class MATRIX, typename number>
+    template <int dim>
+    inline void
+    MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
+      const DoFInfo<dim>& info1,
+      const DoFInfo<dim>& info2)
+    {
+      const unsigned int level1 = info1.cell->level();
+      const unsigned int level2 = info2.cell->level();
+      
+      for (unsigned int i=0;i<matrices.size();++i)
+       {
+                                          // Row and column index of
+                                          // the block we are dealing with
+         const unsigned int row = matrices[i]->row;
+         const unsigned int col = matrices[i]->column;
+
+         if (level1 == level2)
+           {
+             assemble(matrices[i]->matrix[level1], this->M11[i].matrix, row, col, info1.indices, info1.indices, level1, level1);
+             assemble(matrices[i]->matrix[level1], this->M12[i].matrix, row, col, info1.indices, info2.indices, level1, level2);
+             assemble(matrices[i]->matrix[level1], this->M21[i].matrix, row, col, info2.indices, info1.indices, level2, level1);
+             assemble(matrices[i]->matrix[level1], this->M22[i].matrix, row, col, info2.indices, info2.indices, level2, level2);
+           }
+         else
+           {
+             Assert(level1 > level2, ExcNotImplemented());
+             if (flux_up.size() != 0)
+               {
+                                                  // Do not add M22,
+                                                  // which is done by
+                                                  // the coarser cell
+                 assemble(matrices[i]->matrix[level1], this->M11[i].matrix, row, col,
+                          info1.indices, info1.indices, level1, level1);             
+                 assemble(flux_up[i]->matrix[level1], this->M12[i].matrix, row, col,
+                          info1.indices, info2.indices, level1, level2, true);
+                 assemble(flux_down[i]->matrix[level1], this->M21[i].matrix, row, col,
+                          info2.indices, info1.indices, level2, level1);
+               }
+           }
+         
+         this->M11[i].matrix = 0.;
+         this->M12[i].matrix = 0.;
+         this->M21[i].matrix = 0.;
+         this->M22[i].matrix = 0.;
+       }
+    }
+
+//----------------------------------------------------------------------//
+
+    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<SmartPointer<VECTOR> > data;
+      SmartPointer<VECTOR> p = &rhs;
+      data.add(p, "right hand side");
+      
+      MatrixSimple<MATRIX>::initialize(m);
+      ResidualSimple<VECTOR>::initialize(data);
+    }
+
+    
+    template <class MATRIX, class VECTOR>
+    template <int dim>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::initialize_info(DoFInfo<dim>& info,
+                                                bool interior_face) const
+    {
+      MatrixSimple<MATRIX>::initialize_info(info, interior_face);
+      ResidualSimple<VECTOR>::initialize_info(info, interior_face);
+    }
+
+
+    template <class MATRIX, class VECTOR>
+    template<int dim>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::assemble(const DoFInfo<dim>& info)
+    {
+      MatrixSimple<MATRIX>::assemble(info);
+      ResidualSimple<VECTOR>::assemble(info);
+    }
+
+
+    template <class MATRIX, class VECTOR>
+    template<int dim>
+    inline void
+    SystemSimple<MATRIX,VECTOR>::assemble(const DoFInfo<dim>& info1,
+                                         const DoFInfo<dim>& info2)
+    {
+      MatrixSimple<MATRIX>::assemble(info1, info2);
+      ResidualSimple<VECTOR>::assemble(info1, info2);
+    }    
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.h b/deal.II/deal.II/include/numerics/mesh_worker_info.h
new file mode 100644 (file)
index 0000000..c0e8e52
--- /dev/null
@@ -0,0 +1,874 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
+//
+//    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_info_h
+#define __deal2__mesh_worker_info_h
+
+#include <base/config.h>
+#include <lac/matrix_block.h>
+#include <boost/shared_ptr.hpp>
+#include <dofs/block_info.h>
+#include <fe/fe_values.h>
+#include <numerics/mesh_worker_vector_selector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+/**
+ * The class providing the scrapbook to fill with local integration
+ * results. These can be values, local contributions to forms or cell
+ * and face matrices.
+ *
+ * The local matrices initialized by reinit() of the info object and
+ * then assembled into the global system by Assembler classes.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template <typename number>
+  class LocalResults
+  {
+                                      /**
+                                       * Initialize a single local
+                                       * matrix block. A helper
+                                       * function for initialize()
+                                       */
+      void initialize_local(
+       MatrixBlock<FullMatrix<number> >& M,
+       unsigned int row, unsigned int col);
+    public:
+      void initialize_vectors(const unsigned int n_vectors);
+                                      /**
+                                       * Allocate @p n local
+                                       * matrices. Additionally,
+                                       * set their block row and
+                                       * column coordinates to
+                                       * zero. The matrices
+                                       * themselves are resized by
+                                       * reinit().
+                                       *
+                                       * The template parameter @p
+                                       * MatrixPtr should point to
+                                       * a MatrixBlock
+                                       * instantiation in order to
+                                       * provide row and column info.
+                                       */
+      void initialize_matrices(unsigned int n, bool both);
+       
+                                      /**
+                                       * Allocate a local matrix
+                                       * for each of the global
+                                       * ones in @p
+                                       * matrices. Additionally,
+                                       * set their block row and
+                                       * column coordinates. The
+                                       * matrices themselves are
+                                       * resized by reinit().
+                                       *
+                                       * The template parameter @p
+                                       * MatrixPtr should point to
+                                       * a MatrixBlock
+                                       * instantiation in order to
+                                       * provide row and column info.
+                                       */
+      template <class MatrixPtr>
+      void initialize_matrices(std::vector<MatrixPtr>& matrices,
+                              bool both);
+       
+                                      /**
+                                       * Reinitialize matrices for
+                                       * new cell. Resizes the
+                                       * matrices for hp and sets
+                                       * them to zero.
+                                       */
+      void reinit(const BlockIndices& local_sizes);
+       
+                                      /**
+                                       * The local numbers,
+                                       * computed on a cell or on a
+                                       * face.
+                                       */
+      std::vector<number> J;
+       
+                                      /**
+                                       * The local vectors. This
+                                       * field is public, so that
+                                       * local integrators can
+                                       * write it.
+                                       */
+      std::vector<BlockVector<number> > R;
+       
+                                      /**
+                                       * The local matrices
+                                       * coupling degrees of
+                                       * freedom in the cell
+                                       * itself or within the
+                                       * first cell on a face.
+                                       */
+      std::vector<MatrixBlock<FullMatrix<number> > > M1;
+       
+                                      /**
+                                       * The local matrices
+                                       * coupling test functions on
+                                       * the cell with trial
+                                       * functions on the other
+                                       * cell.
+                                       *
+                                       * Only used on interior
+                                       * faces.
+                                       */
+      std::vector<MatrixBlock<FullMatrix<number> > > M2;
+  };
+
+    
+/**
+ * Basic info class only containing information on geometry and
+ * degrees of freedom of the mesh object.
+ *
+ * The information in these objects is usually used by one of the
+ * Assembler classes. It is also the kind of information which is
+ * needed in mesh based matrices (often referred to as matrix free
+ * methods).
+ *
+ * In addition to the information on degrees of freedom stored in this
+ * class, it also provides the local computation space for the worker
+ * object operating on it. This space is provided by the base class
+ * template DATATYPE. This base class will automatically
+ * reinitialized on each cell, but initial setup is up to the user and
+ * should be done when initialize() for this class is called. The
+ * currently available base classes are
+ * <ul>
+ * <li> LocalVectors
+ * <li> LocalMatrices
+ * </ul>
+ *
+ * This class operates in two different modes, corresponding to the
+ * data models discussed in the Assembler namespace documentation.
+ *
+ * The choice of the local data model is triggered by the vector
+ * #BlockInfo::local_renumbering, which in turn is usually filled by
+ * BlockInfo::initialize_local(). If this function has been used, or
+ * the vector has been changed from zero-length, then local dof
+ * indices stored in this object will automatically be renumbered to
+ * reflect local block structure.
+ *
+ * The BlockInfo object is stored as a pointer. Therefore, if the
+ * block structure changes, for instance because of mesh refinement,
+ * the DoFInfo class will automatically use the new structures.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template<int dim, int spacedim = dim>
+  class DoFInfo : public LocalResults<double>
+  {
+    public:    
+                                      /// The current cell
+      typename Triangulation<dim>::cell_iterator cell;
+       
+                                      /// The current face
+      typename Triangulation<dim>::face_iterator face;
+       
+                                      /**
+                                       * The number of the current
+                                       * face on the current cell.
+                                       *
+                                       * This number is
+                                       * deal_II_numbers::invalid_unsigned_int
+                                       * if the info object was
+                                       * initialized with a cell.
+                                       */
+                                         
+      unsigned int face_number;
+                                      /**
+                                       * The number of the current
+                                       * subface on the current
+                                       * face
+                                       *
+                                       * This number is
+                                       * deal_II_numbers::invalid_unsigned_int
+                                       * if the info object was not
+                                       * initialized with a subface.
+                                       */
+      unsigned int sub_number;
+
+                                      /// The DoF indices of the current cell
+      std::vector<unsigned int> indices;
+       
+                                      /**
+                                       * Constructor setting the
+                                       * #block_info pointer.
+                                       */
+      DoFInfo(const BlockInfo& block_info);
+
+                                      /**
+                                       * Default constructor
+                                       * leaving the #block_info
+                                       * pointer empty, but setting
+                                       * the #aux_local_indices.
+                                       */
+      template <class DH>
+      DoFInfo(const DH& dof_handler);
+       
+                                      /**
+                                       * Set the current cell and
+                                       * fill #indices.
+                                       */
+      template <class DHCellIterator>
+      void reinit(const DHCellIterator& c);
+
+                                      /**
+                                       * Set the current face and
+                                       * fill #indices if the #cell
+                                       * changed.
+                                       */
+      template <class DHCellIterator, class DHFaceIterator>
+      void reinit(const DHCellIterator& c,
+                 const DHFaceIterator& f,
+                 unsigned int n);
+       
+                                      /**
+                                       * Set the current subface
+                                       * and fill #indices if the
+                                       * #cell changed.
+                                       */
+      template <class DHCellIterator, class DHFaceIterator>
+      void reinit(const DHCellIterator& c,
+                 const DHFaceIterator& f,
+                 const unsigned int n, const unsigned int s);
+
+      const BlockIndices& local_indices() const;
+      
+       
+                                      /// The block structure of the system
+      SmartPointer<const BlockInfo> block_info;
+    private:
+                                      /// Fill index vector
+      void get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator c);
+       
+                                      /// Fill index vector with level indices
+      void get_indices(const typename MGDoFHandler<dim, spacedim>::cell_iterator c);
+       
+                                      /// Auxiliary vector
+      std::vector<unsigned int> indices_org;
+       
+                                      /**
+                                       * An auxiliary local
+                                       * BlockIndices object created
+                                       * if #block_info is not set.
+                                       * It contains just a single
+                                       * block of the size of
+                                       * degrees of freedom per cell.
+                                       */
+      BlockIndices aux_local_indices;
+  };
+    
+/**
+ * Class for objects handed to local integration functions.
+ *
+ * Objects of this class contain one or more objects of type FEValues,
+ * FEFaceValues or FESubfacevalues to be used in local
+ * integration. They are stored in an array of pointers to the base
+ * classes FEValuesBase for cells and FEFaceValuesBase for faces and
+ * subfaces, respectively. The template parameter VECTOR allows the
+ * use of different data types for the global system.
+ *
+ * The @p FEVALUESBASE template parameter should be either
+ * FEValuesBase or FEFaceValuesBase, depending on whether the object
+ * is used to integrate over cells or faces. The actual type of @p
+ * FEVALUES object is fixed in the constructor and only used to
+ * initialize the pointers in #fevalv.
+ *
+ * Additionally, this function containes space to store the values of
+ * finite element functions stored in #global_data in the
+ * quadrature points. These vectors are initialized automatically on
+ * each cell or face. In order to avoid initializing unused vectors,
+ * you can use initialize_selector() to select the vectors by name
+ * that you actually want to use.
+ *
+ * <h3>Integration models</h3>
+ *
+ * This class supports two local integration models, corresponding to
+ * the data models in the documentation of the Assembler namespace.
+ * One is the
+ * standard model suggested by the use of FESystem. Namely, there is
+ * one FEValuseBase object in this class, containing all shape
+ * functions of the whole system, and having as many components as the
+ * system. Using this model involves loops over all system shape
+ * functions. It requires to identify the system components
+ * for each shape function and to select the correct bilinear form,
+ * usually in an @p if or @p switch statement.
+ *
+ * The second integration model builds one FEValuesBase object per
+ * base element of the system. The degrees of freedom on each cell are
+ * renumbered by block, such that they represent the same block
+ * structure as the global system. Objects performing the integration
+ * can then process each block separately, which improves reusability
+ * of code considerably.
+ *
+ * @note As described in DoFInfo, the use of the local block model is
+ * triggered by calling BlockInfo::initialize_local() before
+ * using initialize() in this class.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template<int dim, class FEVALUESBASE, int spacedim = dim>
+  class IntegrationInfo : public DoFInfo<dim, spacedim>
+  {
+                                      /// vector of FEValues objects
+      std::vector<boost::shared_ptr<FEVALUESBASE> > fevalv;
+    public:
+                                      /**
+                                       * Constructor forwarding
+                                       * information to DoFInfo.
+                                       */
+      IntegrationInfo(const BlockInfo& block_info);
+       
+                                      /**
+                                       * Constructor forwarding
+                                       * information to DoFInfo.
+                                       */
+      template <class DH>
+      IntegrationInfo(const DH& dof_handler);
+       
+                                      /**
+                                       * Build all internal
+                                       * structures, in particular
+                                       * the FEValuesBase objects
+                                       * and allocate space for
+                                       * data vectors.
+                                       *
+                                       * @param fe is only used
+                                       * to determine the type of
+                                       * the element and can be a
+                                       * null pointer to FEValues,
+                                       * FEVaceValues or
+                                       * FESubfaceValues.
+                                       *
+                                       * @param el is the finite
+                                       * element of the DoFHandler.
+                                       *
+                                       * @param mapping is the Mapping
+                                       * object used to map the
+                                       * mesh cells.
+                                       *
+                                       * @param quadrature is a
+                                       * Quadrature formula used in
+                                       * the constructor of the
+                                       * FEVALUES objects.
+                                       *
+                                       * @param flags are the
+                                       * UpdateFlags used in
+                                       * the constructor of the
+                                       * FEVALUES objects.
+                                       */
+      template <class FEVALUES>
+      void initialize(const FiniteElement<dim,spacedim>& el,
+                     const Mapping<dim,spacedim>& mapping,
+                     const Quadrature<FEVALUES::integral_dimension>& quadrature,
+                     const UpdateFlags flags);
+
+                                      /**
+                                       * Initialize the data
+                                       * vector and cache the
+                                       * selector.
+                                       */
+      void initialize_data(const boost::shared_ptr<VectorDataBase<dim,spacedim> > data);
+       
+                                      /**
+                                       * Delete the data created by initialize().
+                                       */
+      void clear();
+       
+                                      /// This is true if we are assembling for multigrid
+      bool multigrid;
+                                      /// Access to finite element
+                                      /**
+                                       * This is the access
+                                       * function being used, if
+                                       * the constructor for a
+                                       * single element was
+                                       * used. It throws an
+                                       * exception, if applied to a
+                                       * vector of elements.
+                                       */
+      const FEVALUESBASE& fe() const;
+       
+                                      /// Access to finite elements
+                                      /**
+                                       * This access function must
+                                       * be used if the constructor
+                                       * for a group of elements
+                                       * was used.
+                                       *
+                                       * @see DGBlockSplitApplication
+                                       */
+      const FEVALUESBASE& fe(unsigned int i) const;
+       
+                                      /**
+                                       * The vector containing the
+                                       * values of finite element
+                                       * functions in the quadrature
+                                       * points.
+                                       *
+                                       * There is one vector per
+                                       * selected finite element
+                                       * function, containing one
+                                       * vector for each component,
+                                       * containing vectors with
+                                       * values for each quadrature
+                                       * point.
+                                       */
+     std::vector<std::vector<std::vector<double> > > values;
+       
+                                      /**
+                                       * The vector containing the
+                                       * derivatives of finite
+                                       * element functions in the
+                                       * quadrature points.
+                                       *
+                                       * There is one vector per
+                                       * selected finite element
+                                       * function, containing one
+                                       * vector for each component,
+                                       * containing vectors with
+                                       * values for each quadrature
+                                       * point.
+                                       */
+      std::vector<std::vector<std::vector<Tensor<1,dim> > > > gradients;
+      
+                                      /**
+                                       * The vector containing the
+                                       * second derivatives of finite
+                                       * element functions in the
+                                       * quadrature points.
+                                       *
+                                       * There is one vector per
+                                       * selected finite element
+                                       * function, containing one
+                                       * vector for each component,
+                                       * containing vectors with
+                                       * values for each quadrature
+                                       * point.
+                                       */
+      std::vector<std::vector<std::vector<Tensor<2,dim> > > > hessians;
+       
+                                      /**
+                                       * Reinitialize internal data
+                                       * structures for use on a cell.
+                                       */
+      template <class DHCellIterator>
+      void reinit(const DHCellIterator& c);
+
+                                      /**
+                                       * Reinitialize internal data
+                                       * structures for use on a face.
+                                       */
+      template <class DHCellIterator, class DHFaceIterator>
+      void reinit(const DHCellIterator& c,
+                 const DHFaceIterator& f,
+                 const unsigned int fn);
+       
+                                      /**
+                                       * Reinitialize internal data
+                                       * structures for use on a subface.
+                                       */
+      template <class DHCellIterator, class DHFaceIterator>
+      void reinit(const DHCellIterator& c,
+                 const DHFaceIterator& f,
+                 const unsigned int fn,
+                 const unsigned int sn);
+
+
+                                      /**
+                                       * @deprecated This is the
+                                       * old version not using
+                                       * VectorSelector.
+                                       *
+                                       * Use the finite element
+                                       * functions in #global_data
+                                       * and fill the vectors
+                                       * #values, #gradients and
+                                       * #hessians.
+                                       */
+      void fill_local_data(bool split_fevalues);
+
+                                      /**
+                                       * The global data vector
+                                       * used to compute function
+                                       * values in quadrature
+                                       * points.
+                                       */
+      boost::shared_ptr<VectorDataBase<dim, spacedim> > global_data;
+    private:
+                                      /**
+                                       * Use the finite element
+                                       * functions in #global_data
+                                       * and fill the vectors
+                                       * #values, #gradients and
+                                       * #hessians with values
+                                       * according to the
+                                       * selector.
+                                       */
+      template <typename TYPE>
+      void fill_local_data(
+       std::vector<std::vector<std::vector<TYPE> > >& data,
+       VectorSelector& selector,
+       bool split_fevalues) const;
+       
+  };
+
+/**
+ * A simple container collecting the five info objects required by the
+ * integration loops.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template <int dim, int spacedim=dim>
+  class IntegrationInfoBox
+  {
+    public:
+      typedef IntegrationInfo<dim, FEValuesBase<dim, spacedim>, spacedim> CellInfo;
+      typedef IntegrationInfo<dim, FEFaceValuesBase<dim, spacedim>, spacedim> FaceInfo;
+
+                                      /**
+                                       * Initialize all members
+                                       * with the same argument.
+                                       */
+      template <typename T>
+      IntegrationInfoBox(const T&);
+       
+      template <class WORKER>
+      void initialize(const WORKER&,
+                     const FiniteElement<dim, spacedim>& el,
+                     const Mapping<dim, spacedim>& mapping);
+       
+      template <class WORKER, class VECTOR>
+      void initialize(const WORKER&,
+                     const FiniteElement<dim, spacedim>& el,
+                     const Mapping<dim, spacedim>& mapping,
+                     const NamedData<SmartPointer<VECTOR> >& data);
+//      private:
+       
+      CellInfo cell_info;
+      FaceInfo bdry_info;
+      FaceInfo face_info;
+      FaceInfo subface_info;
+      FaceInfo neighbor_info;
+  };
+    
+
+//----------------------------------------------------------------------//
+
+  template <typename number>
+  inline void
+  LocalResults<number>::initialize_vectors(unsigned int n)
+  {
+    R.resize(n);
+  }
+
+    
+  template <typename number>
+  template <class MatrixPtr>
+  inline void
+  LocalResults<number>::initialize_matrices(
+    std::vector<MatrixPtr>& matrices,
+    bool both)
+  {
+    M1.resize(matrices.size());
+    if (both)
+      M2.resize(matrices.size());
+    for (unsigned int i=0;i<matrices.size();++i)
+      {
+       const unsigned int row = matrices[i]->row;
+       const unsigned int col = matrices[i]->column;
+
+       M1[i].row = row;
+       M1[i].column = col;
+       if (both)
+         {
+           M2[i].row = row;
+           M2[i].column = col;
+         }
+      }
+  }
+    
+    
+  template <typename number>
+  inline void
+  LocalResults<number>::initialize_matrices(unsigned int n, bool both)
+  {
+    M1.resize(n);
+    if (both)
+      M2.resize(n);
+    for (unsigned int i=0;i<n;++i)
+      {
+       M1[i].row = 0;
+       M1[i].column = 0;
+       if (both)
+         {
+           M2[i].row = 0;
+           M2[i].column = 0;
+         }
+      }
+  }
+    
+    
+  template <typename number>
+  inline void
+  LocalResults<number>::reinit(const BlockIndices& bi)
+  {
+    for (unsigned int i=0;i<R.size();++i)
+      R[i].reinit(bi);
+    for (unsigned int i=0;i<M1.size();++i)
+      M1[i].matrix.reinit(bi.block_size(M1[i].row),
+                         bi.block_size(M1[i].column));
+    for (unsigned int i=0;i<M2.size();++i)
+      M2[i].matrix.reinit(bi.block_size(M2[i].row),
+                         bi.block_size(M2[i].column));
+  }   
+   
+
+//----------------------------------------------------------------------//
+    
+  template <int dim, int spacedim>
+  template <class DH>
+  DoFInfo<dim,spacedim>::DoFInfo(const DH& dof_handler)
+  {
+    std::vector<unsigned int> aux(1);
+    aux[0] = dof_handler.get_fe().dofs_per_cell;
+    aux_local_indices.reinit(aux);
+  }
+  
+  
+  template <int dim, int spacedim>
+  template <class DHCellIterator>
+  inline void
+  DoFInfo<dim,spacedim>::reinit(const DHCellIterator& c)
+  {
+    get_indices(c);
+    cell = static_cast<typename Triangulation<dim,spacedim>::cell_iterator> (c);
+    face_number = deal_II_numbers::invalid_unsigned_int;
+    sub_number = deal_II_numbers::invalid_unsigned_int;
+    if (block_info)
+      LocalResults<double>::reinit(block_info->local());
+    else
+      LocalResults<double>::reinit(aux_local_indices); 
+  }
+  
+  
+  template<int dim, int spacedim>
+  template <class DHCellIterator, class DHFaceIterator>
+  inline void
+  DoFInfo<dim, spacedim>::reinit(const DHCellIterator& c,
+                                const DHFaceIterator& f,
+                                unsigned int n)
+  {  
+    if ((cell.state() != IteratorState::valid)
+       ||  cell != static_cast<typename Triangulation<dim>::cell_iterator> (c))
+      get_indices(c);
+    cell = static_cast<typename Triangulation<dim>::cell_iterator> (c);
+    face = static_cast<typename Triangulation<dim>::face_iterator> (f);
+    face_number = n;
+    sub_number = deal_II_numbers::invalid_unsigned_int;
+    if (block_info)
+      LocalResults<double>::reinit(block_info->local());
+    else
+      LocalResults<double>::reinit(aux_local_indices); 
+  }
+  
+  
+  template<int dim, int spacedim>
+  template <class DHCellIterator, class DHFaceIterator>
+  inline void
+  DoFInfo<dim, spacedim>::reinit(const DHCellIterator& c,
+                                const DHFaceIterator& f,
+                                unsigned int n,
+                                const unsigned int s)
+  {
+    if (cell.state() != IteratorState::valid
+       || cell != static_cast<typename Triangulation<dim>::cell_iterator> (c))
+      get_indices(c);
+    cell = static_cast<typename Triangulation<dim>::cell_iterator> (c);
+    face = static_cast<typename Triangulation<dim>::face_iterator> (f);
+    face_number = n;
+    sub_number = s;
+    if (block_info)
+      LocalResults<double>::reinit(block_info->local());
+    else
+      LocalResults<double>::reinit(aux_local_indices); 
+  }
+  
+  
+  template<int dim, int spacedim>
+  inline const BlockIndices&
+  DoFInfo<dim, spacedim>::local_indices() const
+  {
+    if (block_info)
+      return block_info->local();
+    return aux_local_indices;
+  }
+    
+
+//----------------------------------------------------------------------//
+    
+  template <int dim, class FVB, int spacedim>
+  template <class DH>
+  IntegrationInfo<dim,FVB,spacedim>::IntegrationInfo(const DH& dof_handler)
+                 :
+                 DoFInfo<dim, spacedim>(dof_handler),
+                 fevalv(0),
+                 multigrid(false),
+                 global_data(boost::shared_ptr<VectorDataBase<dim, spacedim> >(new VectorDataBase<dim, spacedim>))
+  {}
+
+    
+  template <int dim, class FVB, int spacedim>
+  inline const FVB&
+  IntegrationInfo<dim,FVB,spacedim>::fe() const
+  {
+    AssertDimension(fevalv.size(), 1);
+    return *fevalv[0];
+  }
+
+
+  template <int dim, class FVB, int spacedim>
+  inline const FVB&
+  IntegrationInfo<dim,FVB,spacedim>::fe(unsigned int i) const
+  {
+    Assert (i<fevalv.size(), ExcIndexRange(i,0,fevalv.size()));
+    return *fevalv[i];
+  }
+
+
+  template <int dim, class FVB, int spacedim>
+  template <class DHCellIterator>
+  inline void
+  IntegrationInfo<dim,FVB,spacedim>::reinit(const DHCellIterator& c)
+  {
+    DoFInfo<dim,spacedim>::reinit(c);
+    for (unsigned int i=0;i<fevalv.size();++i)
+      {
+       FVB& febase = *fevalv[i];
+       FEValues<dim>& fe = dynamic_cast<FEValues<dim>&> (febase);
+       fe.reinit(this->cell);
+      }
+
+    const bool split_fevalues = this->block_info != 0;
+    fill_local_data(split_fevalues);
+  }
+    
+    
+  template <int dim, class FVB, int spacedim>
+  template <class DHCellIterator, class DHFaceIterator>
+  inline void
+  IntegrationInfo<dim,FVB,spacedim>::reinit(
+    const DHCellIterator& c,
+    const DHFaceIterator& f,
+    const unsigned int fn)
+  {
+    DoFInfo<dim,spacedim>::reinit(c, f, fn);
+    for (unsigned int i=0;i<fevalv.size();++i)
+      {
+       FVB& febase = *fevalv[i];
+       FEFaceValues<dim>& fe = dynamic_cast<FEFaceValues<dim>&> (febase);
+       fe.reinit(this->cell, fn);
+      }
+      
+    const bool split_fevalues = this->block_info != 0;
+    fill_local_data(split_fevalues);
+  }
+    
+    
+  template <int dim, class FVB, int spacedim>
+  template <class DHCellIterator, class DHFaceIterator>
+  inline void
+  IntegrationInfo<dim,FVB,spacedim>::reinit(
+    const DHCellIterator& c,
+    const DHFaceIterator& f,
+    const unsigned int fn,
+    const unsigned int sn)
+  {
+    DoFInfo<dim,spacedim>::reinit(c, f, fn, sn);
+    for (unsigned int i=0;i<fevalv.size();++i)
+      {
+       FVB& febase = *fevalv[i];
+       FESubfaceValues<dim>& fe = dynamic_cast<FESubfaceValues<dim>&> (febase);
+       fe.reinit(this->cell, fn, sn);
+      }
+      
+    const bool split_fevalues = this->block_info != 0;
+    fill_local_data(split_fevalues);
+  }
+
+//----------------------------------------------------------------------//
+
+  template <int dim, int sdim>
+  template <typename T>
+  IntegrationInfoBox<dim,sdim>::IntegrationInfoBox(const T& t)
+                 :
+                 cell_info(t),
+                 bdry_info(t),
+                 face_info(t),
+                 subface_info(t),
+                 neighbor_info(t)
+  {}
+    
+  template <int dim, int sdim>
+  template <class WORKER>
+  void
+  IntegrationInfoBox<dim,sdim>::initialize(
+    const WORKER& integrator,
+    const FiniteElement<dim,sdim>& el,
+    const Mapping<dim,sdim>& mapping)
+  {
+    integrator.initialize_info(cell_info, false);
+    integrator.initialize_info(bdry_info, false);
+    integrator.initialize_info(face_info, true);
+    integrator.initialize_info(subface_info, true);
+    integrator.initialize_info(neighbor_info, true);
+      
+    cell_info.initialize<FEValues<dim,sdim> >(el, mapping, integrator.cell_quadrature,
+                        integrator.cell_flags);
+    bdry_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.bdry_quadrature,
+                        integrator.face_flags);
+    face_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
+                        integrator.face_flags);
+    subface_info.initialize<FESubfaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
+                           integrator.face_flags);
+    neighbor_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
+                            integrator.ngbr_flags);
+  }
+
+    
+  template <int dim, int sdim>
+  template <class WORKER, class VECTOR>
+  void
+  IntegrationInfoBox<dim,sdim>::initialize(
+    const WORKER& integrator,
+    const FiniteElement<dim,sdim>& el,
+    const Mapping<dim,sdim>& mapping,
+    const NamedData<SmartPointer<VECTOR> >& data)
+  {
+    cell_info.initialize_data(data);
+    bdry_info.initialize_data(data);
+    face_info.initialize_data(data);
+    subface_info.initialize_data(data);
+    neighbor_info.initialize_data(data);
+
+    initialize(integrator, el, mapping);
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h
new file mode 100644 (file)
index 0000000..c92dbc1
--- /dev/null
@@ -0,0 +1,192 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <numerics/mesh_worker_info.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace MeshWorker
+{
+  template <int dim, int spacedim>
+  DoFInfo<dim,spacedim>::DoFInfo(const BlockInfo& info)
+                 : block_info(&info, typeid(*this).name())
+  {}
+
+
+  template <int dim, int spacedim>
+  void
+  DoFInfo<dim,spacedim>::get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator c)
+  {
+    indices.resize(c->get_fe().dofs_per_cell);
+  
+    if (block_info == 0 || block_info->local().size() == 0)
+      c->get_dof_indices(indices);
+    else
+      {
+       indices_org.resize(c->get_fe().dofs_per_cell);
+       c->get_dof_indices(indices_org);
+       for (unsigned int i=0;i<indices.size();++i)
+         indices[this->block_info->renumber(i)] = indices_org[i];
+      }
+  }
+
+
+  template <int dim, int spacedim>
+  void
+  DoFInfo<dim,spacedim>::get_indices(const typename MGDoFHandler<dim, spacedim>::cell_iterator c)
+  {
+    indices.resize(c->get_fe().dofs_per_cell);
+  
+    if (block_info == 0 || block_info->local().size() == 0)
+      c->get_mg_dof_indices(indices);
+    else
+      {
+       indices_org.resize(c->get_fe().dofs_per_cell);
+       c->get_mg_dof_indices(indices_org);
+       for (unsigned int i=0;i<indices.size();++i)
+         indices[this->block_info->renumber(i)] = indices_org[i];
+      }
+  }
+
+//----------------------------------------------------------------------//
+
+  template<int dim, class FVB, int sdim>
+  IntegrationInfo<dim,FVB,sdim>::IntegrationInfo(const BlockInfo& block_info)
+                 :
+                 DoFInfo<dim,sdim>(block_info),
+                 fevalv(0),
+                 multigrid(false),
+                 global_data(boost::shared_ptr<VectorDataBase<dim, sdim> >(new VectorDataBase<dim, sdim>))
+  {}
+
+
+  template<int dim, class FVB, int sdim>
+  template <class FEVALUES>
+  void
+  IntegrationInfo<dim,FVB,sdim>::initialize(
+    const FiniteElement<dim,sdim>& el,
+    const Mapping<dim,sdim>& mapping,
+    const Quadrature<FEVALUES::integral_dimension>& quadrature,
+    const UpdateFlags flags)
+  {
+    if (this->block_info == 0 || this->block_info->local().size() == 0)
+      {
+       fevalv.resize(1);             
+       fevalv[0] = boost::shared_ptr<FVB> (
+         new FEVALUES (mapping, el, quadrature, flags));
+      }
+    else
+      {
+       fevalv.resize(el.n_base_elements());
+       for (unsigned int i=0;i<fevalv.size();++i)
+         {
+           fevalv[i] = boost::shared_ptr<FVB> (
+             new FEVALUES (mapping, el.base_element(i), quadrature, flags));
+         }
+      }
+
+    values.resize(global_data->n_values());
+                                    // For all selected finite
+                                    // element functions
+    for (unsigned int i=0;i<values.size();++i)
+      {
+       values[i].resize(this->local_indices().size());
+                                        // For all components
+       for (unsigned int j=0;j<values[i].size();++j)
+         {
+           values[i][j].resize(el.n_components());
+         }
+      }
+    
+    gradients.resize(global_data->n_gradients());
+                                    // For all selected finite
+                                    // element functions
+    for (unsigned int i=0;i<gradients.size();++i)
+      {
+       gradients[i].resize(this->local_indices().size());
+                                        // For all components
+       for (unsigned int j=0;j<gradients[i].size();++j)
+         {
+           gradients[i][j].resize(el.n_components());
+         }
+      }
+    
+    hessians.resize(global_data->n_hessians());
+                                    // For all selected finite
+                                    // element functions
+    for (unsigned int i=0;i<hessians.size();++i)
+      {
+       hessians[i].resize(this->local_indices().size());
+                                        // For all components
+       for (unsigned int j=0;j<hessians[i].size();++j)
+         {
+           hessians[i][j].resize(el.n_components());
+         }
+      }
+  }
+  
+
+  template<int dim, class FVB, int sdim>
+  void
+  IntegrationInfo<dim,FVB,sdim>::initialize_data(
+    const boost::shared_ptr<VectorDataBase<dim,sdim> > data)
+  {
+    global_data = data;
+  }
+
+
+  template<int dim, class FVB, int sdim>
+  void
+  IntegrationInfo<dim,FVB,sdim>::clear()
+  {
+    fevalv.resize(0);
+  }
+
+
+
+  template<int dim, class FVB, int sdim>
+  void
+  IntegrationInfo<dim,FVB,sdim>::fill_local_data(bool split_fevalues)
+  {
+     if (split_fevalues)
+       {
+       unsigned int comp = 0;
+                                        // Loop over all blocks
+       for (unsigned int b=0;b<this->block_info->local().size();++b)
+         {
+           const unsigned int fe_no = this->block_info->base_element(b);
+           const FEValuesBase<dim>& fe = this->fe(fe_no);
+           const unsigned int n_comp = fe.get_fe().n_components();
+           const unsigned int block_start = this->block_info->local().block_start(b);
+           const unsigned int block_size = this->block_info->local().block_size(b);
+           
+           this->global_data->fill(values, gradients, hessians, fe, this->indices,
+                                   comp, n_comp, block_start, block_size);
+           comp += n_comp;
+         }
+       }
+     else
+       {
+        const FEValuesBase<dim>& fe = this->fe(0);
+        const unsigned int n_comp = fe.get_fe().n_components();
+        this->global_data->fill(values, gradients, hessians, fe, this->indices,
+                                0, n_comp, 0, this->indices.size());
+       }
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_loop.h b/deal.II/deal.II/include/numerics/mesh_worker_loop.h
new file mode 100644 (file)
index 0000000..c10f499
--- /dev/null
@@ -0,0 +1,188 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
+//
+//    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_loop_h
+#define __deal2__mesh_worker_loop_h
+
+DEAL_II_NAMESPACE_OPEN
+
+template <typename> class TriaActiveIterator;
+
+namespace MeshWorker
+{
+  namespace internal
+  {
+/**
+ * Find out if an iterator supports inactive cells.
+ */
+    template <class DI>
+    inline bool is_active_iterator(const DI&)
+    {
+      return false;
+    }
+
+    template <class ACCESSOR>
+    inline bool is_active_iterator(const TriaActiveIterator<ACCESSOR>&)
+    {
+      return true;
+    }
+  }
+  
+
+/**
+ * The main work function of this namespace. Its action consists of two
+ * loops.
+ *
+ * First, a loop over all cells in the iterator range is performed, in
+ * each step updating the CellInfo object, then calling
+ * LocalWorker::cell() with this object.
+ *
+ * In the second loop, we work through all the faces of cells in the
+ * iterator range. The functions LocalWorker::bdry() and
+ * LocalWorker::face() are called for each boundary and interior face,
+ * respectively. Unilaterally refined interior faces are handled
+ * automatically by the loop.
+ *
+ * The extend of the second loop will be determined by the two control
+ * variables LocalWorker::boundary_fluxes and
+ * LocalWorker::interior_fluxes.
+ */
+  template<class ITERATOR, class ENDITERATOR, class CELLINFO, class FACEINFO, class LOCALWORKER>
+  void loop(ITERATOR begin, ENDITERATOR end,
+           CELLINFO& cellinfo, FACEINFO& bdryinfo,
+           FACEINFO& faceinfo, FACEINFO& subfaceinfo, FACEINFO& ngbrinfo,
+           LOCALWORKER& localworker,
+           bool cells_first = true)
+  {
+    const_cast<Triangulation<ITERATOR::AccessorType::Container::dimension>&>(
+      begin->get_triangulation()).clear_user_flags();
+    
+                                    // Loop over all cells
+    for (ITERATOR cell = begin; cell != end; ++cell)
+      {
+                                        // Execute this, if cells
+                                        // have to be dealt with
+                                        // before faces
+       if (cells_first)
+         {
+           cellinfo.reinit(cell);
+           localworker.cell(cellinfo);
+         }
+       
+       if (localworker.interior_fluxes || localworker.boundary_fluxes)
+         for (unsigned int face_no=0; face_no < GeometryInfo<ITERATOR::AccessorType::Container::dimension>::faces_per_cell; ++face_no)
+           {
+             typename ITERATOR::AccessorType::Container::face_iterator face = cell->face(face_no);
+                                              // Treat every face only once
+             if (face->user_flag_set ()) continue;
+             
+             if (cell->at_boundary(face_no))
+               {
+                 if (localworker.boundary_fluxes)
+                   {
+                     bdryinfo.reinit(cell, face, face_no);
+                     localworker.bdry(bdryinfo);
+                   }
+               }
+             else if (localworker.interior_fluxes)
+               {
+                 if (face->user_flag_set ()) continue;
+                 face->set_user_flag ();
+                                                  // Interior face
+                 typename ITERATOR::AccessorType::Container::cell_iterator
+                   neighbor = cell->neighbor(face_no);
+                 
+                                                  // Deal with
+                                                  // refinement edges
+                                                  // from the refined
+                                                  // side. Assuming
+                                                  // one-irregular
+                                                  // meshes, this
+                                                  // situation should
+                                                  // only occur if
+                                                  // both cells are
+                                                  // active.
+                 if (neighbor->level() < cell->level())
+                   {
+                     Assert(!cell->has_children(), ExcInternalError());
+                     Assert(!neighbor->has_children(), ExcInternalError());
+                     
+                     std::pair<unsigned int, unsigned int> neighbor_face_no
+                       = cell->neighbor_of_coarser_neighbor(face_no);
+                     typename ITERATOR::AccessorType::Container::face_iterator nface
+                       = neighbor->face(neighbor_face_no.first);
+                     faceinfo.reinit(cell, face, face_no);
+                     subfaceinfo.reinit(neighbor, nface, neighbor_face_no.first, neighbor_face_no.second);
+                                                      // Neighbor
+                                                      // first to
+                                                      // conform to
+                                                      // old version
+                     localworker.face(faceinfo, subfaceinfo);
+                   }
+                 else
+                   {
+                                                      // Neighbor is
+                                                      // on same
+                                                      // level
+
+                                                      // If iterator
+                                                      // is active
+                                                      // and neighbor
+                                                      // is refined,
+                                                      // skip
+                                                      // internal face.
+                     if (internal::is_active_iterator(cell) && neighbor->has_children())
+                       continue;
+                     
+                     unsigned int neighbor_face_no = cell->neighbor_of_neighbor(face_no);
+                     Assert (neighbor->face(neighbor_face_no) == face, ExcInternalError());
+                                                      // Regular interior face
+                     faceinfo.reinit(cell, face, face_no);
+                     ngbrinfo.reinit(neighbor, neighbor->face(neighbor_face_no),
+                                     neighbor_face_no);
+                     
+                     localworker.face(faceinfo, ngbrinfo);
+                     neighbor->face(neighbor_face_no)->set_user_flag ();
+                   }
+                 
+               }
+           } // faces
+                                        // Execute this, if faces
+                                        // have to be handled first
+       if (!cells_first)
+         {
+           cellinfo.reinit(cell);
+           localworker.cell(cellinfo);
+         }
+      }
+  }
+
+/**
+ * Simplified interface for loop() if specialized for integration.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template<int dim, class ITERATOR, class ENDITERATOR, class LOCALWORKER>
+  void integration_loop(ITERATOR begin, ENDITERATOR end,
+                       IntegrationInfoBox<dim, dim>& box,
+                       LOCALWORKER& localworker,
+                       bool cells_first = true)
+  {
+    loop(begin, end, box.cell_info, box.bdry_info,
+        box.face_info, box.subface_info, box.neighbor_info,
+        localworker, cells_first);
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
new file mode 100644 (file)
index 0000000..45397d8
--- /dev/null
@@ -0,0 +1,375 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 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_vector_selector_h
+#define __deal2__mesh_worker_vector_selector_h
+
+#include <base/named_data.h>
+#include <base/tensor.h>
+#include <base/smartpointer.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template<int,int> class FEValuesBase;
+
+namespace MeshWorker
+{
+  
+/**
+ * A class that selects vectors from a list of named vectors.
+ *
+ * Since the number of vectors in FEVectors may grow with every
+ * nesting of applications or loops, it is important to be able to
+ * select those, which are actually used in computing residuals etc.
+ * This class organizes the selection.
+ *
+ * It is used for instance in IntegrationWorker to
+ * determine which values, derivatives or second derivatives are
+ * actually computed.
+ *
+ * @author Guido Kanschat 2009
+ */
+  class VectorSelector :
+      public Subscriptor
+  {
+    public:
+                                      /**
+                                       * Add a vector to the
+                                       * selection. The arguments are
+                                       * the name of the vector and
+                                       * indicators, which
+                                       * information is to be
+                                       * extracted from the vector.
+                                       */
+      void add(const std::string& name,
+              bool values = true,
+              bool gradients = false,
+              bool hessians = false);
+
+                                      /**
+                                       * Initialize the selection
+                                       * field with a data vector.
+                                       */
+      template <class DATA>
+      void initialize(const NamedData<DATA>&);
+      
+                                      /**
+                                       * Check whether any vector is selected.
+                                       */
+      bool empty () const;
+
+                                      /**
+                                       * Returns true if values are
+                                       * selected for any vector.
+                                       */
+      bool has_values () const;
+      
+                                      /**
+                                       * Returns true if gradients are
+                                       * selected for any vector.
+                                       */
+      bool has_gradients () const;
+      
+                                      /**
+                                       * Returns true if hessians are
+                                       * selected for any vector.
+                                       */
+      bool has_hessians () const;
+      
+                                      /**
+                                       * Number of vectors for values
+                                       */
+      unsigned int n_values () const;
+      
+                                      /**
+                                       * Number of vectors for gradients
+                                       */
+      unsigned int n_gradients () const;
+      
+                                      /**
+                                       * Number of vectors for Hessians
+                                       */
+      unsigned int n_hessians () const;
+      
+                                      /**
+                                       * The vector index for the ith value
+                                       */
+      unsigned int value_index (unsigned int i) const;
+      
+                                      /**
+                                       * The vector index for the ith gradient
+                                       */
+      unsigned int gradient_index (unsigned int i) const;
+      
+                                      /**
+                                       * The vector index for the ith Hessian
+                                       */
+      unsigned int hessian_index (unsigned int i) const;
+      
+                                      /**
+                                       * Print the contents of the
+                                       * selection to the stream.
+                                       */
+      template <class STREAM, typename DATA>
+      void print (STREAM& s, const NamedData<DATA>& v) const;
+      
+    protected:
+                                      /**
+                                       * Selection of the vectors
+                                       * used to compute values.
+                                       */
+      NamedSelection value_selection;
+      
+                                      /**
+                                       * Selection of the vectors
+                                       * used to compute gradients.
+                                       */
+      NamedSelection gradient_selection;
+      
+                                      /**
+                                       * Selection of the vectors
+                                       * used to compute hessians.
+                                       */
+      NamedSelection hessian_selection;      
+  };
+
+/**
+ * Based on VectorSelector, this is the class used by IntegrationInfo
+ * to compute values of source vectors in quadrature points.
+ *
+ * @guido Kanschat, 2009
+ */
+  template <int dim, int spacedim = dim>
+  class VectorDataBase :
+      public VectorSelector
+  {
+    public:
+                                      /**
+                                       * Virtual, but empty
+                                       * destructor.
+                                       */
+      virtual ~VectorDataBase();
+      
+                                      /**
+                                       * The only function added to
+                                       * VectorSelector is an
+                                       * abstract virtual function
+                                       * implemented in the derived
+                                       * class template and called by
+                                       * IntegrationInfo.
+                                       *
+                                       * Depending on the selections
+                                       * made in our base class, this
+                                       * fills the first three
+                                       * arguments with the local
+                                       * data of the finite element
+                                       * functions. It is usually
+                                       * called either for the whole
+                                       * FESystem, or for each base
+                                       * element separately.
+                                       *
+                                       * @param values is the vector
+                                       * filled with the values of
+                                       * the finite element function
+                                       * in the quadrature points.
+                                       *
+                                       * @param gradients is the vector
+                                       * filled with the derivatives of
+                                       * the finite element function
+                                       * in the quadrature points.
+                                       *
+                                       * @param hessians is the
+                                       * vector filled with the
+                                       * second derivatives of the
+                                       * finite element function in
+                                       * the quadrature points.
+                                       *
+                                       * @param fe is the
+                                       * FEValuesBase object which is
+                                       * used to compute the function
+                                       * values. Its UpdateFlags must
+                                       * have been set appropriately.
+                                       *
+                                       * @param index is the local
+                                       * index vector. If @p fe
+                                       * refers to base elements of
+                                       * the system, this vector
+                                       * should be sorted by block
+                                       * and the arguments @p start
+                                       * and @p size below specify
+                                       * the subset of @p indices
+                                       * used.
+                                       *
+                                       * @param block is the block
+                                       * number processed, or zero if
+                                       * no base elements are used.
+                                       *
+                                       * @param start is the first
+                                       * index of this block in @p
+                                       * indices, or zero if
+                                       * no base elements are used.
+                                       *
+                                       * @param size is the number of
+                                       * dofs per cell of the current
+                                       * element or base element.
+                                       */
+      virtual void fill(
+       std::vector<std::vector<std::vector<double> > >& values,
+       std::vector<std::vector<std::vector<Tensor<1,dim> > > >& gradients,
+       std::vector<std::vector<std::vector<Tensor<2,dim> > > >& hessians,
+       const FEValuesBase<dim,spacedim>& fe,
+       const std::vector<unsigned int>& index,
+       unsigned int component,
+       unsigned int n_comp,
+       unsigned int start,
+       unsigned int size) const;
+  };
+
+/**
+ * Based on VectorSelector, this is the class that implements the
+ * function VectorDataBase::fill() for a certain type of vector.
+ *
+ * @author Guido Kanschat, 2009
+ */
+  template <class VECTOR, int dim, int spacedim = dim>
+  class VectorData :
+      public VectorDataBase<dim, spacedim>
+  {
+    public:
+      void initialize(const NamedData<SmartPointer<VECTOR> >&);
+      
+      virtual void fill(
+       std::vector<std::vector<std::vector<double> > >& values,
+       std::vector<std::vector<std::vector<Tensor<1,dim> > > >& gradients,
+       std::vector<std::vector<std::vector<Tensor<2,dim> > > >& hessians,
+       const FEValuesBase<dim,spacedim>& fe,
+       const std::vector<unsigned int>& index,
+       unsigned int component,
+       unsigned int n_comp,
+       unsigned int start,
+       unsigned int size) const;
+    private:
+      SmartPointer<const NamedData<SmartPointer<VECTOR> > > data;
+  };
+  
+//----------------------------------------------------------------------//
+
+  inline void
+  VectorSelector::add(const std::string& name, bool values, bool gradients, bool hessians)
+  {
+    if (values) value_selection.add(name);
+    if (gradients) gradient_selection.add(name);
+    if (hessians) hessian_selection.add(name);  
+  }
+
+
+  template <typename DATA>
+  inline void
+  VectorSelector::initialize(const NamedData<DATA>& src)
+  {
+    value_selection.initialize(src);
+    gradient_selection.initialize(src);
+    hessian_selection.initialize(src);
+  }
+
+  inline bool
+  VectorSelector::empty() const
+  {
+    return (value_selection.size() == 0 &&
+           gradient_selection.size() == 0 &&
+           hessian_selection.size() == 0);
+  }
+  
+  
+  inline bool
+  VectorSelector::has_values() const
+  {
+    return value_selection.size() != 0;
+  }
+  
+
+  inline bool
+  VectorSelector::has_gradients() const
+  {
+    return gradient_selection.size() != 0;
+  }
+  
+
+  inline bool
+  VectorSelector::has_hessians() const
+  {
+    return hessian_selection.size() != 0;
+  }
+  
+
+  inline unsigned int
+  VectorSelector::n_values() const
+  {
+    return value_selection.size();
+  }
+  
+
+  inline unsigned int
+  VectorSelector::n_gradients() const
+  {
+    return gradient_selection.size();
+  }
+  
+
+  inline unsigned int
+  VectorSelector::n_hessians() const
+  {
+    return hessian_selection.size();
+  }
+  
+
+  inline unsigned int
+  VectorSelector::value_index(unsigned int i) const
+  {
+    return value_selection(i);
+  }
+  
+
+  inline unsigned int
+  VectorSelector::gradient_index(unsigned int i) const
+  {
+    return gradient_selection(i);
+  }
+  
+
+  inline unsigned int
+  VectorSelector::hessian_index(unsigned int i) const
+  {
+    return hessian_selection(i);
+  }
+  
+
+  template <class STREAM, typename DATA>
+  inline void
+  VectorSelector::print(STREAM& s, const NamedData<DATA>& v) const
+  {
+    s << "values:   ";
+    for (unsigned int i=0;i<n_values();++i)
+      s << " '" << v.name(value_selection(i)) << '\'';
+    s << std::endl << "gradients:";
+    for (unsigned int i=0;i<n_gradients();++i)
+      s << " '" << v.name(gradient_selection(i)) << '\'';
+    s << std::endl << "hessians: ";
+    for (unsigned int i=0;i<n_hessians();++i)
+      s << " '" << v.name(hessian_selection(i)) << '\'';
+    s << std::endl;
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.inst.in b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.inst.in
new file mode 100644 (file)
index 0000000..17ba91a
--- /dev/null
@@ -0,0 +1,17 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_field_function.inst.in 19046 2009-07-08 19:30:23Z bangerth $
+//
+//    Copyright (C) 2009 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.
+//
+//---------------------------------------------------------------------------
+
+
+for (VECTOR : SERIAL_VECTORS)
+{
+  template class VectorData<VECTOR,deal_II_dimension>;
+}
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h
new file mode 100644 (file)
index 0000000..b6cf863
--- /dev/null
@@ -0,0 +1,86 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/mesh_worker_vector_selector.h>
+#include <base/vector_slice.h>
+#include <fe/fe_values.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+  template <int dim, int spacedim>
+  VectorDataBase<dim, spacedim>::~VectorDataBase()
+  {}
+  
+  template <int dim, int spacedim>
+  void
+  VectorDataBase<dim, spacedim>::fill(
+       std::vector<std::vector<std::vector<double> > >&,
+       std::vector<std::vector<std::vector<Tensor<1,dim> > > >&,
+       std::vector<std::vector<std::vector<Tensor<2,dim> > > >&,
+       const FEValuesBase<dim,spacedim>&,
+       const std::vector<unsigned int>&,
+       unsigned int,
+       unsigned int,
+       unsigned int,
+       unsigned int) const
+  {}
+//----------------------------------------------------------------------//
+
+  template <class VECTOR, int dim, int spacedim>
+  void
+  VectorData<VECTOR, dim, spacedim>::initialize(const NamedData<SmartPointer<VECTOR> >& d)
+  {
+    data = &d;
+    VectorSelector::initialize(d);
+  }
+  
+  
+  template <class VECTOR, int dim, int spacedim>
+  void
+  VectorData<VECTOR, dim, spacedim>::fill(
+       std::vector<std::vector<std::vector<double> > >& values,
+       std::vector<std::vector<std::vector<Tensor<1,dim> > > >& gradients,
+       std::vector<std::vector<std::vector<Tensor<2,dim> > > >& hessians,
+       const FEValuesBase<dim,spacedim>& fe,
+       const std::vector<unsigned int>& index,
+       unsigned int component,
+       unsigned int n_comp,
+       unsigned int start,
+       unsigned int size) const
+  {
+    for (unsigned int i=0;i<this->n_values();++i)
+      {
+       const VECTOR& src = *(*data)(this->value_index(i));
+       VectorSlice<std::vector<std::vector<double> > > dst(values[i], component, n_comp);
+       fe.get_function_values(src, make_slice(index, start, size), dst, true);
+      }
+    
+    for (unsigned int i=0;i<this->n_gradients();++i)
+      {
+       const VECTOR& src = *(*data)(this->value_index(i));
+       VectorSlice<std::vector<std::vector<Tensor<1,dim> > > > dst(gradients[i], component, n_comp);
+       fe.get_function_gradients(src, make_slice(index, start, size), dst, true);
+      }
+    
+    for (unsigned int i=0;i<this->n_hessians();++i)
+      {
+       const VECTOR& src = *(*data)(this->value_index(i));
+       VectorSlice<std::vector<std::vector<Tensor<2,dim> > > > dst(hessians[i], component, n_comp);
+       fe.get_function_hessians(src, make_slice(index, start, size), dst, true);
+      }
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
index 4082dd26a45f3c505e10bc4dc66c739dfa059715..54a65c46a79b1fb94f603b38cfb7361b32de93c7 100644 (file)
@@ -1,8 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
-//    Version: $Name$
 //
-//    Copyright (C) 2007, 2008 by the deal.II authors
+//    Copyright (C) 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
diff --git a/deal.II/deal.II/source/numerics/mesh_worker.cc b/deal.II/deal.II/source/numerics/mesh_worker.cc
new file mode 100644 (file)
index 0000000..465eabd
--- /dev/null
@@ -0,0 +1,116 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <multigrid/mg_tools.h>
+#include <fe/fe.h>
+#include <fe/fe_tools.h>
+#include <numerics/mesh_worker.h>
+#include <base/quadrature_lib.h>
+
+using namespace dealii;
+using namespace MeshWorker;
+
+template <int dim>
+IntegrationWorker<dim>::IntegrationWorker ()
+{
+  cell_flags = update_JxW_values;
+  bdry_flags = UpdateFlags(update_JxW_values | update_normal_vectors);
+  face_flags = bdry_flags;
+  ngbr_flags = update_default;
+}
+
+
+template<int dim>
+void
+IntegrationWorker<dim>::initialize_selectors(
+  const VectorSelector& cs,
+  const VectorSelector& bs,
+  const VectorSelector& fs)
+{
+  cell_selector = cs;
+  bdry_selector = bs;
+  face_selector = fs;
+  
+  if (cell_selector.has_values() != 0) cell_flags |= update_values;
+  if (cell_selector.has_gradients() != 0) cell_flags |= update_gradients;
+  if (cell_selector.has_hessians() != 0) cell_flags |= update_hessians;
+  
+  if (bdry_selector.has_values() != 0) bdry_flags |= update_values;
+  if (bdry_selector.has_gradients() != 0) bdry_flags |= update_gradients;
+  if (bdry_selector.has_hessians() != 0) bdry_flags |= update_hessians;
+  
+  if (face_selector.has_values() != 0) face_flags |= update_values;
+  if (face_selector.has_gradients() != 0) face_flags |= update_gradients;
+  if (face_selector.has_hessians() != 0) face_flags |= update_hessians;
+  
+  if (face_selector.has_values() != 0) ngbr_flags |= update_values;
+  if (face_selector.has_gradients() != 0) ngbr_flags |= update_gradients;
+  if (face_selector.has_hessians() != 0) ngbr_flags |= update_hessians;  
+}
+
+
+template<int dim>
+void
+IntegrationWorker<dim>::add_selector(
+  const std::string& name, bool values, bool gradients, bool hessians,
+  bool cell, bool bdry, bool face)
+{
+  if (cell) cell_selector.add(name, values, gradients, hessians);
+  if (bdry) bdry_selector.add(name, values, gradients, hessians);
+  if (face) face_selector.add(name, values, gradients, hessians);  
+  
+  if (cell_selector.has_values() != 0) cell_flags |= update_values;
+  if (cell_selector.has_gradients() != 0) cell_flags |= update_gradients;
+  if (cell_selector.has_hessians() != 0) cell_flags |= update_hessians;
+
+  if (bdry_selector.has_values() != 0) bdry_flags |= update_values;
+  if (bdry_selector.has_gradients() != 0) bdry_flags |= update_gradients;
+  if (bdry_selector.has_hessians() != 0) bdry_flags |= update_hessians;
+  
+  if (face_selector.has_values() != 0) face_flags |= update_values;
+  if (face_selector.has_gradients() != 0) face_flags |= update_gradients;
+  if (face_selector.has_hessians() != 0) face_flags |= update_hessians;
+  
+  if (face_selector.has_values() != 0) ngbr_flags |= update_values;
+  if (face_selector.has_gradients() != 0) ngbr_flags |= update_gradients;
+  if (face_selector.has_hessians() != 0) ngbr_flags |= update_hessians;  
+}
+
+
+template<int dim>
+void
+IntegrationWorker<dim>::add_update_flags(
+  const UpdateFlags flags, bool cell, bool bdry, bool face, bool ngbr)
+{
+  if (cell) cell_flags |= flags;
+  if (bdry) bdry_flags |= flags;
+  if (face) face_flags |= flags;
+  if (ngbr) ngbr_flags |= flags;  
+}
+
+  
+template<int dim>
+void
+IntegrationWorker<dim>::initialize_gauss_quadrature(
+  unsigned int cp,
+  unsigned int bp,
+  unsigned int fp)
+{
+  cell_quadrature = QGauss<dim>(cp);
+  bdry_quadrature = QGauss<dim-1>(bp);
+  face_quadrature = QGauss<dim-1>(fp);
+
+}
+
+  
+template class IntegrationWorker<deal_II_dimension>;
+
diff --git a/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc b/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc
new file mode 100644 (file)
index 0000000..c575931
--- /dev/null
@@ -0,0 +1,130 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <numerics/mesh_worker_assembler.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace MeshWorker
+{
+  namespace Assembler
+  {
+//----------------------------------------------------------------------//
+    
+    template <typename number>
+    inline void
+    Functional<number>::initialize(unsigned int n)
+    {
+      results.resize(n);
+    }
+    
+    
+    template <typename number>
+    template<int dim>
+    inline void
+    Functional<number>::assemble(const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       results[i] += info.J[i];
+    }
+    
+
+    template <typename number>
+    template<int dim>
+    inline void
+    Functional<number>::assemble(const DoFInfo<dim>& info1,
+                                const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       {
+         results[i] += info1.J[i];
+         results[i] += info2.J[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", DataVectors::ExcNameMismatch(0, "cells"));
+      if (sep)
+       Assert(r.name(1) == "faces", DataVectors::ExcNameMismatch(1, "faces"));
+      AssertDimension(r(0).n_blocks(), r(1).n_blocks());
+      
+      results = r;
+      separate_faces = sep;
+    }
+    
+    
+    template <typename number>
+    template<int dim>
+    inline void
+    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<info.J.size();++i)
+       {
+         if (separate_faces &&
+             info.face_number != deal_II_numbers::invalid_unsigned_int)
+           results.vector(1).block(i)(info.face->user_index()) += info.J[i];
+         else
+           results.vector(0).block(i)(info.cell->user_index()) += info.J[i];
+       }
+    }
+    
+
+    template <typename number>
+    template<int dim>
+    inline void
+    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
+                                   const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<info1.J.size();++i)
+       {
+         if (separate_faces)
+           {
+             const double J = info1.J[i] + info2.J[i];
+             results.vector(1).block(i)(info1.face->user_index()) += J;
+             if (info2.face != info1.face)
+               results.vector(1).block(i)(info2.face->user_index()) += J;
+           }
+         else
+           {
+             results.vector(0).block(i)(info1.cell->user_index()) += info1.J[i];
+             results.vector(0).block(i)(info2.cell->user_index()) += info2.J[i];
+           }
+       }
+    }
+    
+
+//----------------------------------------------------------------------//
+
+
+
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/mesh_worker_info.cc b/deal.II/deal.II/source/numerics/mesh_worker_info.cc
new file mode 100644 (file)
index 0000000..d7c6e77
--- /dev/null
@@ -0,0 +1,42 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <numerics/mesh_worker_info.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+  template class LocalResults<double>;
+  template class DoFInfo<deal_II_dimension,deal_II_dimension>;
+  
+#include "mesh_worker_info.inst"
+  
+  template void IntegrationInfo<deal_II_dimension, FEValuesBase<deal_II_dimension> >
+  ::initialize<FEValues<deal_II_dimension> >(
+    const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
+    const Quadrature<FEValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+  template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
+  ::initialize<FEFaceValues<deal_II_dimension> >(
+    const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
+    const Quadrature<FEFaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+  template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
+  ::initialize<FESubfaceValues<deal_II_dimension> >(
+    const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
+    const Quadrature<FESubfaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
diff --git a/deal.II/deal.II/source/numerics/mesh_worker_info.inst.in b/deal.II/deal.II/source/numerics/mesh_worker_info.inst.in
new file mode 100644 (file)
index 0000000..8e123d5
--- /dev/null
@@ -0,0 +1,34 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_field_function.inst.in 19046 2009-07-08 19:30:23Z bangerth $
+//
+//    Copyright (C) 2009 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.
+//
+//---------------------------------------------------------------------------
+
+
+for (FEV : FEVALUES_BASES)
+{
+  template class IntegrationInfo<deal_II_dimension, FEV>;
+}
+
+for (FEV : FEVALUES_BASES)
+{
+//   for (TYPE : DERIVATIVE_TENSORS)
+//   {
+//     template void IntegrationInfo<deal_II_dimension, FEV>::fill_local_data(
+//       std::vector<std::vector<std::vector<TYPE > > >&, bool) const;
+//     template void IntegrationInfo<deal_II_dimension, FEV>::fill_local_data(
+//       std::vector<std::vector<std::vector<TYPE > > >&, VectorSelector&, bool) const;
+//   }
+//     for (VECTOR : SERIAL_VECTORS)
+//       {
+//       }
+ }
+
+
+
diff --git a/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc
new file mode 100644 (file)
index 0000000..9296026
--- /dev/null
@@ -0,0 +1,28 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+#include <numerics/mesh_worker_vector_selector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MeshWorker
+{
+  template class VectorDataBase<deal_II_dimension>;
+
+#include "mesh_worker_vector_selector.inst"
+}
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in
new file mode 100644 (file)
index 0000000..17ba91a
--- /dev/null
@@ -0,0 +1,17 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_field_function.inst.in 19046 2009-07-08 19:30:23Z bangerth $
+//
+//    Copyright (C) 2009 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.
+//
+//---------------------------------------------------------------------------
+
+
+for (VECTOR : SERIAL_VECTORS)
+{
+  template class VectorData<VECTOR,deal_II_dimension>;
+}
index 2e79616d83cdb26ca1d504acab0bc0dde19a6432..2feb006baf95ecc07e9cea7c796c8c645f7bd359 100644 (file)
@@ -38,7 +38,8 @@
 
                                 // Here come the new include files
                                 // for using the MeshWorker framework
-
+#include <numerics/mesh_worker.h>
+#include <numerics/mesh_worker_loop.h>
 
 #include <iostream>
 #include <fstream>
@@ -178,280 +179,219 @@ void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
 }
 
 
-                                // @sect3{Class: DGTransportEquation}
+                                // @sect3{Integrating cell and face matrices}
+                                //
+                                // In order to use the MeshWorker
+                                // framework, we need to define a
+                                // worker object, which serves three
+                                // purposes. First, it performs the
+                                // local integration on each mesh
+                                // cell. Second, it controls the loop
+                                // over cells and faces and the data
+                                // being generated. This part is
+                                // easily implemented here, since the
+                                // default structure can be
+                                // used. Finally, through its base
+                                // class from the
+                                // MeshWorker::Assembler namespace,
+                                // this object determins, how locally
+                                // produced data is assembled into
+                                // the global system.
                                 //
-                                // Next we define the
-                                // equation-dependent and
-                                // DG-method-dependent class
-                                // <code>DGTransportEquation</code>. Its
-                                // member functions were already
-                                // mentioned in the Introduction and
-                                // will be explained
-                                // below. Furthermore it includes
-                                // objects of the previously defined
-                                // <code>Beta</code>, <code>RHS</code> and
-                                // <code>BoundaryValues</code> function
-                                // classes.
+                                // The fist base class is
+                                // MeshWorker::IntegrationWorker. It
+                                // has all the control structures the
+                                // generic loop needs. Fortunately,
+                                // most of them havereasonable
+                                // default values. The others will be
+                                // set below.
+                                //
+                                // The second base class is
+                                // MeshWorker::Assembler::SystemSimple,
+                                // which will assemble the matrices
+                                // and residuals on cells and faces,
+                                // respectively, into the global
+                                // system. The "simple" indicates
+                                // that we do not want to use
+                                // sophisticated block structures.
 template <int dim>
-class DGTransportEquation
+class DGIntegrator : public Subscriptor
 {
   public:
-    DGTransportEquation();
-
-    void assemble_cell_term(const FEValues<dim>& fe_v,
-                           FullMatrix<double> &ui_vi_matrix,
-                           Vector<double> &cell_vector) const;
-    
-    void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
-                               FullMatrix<double> &ui_vi_matrix,
-                               Vector<double> &cell_vector) const;
-    
-    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
-                            FullMatrix<double> &ui_vi_matrix,
-                            FullMatrix<double> &ue_vi_matrix) const;
-
-    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
-                            FullMatrix<double> &ui_vi_matrix,
-                            FullMatrix<double> &ue_vi_matrix,
-                            FullMatrix<double> &ui_ve_matrix,
-                            FullMatrix<double> &ue_ve_matrix) const;
+                                    // First, we define the types of
+                                    // the two info objects handed to
+                                    // the local integration
+                                    // functions in order to make our
+                                    // life easier below.
+    typedef typename MeshWorker::IntegrationWorker<dim>::CellInfo CellInfo;
+    typedef typename MeshWorker::IntegrationWorker<dim>::FaceInfo FaceInfo;
+
+                                    // The following three functions
+                                    // are the ones that get called
+                                    // by the generic loop over all
+                                    // ccells and faces. They are the
+                                    // ones doing the actual
+                                    // integration.
+                                    //
+                                    // Note that the arguments are
+                                    // not constant, because the data
+                                    // spae for local matrices and
+                                    // vectors is part of the
+                                    // MeshWorker::IntegrationInfo
+                                    // object (namely, in its base
+                                    // class
+                                    // MeshWorker::LocalResults).
+    void cell(CellInfo& info) const;
+    void bdry(FaceInfo& info) const;
+    void face(FaceInfo& info1, FaceInfo& info2) const;
+
+                                    // Additionally, like in step 12,
+                                    // we have objects of the
+                                    // functions used in this class.
   private:
     const Beta<dim> beta_function;
     const RHS<dim> rhs_function;
     const BoundaryValues<dim> boundary_function;
 };
 
+                                // @sect4{The local integrators}
+                                // These functions are analogous to
+                                // step 12 and differ only in the
+                                // data structures. Instead of
+                                // providing the local matrices
+                                // explicitly in the argument list,
+                                // they are part of the info object.
+
+                                // Note that here we still have the
+                                // local integration loop inside the
+                                // following functions. The program
+                                // would be even shorter, if we used
+                                // premade operators from the
+                                // Operators namespace (which will be
+                                // added soon).
 
 template <int dim>
-DGTransportEquation<dim>::DGTransportEquation ()
-               :
-               beta_function (),
-               rhs_function (),
-               boundary_function ()
-{}
-
-
-                                // @sect4{Function: assemble_cell_term}
-                                //
-                                // The <code>assemble_cell_term</code>
-                                // function assembles the cell terms
-                                // of the discretization.
-                                // <code>ui_vi_matrix</code> is a cell matrix,
-                                // i.e. for a DG method of degree 1,
-                                // it is of size 4 times 4, and
-                                // <code>cell_vector</code> is of size 4.
-                                // When this function is invoked,
-                                // <code>fe_v</code> is already reinit'ed with the
-                                // current cell before and includes
-                                // all shape values needed.
-template <int dim>
-void DGTransportEquation<dim>::assemble_cell_term(
-  const FEValues<dim> &fe_v,
-  FullMatrix<double> &ui_vi_matrix,
-  Vector<double> &cell_vector) const
+void DGIntegrator<dim>::cell(CellInfo& info) const
 {
-                                  // First we ask <code>fe_v</code> for the
-                                  // quadrature weights,
+                                  // First, let us retrieve some of
+                                  // the objects used here from
+                                  // @p info. Note that these objects
+                                  // can handle much more complex
+                                  // structures, thus the access here
+                                  // looks more complicated than
+                                  // might seem necessary.
+  const FEValuesBase<dim>& fe_v = info.fe();
+  FullMatrix<double>& local_matrix = info.M1[0].matrix;
+  Vector<double>& local_vector = info.R[0].block(0);
+
+                                  // With these objects, we continue
+                                  // local integration like in step 12.
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
-
-                                  // Then the flow field beta and the
-                                  // <code>rhs_function</code> are evaluated at
-                                  // the quadrature points,
+  
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> rhs (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
   
-                                  // and the cell matrix and cell
-                                  // vector are assembled due to the
-                                  // terms $-(u,\beta\cdot\nabla
-                                  // v)_K$ and $(f,v)_K$.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
       {
        for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-         ui_vi_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
+         local_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
                              fe_v.shape_value(j,point) *
                              JxW[point];
        
-       cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
+       local_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
       }
 }
 
-
-                                // @sect4{Function: assemble_boundary_term}
-                                //
-                                // The <code>assemble_boundary_term</code>
-                                // function assembles the face terms
-                                // at boundary faces.  When this
-                                // function is invoked, <code>fe_v</code> is
-                                // already reinit'ed with the current
-                                // cell and current face. Hence it
-                                // provides the shape values on that
-                                // boundary face.
+                                // Now the same for the boundary
+                                // terms. Note that now we use
+                                // FEFaceValuesBase in order to get
+                                // normal vectors.
 template <int dim>
-void DGTransportEquation<dim>::assemble_boundary_term(
-  const FEFaceValues<dim>& fe_v,    
-  FullMatrix<double> &ui_vi_matrix,
-  Vector<double> &cell_vector) const
+void DGIntegrator<dim>::bdry(FaceInfo& info) const
 {
-                                  // Again, as in the previous
-                                  // function, we ask the
-                                  // <code>FEValues</code> object for the
-                                  // quadrature weights
+  const FEFaceValuesBase<dim>& fe_v = info.fe();
+  FullMatrix<double>& local_matrix = info.M1[0].matrix;
+  Vector<double>& local_vector = info.R[0].block(0);
+
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
-                                  // but here also for the normals.
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
-                                  // We evaluate the flow field
-                                  // and the boundary values at the
-                                  // quadrature points.
+  
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> g(fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   boundary_function.value_list (fe_v.get_quadrature_points(), g);
 
-                                  // Then we assemble cell vector and
-                                  // cell matrix according to the DG
-                                  // method given in the
-                                  // introduction.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];      
-                                        // We assemble the term
-                                        // $(\beta\cdot n
-                                        // u,v)_{\partial\kappa_+}$,
       if (beta_n>0)
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
          for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-           ui_vi_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
-                              JxW[point];
-      else
-                                        // and the term $(\beta\cdot
-                                        // n g,v)_{\partial
-                                        // \kappa_-\cap\partial\Omega}$,
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         cell_vector(i) -= beta_n *
-                           g[point] *
-                           fe_v.shape_value(i,point) *
-                           JxW[point];
-    }
-}
-
-
-                                // @sect4{Function: assemble_face_term1}
-                                //
-                                // The <code>assemble_face_term1</code>
-                                // function assembles the face terms
-                                // corresponding to the first version
-                                // of the DG method, cf. above. For
-                                // that case, the face terms are
-                                // given as a sum of integrals over
-                                // all cell boundaries.
-                                //
-                                // When this function is invoked,
-                                // <code>fe_v</code> and <code>fe_v_neighbor</code> are
-                                // already reinit'ed with the current
-                                // cell and the neighoring cell,
-                                // respectively, as well as with the
-                                // current face. Hence they provide
-                                // the inner and outer shape values
-                                // on the face.
-                                //
-                                // In addition to the cell matrix
-                                // <code>ui_vi_matrix</code> this function
-                                // gets a new argument
-                                // <code>ue_vi_matrix</code>, that stores
-                                // contributions to the system matrix
-                                // that are based on exterior values
-                                // of $u$ and interior values of
-                                // $v$. Here we note that <code>ue</code> is
-                                // the short notation for <code>u
-                                // exterior</code> and represents $u_h^-$,
-                                // see the introduction.
-template <int dim>
-void DGTransportEquation<dim>::assemble_face_term1(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,      
-  FullMatrix<double> &ui_vi_matrix,
-  FullMatrix<double> &ue_vi_matrix) const
-{
-                                  // Again, as in the previous
-                                  // function, we ask the FEValues
-                                  // objects for the quadrature
-                                  // weights and the normals
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
-                                  // and we evaluate the flow field
-                                  // at the quadrature points.
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-
-                                  // Then we assemble the cell
-                                  // matrices according to the DG
-                                  // method given in the
-                                  // introduction.
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    {
-      const double beta_n=beta[point] * normals[point];
-                                        // We assemble the term
-                                        // $(\beta\cdot n
-                                        // u,v)_{\partial\kappa_+}$,
-      if (beta_n>0)
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-           ui_vi_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
-                              JxW[point];
+           local_matrix(i,j) += beta_n *
+                                fe_v.shape_value(j,point) *
+                                fe_v.shape_value(i,point) *
+                                JxW[point];
       else
-                                        // and the
-                                        // term $(\beta\cdot n
-                                        // \hat u,v)_{\partial
-                                        // \kappa_-}$.
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           ue_vi_matrix(i,k) += beta_n *
-                               fe_v_neighbor.shape_value(k,point) *
-                               fe_v.shape_value(i,point) *
-                               JxW[point];
+         local_vector(i) -= beta_n *
+                            g[point] *
+                            fe_v.shape_value(i,point) *
+                            JxW[point];
     }
 }
 
-
-                                // @sect4{Function: assemble_face_term2}
-                                //
-                                // Now we look at the
-                                // <code>assemble_face_term2</code> function
-                                // that assembles the face terms
-                                // corresponding to the second
-                                // version of the DG method,
-                                // cf. above. For that case the face
-                                // terms are given as a sum of
-                                // integrals over all faces.  Here we
-                                // need two additional cell matrices
-                                // <code>ui_ve_matrix</code> and
-                                // <code>ue_ve_matrix</code> that will store
-                                // contributions due to terms
-                                // involving ui and ve as well as ue
-                                // and ve.
+                                // Finally, the interior face
+                                // terms. The difference here is that
+                                // we receive two info objects, one
+                                // for each cell adjacent to the face
+                                // and we assemble four matrices, one
+                                // for each cell and two for coupling
+                                // back and forth.
 template <int dim>
-void DGTransportEquation<dim>::assemble_face_term2(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,      
-  FullMatrix<double> &ui_vi_matrix,
-  FullMatrix<double> &ue_vi_matrix,
-  FullMatrix<double> &ui_ve_matrix,
-  FullMatrix<double> &ue_ve_matrix) const
+void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
 {
-                                  // the first few lines are the same
+                                  // For quadrature points, weights,
+                                  // etc., we use the
+                                  // FEFaceValuesBase object of the
+                                  // first argument.
+  const FEFaceValuesBase<dim>& fe_v = info1.fe();
+
+                                  // For additional shape functions,
+                                  // we have to ask the neighbors
+                                  // FEFaceValuseBase.
+  const FEFaceValuesBase<dim>& fe_v_neighbor = info2.fe();
+  
+                                  // Then we get references to the
+                                  // four local matrices. The letters
+                                  // u and v refer to trial and test
+                                  // functions, respectively. The
+                                  // numbers indicate the cells
+                                  // provided by info1 and info2. By
+                                  // convention, the two matrices in
+                                  // each info object refer to the
+                                  // test functions on the respective
+                                  // cell. The first matrix contains the
+                                  // interior couplings of that cell,
+                                  // while the second contains the
+                                  // couplings between cells.
+  FullMatrix<double>& u1_v1_matrix = info1.M1[0].matrix;
+  FullMatrix<double>& u2_v1_matrix = info1.M2[0].matrix;
+  FullMatrix<double>& u1_v2_matrix = info2.M2[0].matrix;
+  FullMatrix<double>& u2_v2_matrix = info2.M1[0].matrix;
+  
+                                  // Here, following the previous
+                                  // functions, we would have the
+                                  // local right hand side
+                                  // vectors. Fortunately, the
+                                  // interface terms only involve the
+                                  // solution and the right hand side
+                                  // does not obtain a contribution.
+  
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
@@ -468,7 +408,7 @@ void DGTransportEquation<dim>::assemble_face_term2(
                                           // seen.
          for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
            for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-             ui_vi_matrix(i,j) += beta_n *
+             u1_v1_matrix(i,j) += beta_n *
                                 fe_v.shape_value(j,point) *
                                 fe_v.shape_value(i,point) *
                                 JxW[point];
@@ -479,7 +419,7 @@ void DGTransportEquation<dim>::assemble_face_term2(
                                           // \kappa_+}$,
          for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
            for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-             ui_ve_matrix(k,j) -= beta_n *
+             u1_v2_matrix(k,j) -= beta_n *
                                  fe_v.shape_value(j,point) *
                                  fe_v_neighbor.shape_value(k,point) *
                                  JxW[point];
@@ -490,7 +430,7 @@ void DGTransportEquation<dim>::assemble_face_term2(
                                           // seen, too.
          for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
            for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
-             ue_vi_matrix(i,l) += beta_n *
+             u2_v1_matrix(i,l) += beta_n *
                                  fe_v_neighbor.shape_value(l,point) *
                                  fe_v.shape_value(i,point) *
                                  JxW[point];
@@ -501,7 +441,7 @@ void DGTransportEquation<dim>::assemble_face_term2(
                                           // \kappa_-}$.
          for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
            for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
-             ue_ve_matrix(k,l) -= beta_n *
+             u2_v2_matrix(k,l) -= beta_n *
                                   fe_v_neighbor.shape_value(l,point) *
                                   fe_v_neighbor.shape_value(k,point) *
                                   JxW[point];
@@ -532,8 +472,7 @@ class DGMethod
     
   private:
     void setup_system ();
-    void assemble_system1 ();
-    void assemble_system2 ();
+    void assemble_system ();
     void solve (Vector<double> &solution);
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
@@ -556,13 +495,6 @@ class DGMethod
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
-                                    // We define the quadrature
-                                    // formulae for the cell and the
-                                    // face terms of the
-                                    // discretization.
-    const QGauss<dim>   quadrature;
-    const QGauss<dim-1> face_quadrature;
-    
                                     // And there are two solution
                                     // vectors, that store the
                                     // solutions to the problems
@@ -570,30 +502,19 @@ class DGMethod
                                     // different assembling routines
                                     // <code>assemble_system1</code> and
                                     // <code>assemble_system2</code>;
-    Vector<double>       solution1;
-    Vector<double>       solution2;
-    Vector<double>       right_hand_side;
-    
-                                    // Finally this class includes an
-                                    // object of the
-                                    // DGTransportEquations class
-                                    // described above.
-    const DGTransportEquation<dim> dg;
+    Vector<double>       solution;
+    Vector<double>       right_hand_side;    
 };
 
 
 template <int dim>
 DGMethod<dim>::DGMethod ()
                :
-               mapping (),
                                                 // Change here for DG
                                                 // methods of
                                                 // different degrees.
                 fe (1),
-               dof_handler (triangulation),
-               quadrature (4),
-               face_quadrature (4),
-               dg ()
+               dof_handler (triangulation)
 {}
 
 
@@ -635,723 +556,97 @@ void DGMethod<dim>::setup_system ()
   
   system_matrix.reinit (sparsity_pattern);
 
-  solution1.reinit (dof_handler.n_dofs());
-  solution2.reinit (dof_handler.n_dofs());
+  solution.reinit (dof_handler.n_dofs());
   right_hand_side.reinit (dof_handler.n_dofs());
 }
 
-
-                                // @sect4{Function: assemble_system1}
-                                //
-                                // We proceed with the
-                                // <code>assemble_system1</code> function that
-                                // implements the DG discretization
-                                // in its first version. This
-                                // function repeatedly calls the
-                                // <code>assemble_cell_term</code>,
-                                // <code>assemble_boundary_term</code> and
-                                // <code>assemble_face_term1</code> functions
-                                // of the <code>DGTransportEquation</code>
-                                // object.  The
-                                // <code>assemble_boundary_term</code> covers
-                                // the first case mentioned in the
-                                // introduction.
-                                //
-                                // 1. face is at boundary
-                                //
-                                // This function takes a
-                                // <code>FEFaceValues</code> object as
-                                // argument.  In contrast to that
-                                // <code>assemble_face_term1</code>
-                                // takes two <code>FEFaceValuesBase</code>
-                                // objects; one for the shape
-                                // functions on the current cell and
-                                // the other for shape functions on
-                                // the neighboring cell under
-                                // consideration. Both objects are
-                                // either of class <code>FEFaceValues</code>
-                                // or of class <code>FESubfaceValues</code>
-                                // (both derived from
-                                // <code>FEFaceValuesBase</code>) according to
-                                // the remaining cases mentioned
-                                // in the introduction:
-                                //
-                                // 2. neighboring cell is finer
-                                // (current cell: <code>FESubfaceValues</code>,
-                                // neighboring cell: <code>FEFaceValues</code>);
-                                //
-                                // 3. neighboring cell is of the same
-                                // refinement level (both, current
-                                // and neighboring cell:
-                                // <code>FEFaceValues</code>);
-                                //
-                                // 4. neighboring cell is coarser
-                                // (current cell: <code>FEFaceValues</code>,
-                                // neighboring cell:
-                                // <code>FESubfaceValues</code>).
-                                //
-                                // If we considered globally refined
-                                // meshes then only case 3 would
-                                // occur. But as we consider also
-                                // locally refined meshes we need to
-                                // distinguish all four cases making
-                                // the following assembling function
-                                // a bit longish.
-template <int dim>
-void DGMethod<dim>::assemble_system1 () 
-{
-  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-  std::vector<unsigned int> dofs (dofs_per_cell);
-  std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
-
-                                  // First we create the
-                                  // <code>update_flags</code> for the
-                                  // <code>FEValues</code> and the
-                                  // <code>FEFaceValues</code> objects.
-  const UpdateFlags update_flags = update_values
-                                   | update_gradients
-                                   | update_quadrature_points
-                                   | update_JxW_values;
-
-                                  // Note, that on faces we do not
-                                  // need gradients but we need
-                                  // normal vectors.
-  const UpdateFlags face_update_flags = update_values
-                                        | update_quadrature_points
-                                        | update_JxW_values
-                                        | update_normal_vectors;
-  
-                                  // On the neighboring cell we only
-                                  // need the shape values. Given a
-                                  // specific face, the quadrature
-                                  // points and `JxW values' are the
-                                  // same as for the current cells,
-                                  // the normal vectors are known to
-                                  // be the negative of the normal
-                                  // vectors of the current cell.
-  const UpdateFlags neighbor_face_update_flags = update_values;
-   
-                                  // Then we create the <code>FEValues</code>
-                                  // object. Note, that since version
-                                  // 3.2.0 of deal.II the constructor
-                                  // of this class takes a
-                                  // <code>Mapping</code> object as first
-                                  // argument. Although the
-                                  // constructor without <code>Mapping</code>
-                                  // argument is still supported it
-                                  // is recommended to use the new
-                                  // constructor. This reduces the
-                                  // effect of `hidden magic' (the
-                                  // old constructor implicitely
-                                  // assumes a <code>MappingQ1</code> mapping)
-                                  // and makes it easier to change
-                                  // the mapping object later.
-  FEValues<dim> fe_v (
-    mapping, fe, quadrature, update_flags);
-  
-                                  // Similarly we create the
-                                  // <code>FEFaceValues</code> and
-                                  // <code>FESubfaceValues</code> objects for
-                                  // both, the current and the
-                                  // neighboring cell. Within the
-                                  // following nested loop over all
-                                  // cells and all faces of the cell
-                                  // they will be reinited to the
-                                  // current cell and the face (and
-                                  // subface) number.
-  FEFaceValues<dim> fe_v_face (
-    mapping, fe, face_quadrature, face_update_flags);
-  FESubfaceValues<dim> fe_v_subface (
-    mapping, fe, face_quadrature, face_update_flags);
-  FEFaceValues<dim> fe_v_face_neighbor (
-    mapping, fe, face_quadrature, neighbor_face_update_flags);
-  FESubfaceValues<dim> fe_v_subface_neighbor (
-    mapping, fe, face_quadrature, neighbor_face_update_flags);
-
-                                  // Now we create the cell matrices
-                                  // and vectors. Here we need two
-                                  // cell matrices, both for face
-                                  // terms that include test
-                                  // functions <code>vi</code> (internal shape
-                                  // functions, i.e. shape functions
-                                  // of the current cell). To be more
-                                  // precise, the first matrix will
-                                  // include the `ui and vi terms'
-                                  // and the second will include the
-                                  // `ue and vi terms'. Here we
-                                  // recall the convention that `ui'
-                                  // is the shortcut for $u_h^+$ and
-                                  // `ue' represents $u_h^-$, see the
-                                  // introduction.
-  FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);
-
-  Vector<double>  cell_vector (dofs_per_cell);
-
-                                  // Furthermore we need some cell
-                                  // iterators.
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-
-                                  // Now we start the loop over all
-                                  // active cells.
-  for (;cell!=endc; ++cell) 
-    {
-                                      // In the
-                                      // <code>assemble_face_term1</code>
-                                      // function contributions to
-                                      // the cell matrices and the
-                                      // cell vector are only
-                                      // ADDED. Therefore on each
-                                      // cell we need to reset the
-                                      // <code>ui_vi_matrix</code> and
-                                      // <code>cell_vector</code> to zero,
-                                      // before assembling the cell terms.
-      ui_vi_matrix = 0;
-      cell_vector = 0;
-      
-                                      // Now we reinit the <code>FEValues</code>
-                                      // object for the current cell
-      fe_v.reinit (cell);
-
-                                      // and call the function
-                                      // that assembles the cell
-                                      // terms. The first argument is
-                                      // the <code>FEValues</code> that was
-                                      // previously reinit'ed on the
-                                      // current cell.
-      dg.assemble_cell_term(fe_v,
-                           ui_vi_matrix,
-                           cell_vector);
-
-                                      // As in previous examples the
-                                      // vector `dofs' includes the
-                                      // dof_indices.
-      cell->get_dof_indices (dofs);
-
-                                      // This is the start of the
-                                      // nested loop over all faces.
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-       {
-                                          // First we set the face
-                                          // iterator
-         typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
-         
-                                          // and clear the
-                                          // <code>ue_vi_matrix</code> on each
-                                          // face.
-         ue_vi_matrix = 0;
-                 
-                                          // Now we distinguish the
-                                          // four different cases in
-                                          // the ordering mentioned
-                                          // above. We start with
-                                          // faces belonging to the
-                                          // boundary of the domain.
-         if (face->at_boundary())
-           {
-                                              // We reinit the
-                                              // <code>FEFaceValues</code>
-                                              // object to the
-                                              // current face
-             fe_v_face.reinit (cell, face_no);
-
-                                              // and assemble the
-                                              // corresponding face
-                                              // terms.
-             dg.assemble_boundary_term(fe_v_face,
-                                       ui_vi_matrix,
-                                       cell_vector);
-           }
-         else
-           {
-                                              // Now we are not on
-                                              // the boundary of the
-                                              // domain, therefore
-                                              // there must exist a
-                                              // neighboring cell.
-             typename DoFHandler<dim>::cell_iterator neighbor=
-               cell->neighbor(face_no);;
-
-                                              // We proceed with the
-                                              // second and most
-                                              // complicated case:
-                                              // the neighboring cell
-                                              // is more refined than
-                                              // the current cell. As
-                                              // in deal.II
-                                              // neighboring cells
-                                              // are restricted to
-                                              // have a level
-                                              // difference of not
-                                              // more than one, the
-                                              // neighboring cell is
-                                              // known to be at most
-                                              // ONCE more refined
-                                              // than the current
-                                              // cell. Furthermore
-                                              // also the face is
-                                              // more refined,
-                                              // i.e. it has
-                                              // children. Here we
-                                              // note that the
-                                              // following part of
-                                              // code will not work
-                                              // for <code>dim==1</code>.
-             if (face->has_children())
-               {
-                                                  // First we store
-                                                  // which number the
-                                                  // current cell has
-                                                  // in the list of
-                                                  // neighbors of the
-                                                  // neighboring
-                                                  // cell. Hence,
-                                                  // neighbor-@>neighbor(neighbor2)
-                                                  // equals the
-                                                  // current cell
-                                                  // <code>cell</code>.
-                 const unsigned int neighbor2=
-                   cell->neighbor_of_neighbor(face_no);
-                 
-                 
-                                                  // We loop over
-                                                  // subfaces
-                 for (unsigned int subface_no=0;
-                      subface_no<face->n_children(); ++subface_no)
-                   {
-                                                      // and set the
-                                                      // cell
-                                                      // iterator
-                                                      // <code>neighbor_child</code>
-                                                      // to the cell
-                                                      // placed
-                                                      // `behind' the
-                                                      // current
-                                                      // subface.
-                     typename DoFHandler<dim>::active_cell_iterator
-                        neighbor_child
-                        = cell->neighbor_child_on_subface (face_no, subface_no);
-                     
-                     Assert (!neighbor_child->has_children(), ExcInternalError());
-
-                                                      // We need to
-                                                      // reset the
-                                                      // <code>ue_vi_matrix</code>
-                                                      // on each
-                                                      // subface
-                                                      // because on
-                                                      // each subface
-                                                      // the <code>un</code>
-                                                      // belong to
-                                                      // different
-                                                      // neighboring
-                                                      // cells.
-                     ue_vi_matrix = 0;
-                     
-                                                      // As already
-                                                      // mentioned
-                                                      // above for
-                                                      // the current
-                                                      // case (case
-                                                      // 2) we employ
-                                                      // the
-                                                      // <code>FESubfaceValues</code>
-                                                      // of the
-                                                      // current
-                                                      // cell (here
-                                                      // reinited for
-                                                      // the current
-                                                      // cell, face
-                                                      // and subface)
-                                                      // and we
-                                                      // employ the
-                                                      // FEFaceValues
-                                                      // of the
-                                                      // neighboring
-                                                      // child cell.
-                     fe_v_subface.reinit (cell, face_no, subface_no);
-                     fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
-
-                     dg.assemble_face_term1(fe_v_subface,
-                                            fe_v_face_neighbor,
-                                            ui_vi_matrix,
-                                            ue_vi_matrix);
-                     
-                                                      // Then we get
-                                                      // the dof
-                                                      // indices of
-                                                      // the
-                                                      // neighbor_child
-                                                      // cell
-                     neighbor_child->get_dof_indices (dofs_neighbor);
-                                                               
-                                                      // and
-                                                      // distribute
-                                                      // <code>ue_vi_matrix</code>
-                                                      // to the
-                                                      // system_matrix
-                     for (unsigned int i=0; i<dofs_per_cell; ++i)
-                       for (unsigned int k=0; k<dofs_per_cell; ++k)
-                         system_matrix.add(dofs[i], dofs_neighbor[k],
-                                           ue_vi_matrix(i,k));
-                   }
-                                                  // End of <code>if
-                                                  // (face-@>has_children())</code>
-               }
-             else
-               {
-                                                  // We proceed with
-                                                  // case 3,
-                                                  // i.e. neighboring
-                                                  // cell is of the
-                                                  // same refinement
-                                                  // level as the
-                                                  // current cell.
-                 if (neighbor->level() == cell->level()) 
-                   {
-                                                      // Like before
-                                                      // we store
-                                                      // which number
-                                                      // the current
-                                                      // cell has in
-                                                      // the list of
-                                                      // neighbors of
-                                                      // the
-                                                      // neighboring
-                                                      // cell.
-                     const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
-                                                      // We reinit
-                                                      // the
-                                                      // <code>FEFaceValues</code>
-                                                      // of the
-                                                      // current and
-                                                      // neighboring
-                                                      // cell to the
-                                                      // current face
-                                                      // and assemble
-                                                      // the
-                                                      // corresponding
-                                                      // face terms.
-                     fe_v_face.reinit (cell, face_no);
-                     fe_v_face_neighbor.reinit (neighbor, neighbor2);
-                     
-                     dg.assemble_face_term1(fe_v_face,
-                                            fe_v_face_neighbor,
-                                            ui_vi_matrix,
-                                            ue_vi_matrix);
-                                                      // End of <code>if
-                                                      // (neighbor-@>level()
-                                                      // ==
-                                                      // cell-@>level())</code>
-                   }
-                 else
-                   {
-                                                      // Finally we
-                                                      // consider
-                                                      // case 4. When
-                                                      // the
-                                                      // neighboring
-                                                      // cell is not
-                                                      // finer and
-                                                      // not of the
-                                                      // same
-                                                      // refinement
-                                                      // level as the
-                                                      // current cell
-                                                      // it must be
-                                                      // coarser.
-                     Assert(neighbor->level() < cell->level(), ExcInternalError());
-
-                                                      // Find out the
-                                                      // how many'th
-                                                      // face_no and
-                                                      // subface_no
-                                                      // the current
-                                                      // face is
-                                                      // w.r.t. the
-                                                      // neighboring
-                                                      // cell.
-                     const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                       cell->neighbor_of_coarser_neighbor(face_no);
-                     const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                     neighbor_subface_no=faceno_subfaceno.second;
-
-                     Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
-                                                                   neighbor_subface_no)
-                              == cell,
-                              ExcInternalError());
-
-                                                      // Reinit the
-                                                      // appropriate
-                                                      // <code>FEFaceValues</code>
-                                                      // and assemble
-                                                      // the face
-                                                      // terms.
-                     fe_v_face.reinit (cell, face_no);
-                     fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
-                                                   neighbor_subface_no);
-                     
-                     dg.assemble_face_term1(fe_v_face,
-                                            fe_v_subface_neighbor,
-                                            ui_vi_matrix,
-                                            ue_vi_matrix);
-                   }
-
-                                                  // Now we get the
-                                                  // dof indices of
-                                                  // the
-                                                  // <code>neighbor_child</code>
-                                                  // cell,
-                 neighbor->get_dof_indices (dofs_neighbor);
-                                                               
-                                                  // and distribute the
-                                                  // <code>ue_vi_matrix</code>.
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   for (unsigned int k=0; k<dofs_per_cell; ++k)
-                     system_matrix.add(dofs[i], dofs_neighbor[k],
-                                       ue_vi_matrix(i,k));
-               }
-                                              // End of <code>face not at boundary</code>:
-           }
-                                          // End of loop over all faces:
-       }
-      
-                                      // Finally we distribute the
-                                      // <code>ui_vi_matrix</code>
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
-      
-                                      // and the cell vector.
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       right_hand_side(dofs[i]) += cell_vector(i);
-    }
-}
-
-
-                                // @sect4{Function: assemble_system2}
-                                //
-                                // We proceed with the
-                                // <code>assemble_system2</code> function that
-                                // implements the DG discretization
-                                // in its second version. This
-                                // function is very similar to the
-                                // <code>assemble_system1</code>
-                                // function. Therefore, here we only
-                                // discuss the differences between
-                                // the two functions. This function
-                                // repeatedly calls the
-                                // <code>assemble_face_term2</code> function
-                                // of the DGTransportEquation object,
-                                // that assembles the face terms
-                                // written as a sum of integrals over
-                                // all faces. Therefore, we need to
-                                // make sure that each face is
-                                // treated only once. This is achieved
-                                // by introducing the rule:
-                                // 
-                                // a) If the current and the
-                                // neighboring cells are of the same
-                                // refinement level we access and
-                                // treat the face from the cell with
-                                // lower index.
-                                //
-                                // b) If the two cells are of
-                                // different refinement levels we
-                                // access and treat the face from the
-                                // coarser cell.
-                                //
-                                // Due to rule b) we do not need to
-                                // consider case 4 (neighboring cell
-                                // is coarser) any more.
-
+                                // @sect4{Function: assemble_system}
+                                // Here we see the major difference
+                                // to assembling by hand. Instead of
+                                // writing loops over cells and
+                                // faces, we leave all this to the
+                                // MeshWorker framework. In order to
+                                // do so, we just have to define
+                                // local integration objects and use
+                                // one of the classes in Assembler to
+                                // build the global system. 
 template <int dim>
-void DGMethod<dim>::assemble_system2 () 
+void DGMethod<dim>::assemble_system () 
 {
-  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-  std::vector<unsigned int> dofs (dofs_per_cell);
-  std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
-
-  const UpdateFlags update_flags = update_values
-                                   | update_gradients
-                                   | update_quadrature_points
-                                   | update_JxW_values;
-  
-  const UpdateFlags face_update_flags = update_values
-                                        | update_quadrature_points
-                                        | update_JxW_values
-                                        | update_normal_vectors;
-  
-  const UpdateFlags neighbor_face_update_flags = update_values;
-
-                                  // Here we do not need
-                                  // <code>fe_v_face_neighbor</code> as case 4
-                                  // does not occur.
-  FEValues<dim> fe_v (
-    mapping, fe, quadrature, update_flags);
-  FEFaceValues<dim> fe_v_face (
-    mapping, fe, face_quadrature, face_update_flags);
-  FESubfaceValues<dim> fe_v_subface (
-    mapping, fe, face_quadrature, face_update_flags);
-  FEFaceValues<dim> fe_v_face_neighbor (
-    mapping, fe, face_quadrature, neighbor_face_update_flags);
-
-
-  FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);
-  
-                                  // Additionally we need the
-                                  // following two cell matrices,
-                                  // both for face term that include
-                                  // test function <code>ve</code> (external
-                                  // shape functions, i.e. shape
-                                  // functions of the neighboring
-                                  // cell). To be more precise, the
-                                  // first matrix will include the `u
-                                  // and vn terms' and the second
-                                  // that will include the `un and vn
-                                  // terms'.
-  FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> ue_ve_matrix (dofs_per_cell, dofs_per_cell);
+                                  // Here we generate an object of
+                                  // our own integration class, which
+                                  // knows how to compute cell and
+                                  // face contributions for the
+                                  // matrix and the residual.
+  const DGIntegrator<dim> dg;
+
+                                  // This is the magic object, which
+                                  // knows everything about the data
+                                  // structures and local
+                                  // integration (the latter through
+                                  // our object @p dg). This is the
+                                  // object doing the work in the
+                                  // function MeshWorker::loop, which
+                                  // is implicitly called
+                                  // below. After @p dg did the local
+                                  // integration, the
+                                  // MeshWorker::Assembler::SystemSimple
+                                  // object distributes these into
+                                  // the global sparse matrix and the
+                                  // right hand side vector.
+  MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> >, DGIntegrator<dim> >
+    integrator(dg);
   
-  Vector<double>  cell_vector (dofs_per_cell);
-
-                                  // The following lines are roughly
-                                  // the same as in the previous
-                                  // function.
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (;cell!=endc; ++cell) 
-    {
-      ui_vi_matrix = 0;
-      cell_vector = 0;
-
-      fe_v.reinit (cell);
-
-      dg.assemble_cell_term(fe_v,
-                           ui_vi_matrix,
-                           cell_vector);
-      
-      cell->get_dof_indices (dofs);
-
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-       {
-         typename DoFHandler<dim>::face_iterator face=
-           cell->face(face_no);
-
-                                          // Case 1:
-         if (face->at_boundary())
-           {
-             fe_v_face.reinit (cell, face_no);
-
-             dg.assemble_boundary_term(fe_v_face,
-                                       ui_vi_matrix,
-                                       cell_vector);
-           }
-         else
-           {
-             Assert (cell->neighbor(face_no).state() == IteratorState::valid,
-                     ExcInternalError());
-             typename DoFHandler<dim>::cell_iterator neighbor=
-               cell->neighbor(face_no);
-                                              // Case 2:
-             if (face->has_children())
-               {
-                 const unsigned int neighbor2=
-                   cell->neighbor_of_neighbor(face_no);
-                 
-                 for (unsigned int subface_no=0;
-                      subface_no<face->n_children(); ++subface_no)
-                   {
-                     typename DoFHandler<dim>::cell_iterator neighbor_child
-                        = cell->neighbor_child_on_subface (face_no, subface_no);
-                     Assert (!neighbor_child->has_children(), ExcInternalError());
-                     
-                     ue_vi_matrix = 0;
-                     ui_ve_matrix = 0;
-                     ue_ve_matrix = 0;
-                     
-                     fe_v_subface.reinit (cell, face_no, subface_no);
-                     fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
-
-                     dg.assemble_face_term2(fe_v_subface,
-                                            fe_v_face_neighbor,
-                                            ui_vi_matrix,
-                                            ue_vi_matrix,
-                                            ui_ve_matrix,
-                                            ue_ve_matrix);
-                 
-                     neighbor_child->get_dof_indices (dofs_neighbor);
-                                                               
-                     for (unsigned int i=0; i<dofs_per_cell; ++i)
-                       for (unsigned int j=0; j<dofs_per_cell; ++j)
-                         {
-                           system_matrix.add(dofs[i], dofs_neighbor[j],
-                                             ue_vi_matrix(i,j));
-                           system_matrix.add(dofs_neighbor[i], dofs[j],
-                                             ui_ve_matrix(i,j));
-                           system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
-                                             ue_ve_matrix(i,j));
-                         }
-                   }
-               }
-             else
-               {
-                                                  // Case 3, with the
-                                                  // additional rule
-                                                  // a)
-                 if (neighbor->level() == cell->level() &&
-                     neighbor->index() > cell->index()) 
-                   {
-                     const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-                     
-                     ue_vi_matrix = 0;
-                     ui_ve_matrix = 0;
-                     ue_ve_matrix = 0;
-                     
-                     fe_v_face.reinit (cell, face_no);
-                     fe_v_face_neighbor.reinit (neighbor, neighbor2);
-                     
-                     dg.assemble_face_term2(fe_v_face,
-                                            fe_v_face_neighbor,
-                                            ui_vi_matrix,
-                                            ue_vi_matrix,
-                                            ui_ve_matrix,
-                                            ue_ve_matrix);
-
-                     neighbor->get_dof_indices (dofs_neighbor);
-
-                     for (unsigned int i=0; i<dofs_per_cell; ++i)
-                       for (unsigned int j=0; j<dofs_per_cell; ++j)
-                         {
-                           system_matrix.add(dofs[i], dofs_neighbor[j],
-                                             ue_vi_matrix(i,j));
-                           system_matrix.add(dofs_neighbor[i], dofs[j],
-                                             ui_ve_matrix(i,j));
-                           system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
-                                             ue_ve_matrix(i,j));
-                         }
-                   }
-
-                                                  // Due to rule b)
-                                                  // we do not need
-                                                  // to consider case
-                                                  // 4.
-               }
-           }
-       }
-      
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
-      
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       right_hand_side(dofs[i]) += cell_vector(i);
-    }
+                                  // First, we initialize the
+                                  // quadrature formulae and the
+                                  // update flags in the worker base
+                                  // class. For quadrature, we play
+                                  // safe and use a QGauss forumla
+                                  // with number of points one higher
+                                  // than the polynomial degree
+                                  // used. Since the quadratures for
+                                  // cells, boundary and interior
+                                  // faces can be selected
+                                  // independently, we have to hand
+                                  // over this value three times.
+  const unsigned int n_gauss_points = dof_handler.get_fe().degree+1;
+  integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+
+                                  // These are the types of values we
+                                  // need for integrating our
+                                  // system. They are added to the
+                                  // flags used on cells, boundary
+                                  // and interior faces, as well as
+                                  // interior neighbor faces, which is
+                                  // forced by the four @p true values.
+  UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
+  integrator.add_update_flags(update_flags, true, true, true, true);
+
+                                  // Finally, we have to tell the
+                                  // assembler base class, where to
+                                  // put the local data. These will
+                                  // be our system matrix and the
+                                  // right hand side.
+  integrator.initialize(system_matrix, right_hand_side);
+
+                                  // Finally, we get to the
+                                  // integration loop. @p info_box is
+                                  // an object, that generates the
+                                  // extended iterators for cells and
+                                  // faces of type
+                                  // MeshWorker::IntegrationInfo. Since
+                                  // we need five different of them,
+                                  // this is a handy shortcut. It
+                                  // receives all the stuff we
+                                  // created so far.
+  MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
+  info_box.initialize(integrator, fe, mapping);
+
+                                  // Finally, the integration loop
+                                  // over all active cells
+                                  // (determined by the first
+                                  // argument, which is an active iterator).
+  MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator);
 }
 
 
@@ -1453,7 +748,7 @@ void DGMethod<dim>::refine_grid ()
                                   // are computed
   DerivativeApproximation::approximate_gradient (mapping,
                                                 dof_handler,
-                                                solution2,
+                                                solution,
                                                 gradient_indicator);
 
                                   // and they are cell-wise scaled by
@@ -1509,7 +804,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
   
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution2, "u");
+  data_out.add_data_vector (solution, "u");
 
   data_out.build_patches ();
   
@@ -1557,53 +852,9 @@ void DGMethod<dim>::run ()
                << dof_handler.n_dofs()
                << std::endl;
 
-                                      // The constructor of the Timer
-                                      // class automatically starts
-                                      // the time measurement.
-      Timer assemble_timer;
-                                      // First assembling routine.
-      assemble_system1 ();
-                                      // The operator () accesses the
-                                      // current time without
-                                      // disturbing the time
-                                      // measurement.
-      std::cout << "Time of assemble_system1: "
-               << assemble_timer()
-               << std::endl;
-      solve (solution1);
-
-                                      // As preparation for the
-                                      // second assembling routine we
-                                      // reinit the system matrix, the
-                                      // right hand side vector and
-                                      // the Timer object.
-      system_matrix = 0;
-      right_hand_side = 0;
-      assemble_timer.reset();
-
-                                      // We start the Timer,
-      assemble_timer.start();
-                                      // call the second assembling routine
-      assemble_system2 ();
-                                      // and access the current time.
-      std::cout << "Time of assemble_system2: "
-               << assemble_timer()
-               << std::endl;
-      solve (solution2);
-
-                                      // To make sure that both
-                                      // versions of the DG method
-                                      // yield the same
-                                      // discretization and hence the
-                                      // same solution we check the
-                                      // two solutions for equality.
-      solution1-=solution2;
-      const double difference=solution1.linfty_norm();
-      if (difference>1e-13)
-       std::cout << "solution1 and solution2 differ!!" << std::endl;
-      else
-       std::cout << "solution1 and solution2 coincide." << std::endl;
-       
+      assemble_system ();
+      solve (solution);
+
                                       // Finally we perform the
                                       // output.
       output_results (cycle);

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.