--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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>;
+}
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
//---------------------------------------------------------------------------
// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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>;
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $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)
+// {
+// }
+ }
+
+
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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>;
+}
// 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>
}
- // @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 ();
// 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];
// \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];
// 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];
// \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];
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;
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
// 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)
{}
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);
}
// are computed
DerivativeApproximation::approximate_gradient (mapping,
dof_handler,
- solution2,
+ solution,
gradient_indicator);
// and they are cell-wise scaled by
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 ();
<< 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);