From 44362bf8fb97e8cd407e95bb88ee2a791269d21f Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Fri, 2 Oct 2009 00:08:27 +0000 Subject: [PATCH] first example for the generic MeshWorker::loop() is up and running git-svn-id: https://svn.dealii.org/trunk@19654 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/numerics/mesh_worker.h | 411 ++++++ .../include/numerics/mesh_worker_assembler.h | 1308 +++++++++++++++++ .../include/numerics/mesh_worker_info.h | 874 +++++++++++ .../numerics/mesh_worker_info.templates.h | 192 +++ .../include/numerics/mesh_worker_loop.h | 188 +++ .../numerics/mesh_worker_vector_selector.h | 375 +++++ .../mesh_worker_vector_selector.inst.in | 17 + .../mesh_worker_vector_selector.templates.h | 86 ++ .../source/numerics/fe_field_function.inst.in | 3 +- .../deal.II/source/numerics/mesh_worker.cc | 116 ++ .../mesh_worker_assembler.all_dimensions.cc | 130 ++ .../source/numerics/mesh_worker_info.cc | 42 + .../source/numerics/mesh_worker_info.inst.in | 34 + .../numerics/mesh_worker_vector_selector.cc | 28 + .../mesh_worker_vector_selector.inst.in | 17 + deal.II/examples/step-38/step-38.cc | 1269 ++++------------ 16 files changed, 4079 insertions(+), 1011 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/mesh_worker.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_assembler.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_info.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_info.templates.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_loop.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_vector_selector.inst.in create mode 100644 deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h create mode 100644 deal.II/deal.II/source/numerics/mesh_worker.cc create mode 100644 deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc create mode 100644 deal.II/deal.II/source/numerics/mesh_worker_info.cc create mode 100644 deal.II/deal.II/source/numerics/mesh_worker_info.inst.in create mode 100644 deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc create mode 100644 deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in diff --git a/deal.II/deal.II/include/numerics/mesh_worker.h b/deal.II/deal.II/include/numerics/mesh_worker.h new file mode 100644 index 0000000000..df1e06709b --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker.h @@ -0,0 +1,411 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__mesh_worker_h +#define __deal2__mesh_worker_h + +#include +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +template class DoFHandler; +template 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. + * + *

Simplified interfaces

+ * + * 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 + class LocalWorker + { + public: + /** + * Constructor, setting + * #interior_fluxes and + * #boundary_fluxes to @p + * true. + */ + LocalWorker (); + + /** + * Do the work on a cell. + */ + void cell(DoFInfo& cell); + + /** + * Do the work on a boundary face. + */ + void bdry(DoFInfo& face); + + /** + * Do the work on an interior face. + */ + void face(DoFInfo& face1, DoFInfo& 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 + class IntegrationWorker : public LocalWorker + { + public: + /** + * The info type expected by a + * cell integrator. + */ + typedef IntegrationInfo, dim> CellInfo; + + /** + * The info type expected by a + * face integrator. + */ + typedef IntegrationInfo, 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 cell_quadrature; + + /** + * The quadrature rule used + * on boundary faces. + */ + Quadrature bdry_quadrature; + + /** + * The quadrature rule used + * on interior faces. + */ + Quadrature 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 + * class MyLocalOperator : public Subscriptor + * { + * public: + * void cell(typename IntegrationWorker::CellInfo& info) const; + * void bdry(typename IntegrationWorker::FaceInfo& info) const; + * void face(typename IntegrationWorker::FaceInfo& info1, + * typename IntegrationWorker::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 myop; + * + * AssemblingIntegrator > + * 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 + class AssemblingIntegrator : + public IntegrationWorker, + public ASSEMBLER + { + public: + /** + * Constructor, initializing + * the local data. + */ + AssemblingIntegrator(const INTEGRATOR& local); + /** + * The cell operator called by + * integration_loop(). + */ + void cell(typename IntegrationWorker::CellInfo& info); + + /** + * The local boundary operator + * called by + * integration_loop(). + */ + void bdry(typename IntegrationWorker::FaceInfo& info); + + /** + * The interior face operator + * called by + * integration_loop(). + */ + void face(typename IntegrationWorker::FaceInfo& info1, + typename IntegrationWorker::FaceInfo& info2); + + private: + /** + * Pointer to the object doing + * local integration. + */ + SmartPointer local; + }; + +//----------------------------------------------------------------------// + template + inline + LocalWorker::LocalWorker() + : + interior_fluxes(true), boundary_fluxes(true) + {} + +//----------------------------------------------------------------------// + + template + inline + AssemblingIntegrator::AssemblingIntegrator(const INTEGRATOR& local) + : + local(&local, typeid(*this).name()) + {} + + + template + inline void + AssemblingIntegrator::cell( + typename IntegrationWorker::CellInfo& info) + { + local->cell(info); + ASSEMBLER::assemble(info); + } + + + template + inline void + AssemblingIntegrator::bdry( + typename IntegrationWorker::FaceInfo& info) + { + local->bdry(info); + ASSEMBLER::assemble(info); + } + + + template + inline void + AssemblingIntegrator::face( + typename IntegrationWorker::FaceInfo& info1, + typename IntegrationWorker::FaceInfo& info2) + { + local->face(info1, info2); + ASSEMBLER::assemble(info1, info2); + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h new file mode 100644 index 0000000000..b71701d1a8 --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -0,0 +1,1308 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__mesh_worker_assembler_h +#define __deal2__mesh_worker_assembler_h + +#include +#include +#include +#include +#include + + +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. + * + *

Data models

+ * + * 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: + * + *

The comprehensive data model

+ * + * 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. + * + *

The blocked data model

+ * + * 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 + 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 + void initialize_info(DoFInfo& info, + bool interior_face); + + /** + * Assemble the local values + * into the global vectors. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local values + * into the global vectors. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + + /** + * The value of the ith entry + * in #results. + */ + number operator() (unsigned int i) const; + private: + /** + * The values into which the + * results are added. + */ + std::vector results; + }; + +/** + * Compute cell and face contributions of one or several functionals, + * typically for error estimates. + * + * @author Guido Kanschat, 2009 + */ + template + class CellsAndFaces + { + public: + /** + * The data type for + * communicating the cell and + * face vectors. + */ + typedef NamedData > > 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 + void initialize_info(DoFInfo& info, + bool interior_face); + + /** + * Assemble the local values + * into the global vectors. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local values + * into the global vectors. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + + /** + * The value of the ith entry + * in #results. + */ + number operator() (unsigned int i) const; + private: + DataVectors results; + bool separate_faces; + }; + + +/** + * Assemble residuals without block structure. + * + * The data structure for this Assembler class is a simple vector on + * each cell with entries from zero to + * FiniteElementData::dofs_per_cell and a simple global vector with + * entries numbered from zero to DoFHandler::n_dofs(). No BlockInfo is + * required and the global vector may be any type of vector having + * element access through operator() (unsigned int) + */ + template + class ResidualSimple + { + public: + void initialize(NamedData >& 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 + void initialize_info(DoFInfo& info, + bool interior_face) const; + + /** + * Assemble the local residuals + * into the global residuals. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local residuals + * into the global residuals. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + private: + /** + * The global residal vectors + * filled by assemble(). + */ + NamedData > 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 ResidualLocalBlocksToGlobalBlocks + { + public: + /** + * Copy the BlockInfo and the + * matrix pointers into local + * variables. + */ + void initialize(NamedData >& 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 + void initialize_info(DoFInfo& info, + bool interior_face) const; + + + /** + * Assemble the local residuals + * into the global residuals. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local residuals + * into the global residuals. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + private: + /** + * Assemble a single local + * residual into the global. + */ + void assemble(VECTOR global, + const BlockVector& local, + const BlockIndices& bi, + const std::vector& dof); + + /** + * The global matrices, + * stored as a vector of + * pointers. + */ + NamedData > residuals; + }; + + +/** + * Assemble local matrices into a single global matrix without using + * block structure. + * + * @author Guido Kanschat, 2009 + */ + template + 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 + void initialize_info(DoFInfo& info, + bool interior_face) const; + + /** + * Assemble the matrix + * DoFInfo::M1[0] + * into the global matrix. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local + * matrices in the info + * objects into the global + * matrix. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + private: + /** + * Assemble a single matrix + * into #matrix. + */ + void assemble(const FullMatrix& M, + const std::vector& i1, + const std::vector& i2); + + /** + * The global matrix being + * assembled. + */ + SmartPointer 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 MatrixLocalBlocksToGlobalBlocks + { + public: + /// The object that is stored + typedef boost::shared_ptr > 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& 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 + void initialize_info(DoFInfo& info, + bool interior_face) const; + + + /** + * Assemble the local matrices + * into the global matrices. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble all local matrices + * into the global matrices. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + + private: + /** + * Assemble a single local + * matrix into a global one. + */ + void assemble( + MATRIX& global, + const FullMatrix& local, + unsigned int block_row, + unsigned int block_col, + const std::vector& dof1, + const std::vector& dof2); + + /** + * The global matrices, + * stored as a vector of + * pointers. + */ + std::vector 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 MGMatrixLocalBlocksToGlobalBlocks + { + public: + /// The object that is stored + typedef boost::shared_ptr > > 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& 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& up, std::vector& 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 + void initialize_info(DoFInfo& info, + bool interior_face) const; + + + /** + * Assemble the local matrices + * into the global matrices. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble all local matrices + * into the global matrices. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + + private: + /** + * Assemble a single local + * matrix into a global one. + */ + void assemble( + MATRIX& global, + const FullMatrix& local, + unsigned int block_row, + unsigned int block_col, + const std::vector& dof1, + const std::vector& dof2, + unsigned int level1, + unsigned int level2, + bool transpose = false); + + /** + * The level matrices, + * stored as a vector of + * pointers. + */ + std::vector matrices; + + /** + * The flux matrix between + * the fine and the coarse + * level at refinement edges. + */ + std::vector flux_down; + + /** + * The flux matrix between + * the coarse and the fine + * level at refinement edges. + */ + std::vector 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 SystemSimple : + private MatrixSimple, + private ResidualSimple + { + 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 + void initialize_info(DoFInfo& info, bool) const; + + /** + * Assemble the matrix + * DoFInfo::M1[0] + * into the global matrix. + */ + template + void assemble(const DoFInfo& info); + + /** + * Assemble both local + * matrices in the info + * objects into the global + * matrix. + */ + template + void assemble(const DoFInfo& info1, + const DoFInfo& info2); + }; + + +//----------------------------------------------------------------------// + + template + inline void + ResidualSimple::initialize(NamedData >& results) + { + residuals = results; + } + + + template + template + inline void + ResidualSimple::initialize_info(DoFInfo& info, bool) const + { + info.initialize_vectors(residuals.size()); + } + + + template + template + inline void + ResidualSimple::assemble(const DoFInfo& info) + { + for (unsigned int k=0;k + template + inline void + ResidualSimple::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int k=0;k + inline void + ResidualLocalBlocksToGlobalBlocks::initialize(NamedData >& m) + { + residuals = m; + } + + + template + template + inline void + ResidualLocalBlocksToGlobalBlocks::initialize_info( + DoFInfo& info, bool) const + { + info.initialize_vectors(residuals.size()); + } + + template + inline void + ResidualLocalBlocksToGlobalBlocks::assemble( + VECTOR global, + const BlockVector& local, + const BlockIndices& bi, + const std::vector& dof) + { + for (unsigned int b=0;bbi.local_to_global(b, j); + (*global)(dof[jcell]) += local.block(b)(j); + } + } + + + template + template + inline void + ResidualLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info) + { + for (unsigned int i=0;i + template + inline void + ResidualLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;i + inline + MatrixSimple::MatrixSimple(double threshold) + : + threshold(threshold) + {} + + + template + inline void + MatrixSimple::initialize(MATRIX& m) + { + matrix = &m; + } + + + template + template + inline void + MatrixSimple::initialize_info(DoFInfo& info, + bool interior_face) const + { + info.initialize_matrices(1, interior_face); + } + + + + template + inline void + MatrixSimple::assemble(const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + for (unsigned int j=0; j= threshold) + matrix->add(i1[j], i2[k], M(j,k)); + } + + + template + template + inline void + MatrixSimple::assemble(const DoFInfo& info) + { + assemble(info.M1[0].matrix, info.indices, info.indices); + } + + + template + template + inline void + MatrixSimple::assemble(const DoFInfo& info1, + const DoFInfo& 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 + inline + MatrixLocalBlocksToGlobalBlocks::MatrixLocalBlocksToGlobalBlocks( + double threshold) + : + threshold(threshold) + {} + + + template + inline void + MatrixLocalBlocksToGlobalBlocks::initialize(std::vector& m) + { + matrices = m; + } + + + + template + inline void + MatrixLocalBlocksToGlobalBlocks::assemble( + MATRIX& global, + const FullMatrix& local, + unsigned int block_row, + unsigned int block_col, + const std::vector& dof1, + const std::vector& dof2) + { + for (unsigned int j=0;j= 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 + template + inline void + MatrixLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info) + { + for (unsigned int i=0;irow; + const unsigned int col = matrices[i]->column; + + assemble(matrices[i]->matrix, info->M1[i].matrix, row, col, info.indices, info.indices); + } + } + + + template + template + inline void + MatrixLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;irow; + 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 + inline + MGMatrixLocalBlocksToGlobalBlocks::MGMatrixLocalBlocksToGlobalBlocks( + double threshold) + : + threshold(threshold) + {} + + + template + inline void + MGMatrixLocalBlocksToGlobalBlocks::initialize(std::vector& m) + { + matrices = m; + } + + + template + inline void + MGMatrixLocalBlocksToGlobalBlocks::initialize_edge_flux( + std::vector& up, + std::vector& down) + { + flux_up = up; + flux_down = down; + } + + + template + inline void + MGMatrixLocalBlocksToGlobalBlocks::assemble( + MATRIX& global, + const FullMatrix& local, + unsigned int block_row, + unsigned int block_col, + const std::vector& dof1, + const std::vector& dof2, + unsigned int level1, + unsigned int level2, + bool transpose) + { + for (unsigned int j=0;j= 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 + template + inline void + MGMatrixLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info) + { + const unsigned int level = info.cell->level(); + + for (unsigned int i=0;irow; + 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 + template + inline void + MGMatrixLocalBlocksToGlobalBlocks::assemble( + const DoFInfo& info1, + const DoFInfo& info2) + { + const unsigned int level1 = info1.cell->level(); + const unsigned int level2 = info2.cell->level(); + + for (unsigned int i=0;irow; + 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 + SystemSimple::SystemSimple(double t) + : + MatrixSimple(t) + {} + + + template + inline void + SystemSimple::initialize(MATRIX& m, VECTOR& rhs) + { + NamedData > data; + SmartPointer p = &rhs; + data.add(p, "right hand side"); + + MatrixSimple::initialize(m); + ResidualSimple::initialize(data); + } + + + template + template + inline void + SystemSimple::initialize_info(DoFInfo& info, + bool interior_face) const + { + MatrixSimple::initialize_info(info, interior_face); + ResidualSimple::initialize_info(info, interior_face); + } + + + template + template + inline void + SystemSimple::assemble(const DoFInfo& info) + { + MatrixSimple::assemble(info); + ResidualSimple::assemble(info); + } + + + template + template + inline void + SystemSimple::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + MatrixSimple::assemble(info1, info2); + ResidualSimple::assemble(info1, info2); + } + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.h b/deal.II/deal.II/include/numerics/mesh_worker_info.h new file mode 100644 index 0000000000..c0e8e5213c --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -0,0 +1,874 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__mesh_worker_info_h +#define __deal2__mesh_worker_info_h + +#include +#include +#include +#include +#include +#include + +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 + class LocalResults + { + /** + * Initialize a single local + * matrix block. A helper + * function for initialize() + */ + void initialize_local( + MatrixBlock >& 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 + void initialize_matrices(std::vector& 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 J; + + /** + * The local vectors. This + * field is public, so that + * local integrators can + * write it. + */ + std::vector > R; + + /** + * The local matrices + * coupling degrees of + * freedom in the cell + * itself or within the + * first cell on a face. + */ + std::vector > > M1; + + /** + * The local matrices + * coupling test functions on + * the cell with trial + * functions on the other + * cell. + * + * Only used on interior + * faces. + */ + std::vector > > 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 + *
    + *
  • LocalVectors + *
  • LocalMatrices + *
+ * + * 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 + class DoFInfo : public LocalResults + { + public: + /// The current cell + typename Triangulation::cell_iterator cell; + + /// The current face + typename Triangulation::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 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 + DoFInfo(const DH& dof_handler); + + /** + * Set the current cell and + * fill #indices. + */ + template + void reinit(const DHCellIterator& c); + + /** + * Set the current face and + * fill #indices if the #cell + * changed. + */ + template + void reinit(const DHCellIterator& c, + const DHFaceIterator& f, + unsigned int n); + + /** + * Set the current subface + * and fill #indices if the + * #cell changed. + */ + template + 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 block_info; + private: + /// Fill index vector + void get_indices(const typename DoFHandler::cell_iterator c); + + /// Fill index vector with level indices + void get_indices(const typename MGDoFHandler::cell_iterator c); + + /// Auxiliary vector + std::vector 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. + * + *

Integration models

+ * + * 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 + class IntegrationInfo : public DoFInfo + { + /// vector of FEValues objects + std::vector > fevalv; + public: + /** + * Constructor forwarding + * information to DoFInfo. + */ + IntegrationInfo(const BlockInfo& block_info); + + /** + * Constructor forwarding + * information to DoFInfo. + */ + template + 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 + void initialize(const FiniteElement& el, + const Mapping& mapping, + const Quadrature& quadrature, + const UpdateFlags flags); + + /** + * Initialize the data + * vector and cache the + * selector. + */ + void initialize_data(const boost::shared_ptr > 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 > > 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 > > > 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 > > > hessians; + + /** + * Reinitialize internal data + * structures for use on a cell. + */ + template + void reinit(const DHCellIterator& c); + + /** + * Reinitialize internal data + * structures for use on a face. + */ + template + void reinit(const DHCellIterator& c, + const DHFaceIterator& f, + const unsigned int fn); + + /** + * Reinitialize internal data + * structures for use on a subface. + */ + template + 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 > 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 + void fill_local_data( + std::vector > >& 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 + class IntegrationInfoBox + { + public: + typedef IntegrationInfo, spacedim> CellInfo; + typedef IntegrationInfo, spacedim> FaceInfo; + + /** + * Initialize all members + * with the same argument. + */ + template + IntegrationInfoBox(const T&); + + template + void initialize(const WORKER&, + const FiniteElement& el, + const Mapping& mapping); + + template + void initialize(const WORKER&, + const FiniteElement& el, + const Mapping& mapping, + const NamedData >& data); +// private: + + CellInfo cell_info; + FaceInfo bdry_info; + FaceInfo face_info; + FaceInfo subface_info; + FaceInfo neighbor_info; + }; + + +//----------------------------------------------------------------------// + + template + inline void + LocalResults::initialize_vectors(unsigned int n) + { + R.resize(n); + } + + + template + template + inline void + LocalResults::initialize_matrices( + std::vector& matrices, + bool both) + { + M1.resize(matrices.size()); + if (both) + M2.resize(matrices.size()); + for (unsigned int i=0;irow; + 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 + inline void + LocalResults::initialize_matrices(unsigned int n, bool both) + { + M1.resize(n); + if (both) + M2.resize(n); + for (unsigned int i=0;i + inline void + LocalResults::reinit(const BlockIndices& bi) + { + for (unsigned int i=0;i + template + DoFInfo::DoFInfo(const DH& dof_handler) + { + std::vector aux(1); + aux[0] = dof_handler.get_fe().dofs_per_cell; + aux_local_indices.reinit(aux); + } + + + template + template + inline void + DoFInfo::reinit(const DHCellIterator& c) + { + get_indices(c); + cell = static_cast::cell_iterator> (c); + face_number = deal_II_numbers::invalid_unsigned_int; + sub_number = deal_II_numbers::invalid_unsigned_int; + if (block_info) + LocalResults::reinit(block_info->local()); + else + LocalResults::reinit(aux_local_indices); + } + + + template + template + inline void + DoFInfo::reinit(const DHCellIterator& c, + const DHFaceIterator& f, + unsigned int n) + { + if ((cell.state() != IteratorState::valid) + || cell != static_cast::cell_iterator> (c)) + get_indices(c); + cell = static_cast::cell_iterator> (c); + face = static_cast::face_iterator> (f); + face_number = n; + sub_number = deal_II_numbers::invalid_unsigned_int; + if (block_info) + LocalResults::reinit(block_info->local()); + else + LocalResults::reinit(aux_local_indices); + } + + + template + template + inline void + DoFInfo::reinit(const DHCellIterator& c, + const DHFaceIterator& f, + unsigned int n, + const unsigned int s) + { + if (cell.state() != IteratorState::valid + || cell != static_cast::cell_iterator> (c)) + get_indices(c); + cell = static_cast::cell_iterator> (c); + face = static_cast::face_iterator> (f); + face_number = n; + sub_number = s; + if (block_info) + LocalResults::reinit(block_info->local()); + else + LocalResults::reinit(aux_local_indices); + } + + + template + inline const BlockIndices& + DoFInfo::local_indices() const + { + if (block_info) + return block_info->local(); + return aux_local_indices; + } + + +//----------------------------------------------------------------------// + + template + template + IntegrationInfo::IntegrationInfo(const DH& dof_handler) + : + DoFInfo(dof_handler), + fevalv(0), + multigrid(false), + global_data(boost::shared_ptr >(new VectorDataBase)) + {} + + + template + inline const FVB& + IntegrationInfo::fe() const + { + AssertDimension(fevalv.size(), 1); + return *fevalv[0]; + } + + + template + inline const FVB& + IntegrationInfo::fe(unsigned int i) const + { + Assert (i + template + inline void + IntegrationInfo::reinit(const DHCellIterator& c) + { + DoFInfo::reinit(c); + for (unsigned int i=0;i& fe = dynamic_cast&> (febase); + fe.reinit(this->cell); + } + + const bool split_fevalues = this->block_info != 0; + fill_local_data(split_fevalues); + } + + + template + template + inline void + IntegrationInfo::reinit( + const DHCellIterator& c, + const DHFaceIterator& f, + const unsigned int fn) + { + DoFInfo::reinit(c, f, fn); + for (unsigned int i=0;i& fe = dynamic_cast&> (febase); + fe.reinit(this->cell, fn); + } + + const bool split_fevalues = this->block_info != 0; + fill_local_data(split_fevalues); + } + + + template + template + inline void + IntegrationInfo::reinit( + const DHCellIterator& c, + const DHFaceIterator& f, + const unsigned int fn, + const unsigned int sn) + { + DoFInfo::reinit(c, f, fn, sn); + for (unsigned int i=0;i& fe = dynamic_cast&> (febase); + fe.reinit(this->cell, fn, sn); + } + + const bool split_fevalues = this->block_info != 0; + fill_local_data(split_fevalues); + } + +//----------------------------------------------------------------------// + + template + template + IntegrationInfoBox::IntegrationInfoBox(const T& t) + : + cell_info(t), + bdry_info(t), + face_info(t), + subface_info(t), + neighbor_info(t) + {} + + template + template + void + IntegrationInfoBox::initialize( + const WORKER& integrator, + const FiniteElement& el, + const Mapping& 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 >(el, mapping, integrator.cell_quadrature, + integrator.cell_flags); + bdry_info.initialize >(el, mapping, integrator.bdry_quadrature, + integrator.face_flags); + face_info.initialize >(el, mapping, integrator.face_quadrature, + integrator.face_flags); + subface_info.initialize >(el, mapping, integrator.face_quadrature, + integrator.face_flags); + neighbor_info.initialize >(el, mapping, integrator.face_quadrature, + integrator.ngbr_flags); + } + + + template + template + void + IntegrationInfoBox::initialize( + const WORKER& integrator, + const FiniteElement& el, + const Mapping& mapping, + const NamedData >& data) + { + cell_info.initialize_data(data); + bdry_info.initialize_data(data); + face_info.initialize_data(data); + subface_info.initialize_data(data); + neighbor_info.initialize_data(data); + + initialize(integrator, el, mapping); + } +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h new file mode 100644 index 0000000000..c92dbc15e5 --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h @@ -0,0 +1,192 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace MeshWorker +{ + template + DoFInfo::DoFInfo(const BlockInfo& info) + : block_info(&info, typeid(*this).name()) + {} + + + template + void + DoFInfo::get_indices(const typename DoFHandler::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;iblock_info->renumber(i)] = indices_org[i]; + } + } + + + template + void + DoFInfo::get_indices(const typename MGDoFHandler::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;iblock_info->renumber(i)] = indices_org[i]; + } + } + +//----------------------------------------------------------------------// + + template + IntegrationInfo::IntegrationInfo(const BlockInfo& block_info) + : + DoFInfo(block_info), + fevalv(0), + multigrid(false), + global_data(boost::shared_ptr >(new VectorDataBase)) + {} + + + template + template + void + IntegrationInfo::initialize( + const FiniteElement& el, + const Mapping& mapping, + const Quadrature& quadrature, + const UpdateFlags flags) + { + if (this->block_info == 0 || this->block_info->local().size() == 0) + { + fevalv.resize(1); + fevalv[0] = boost::shared_ptr ( + new FEVALUES (mapping, el, quadrature, flags)); + } + else + { + fevalv.resize(el.n_base_elements()); + for (unsigned int i=0;i ( + 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;ilocal_indices().size()); + // For all components + for (unsigned int j=0;jn_gradients()); + // For all selected finite + // element functions + for (unsigned int i=0;ilocal_indices().size()); + // For all components + for (unsigned int j=0;jn_hessians()); + // For all selected finite + // element functions + for (unsigned int i=0;ilocal_indices().size()); + // For all components + for (unsigned int j=0;j + void + IntegrationInfo::initialize_data( + const boost::shared_ptr > data) + { + global_data = data; + } + + + template + void + IntegrationInfo::clear() + { + fevalv.resize(0); + } + + + + template + void + IntegrationInfo::fill_local_data(bool split_fevalues) + { + if (split_fevalues) + { + unsigned int comp = 0; + // Loop over all blocks + for (unsigned int b=0;bblock_info->local().size();++b) + { + const unsigned int fe_no = this->block_info->base_element(b); + const FEValuesBase& 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& fe = this->fe(0); + const unsigned int n_comp = fe.get_fe().n_components(); + this->global_data->fill(values, gradients, hessians, fe, this->indices, + 0, n_comp, 0, this->indices.size()); + } + } +} + + +DEAL_II_NAMESPACE_CLOSE + diff --git a/deal.II/deal.II/include/numerics/mesh_worker_loop.h b/deal.II/deal.II/include/numerics/mesh_worker_loop.h new file mode 100644 index 0000000000..c10f4990a4 --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_loop.h @@ -0,0 +1,188 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__mesh_worker_loop_h +#define __deal2__mesh_worker_loop_h + +DEAL_II_NAMESPACE_OPEN + +template class TriaActiveIterator; + +namespace MeshWorker +{ + namespace internal + { +/** + * Find out if an iterator supports inactive cells. + */ + template + inline bool is_active_iterator(const DI&) + { + return false; + } + + template + inline bool is_active_iterator(const TriaActiveIterator&) + { + 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 + void loop(ITERATOR begin, ENDITERATOR end, + CELLINFO& cellinfo, FACEINFO& bdryinfo, + FACEINFO& faceinfo, FACEINFO& subfaceinfo, FACEINFO& ngbrinfo, + LOCALWORKER& localworker, + bool cells_first = true) + { + const_cast&>( + 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::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 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 + void integration_loop(ITERATOR begin, ENDITERATOR end, + IntegrationInfoBox& box, + LOCALWORKER& localworker, + bool cells_first = true) + { + loop(begin, end, box.cell_info, box.bdry_info, + box.face_info, box.subface_info, box.neighbor_info, + localworker, cells_first); + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h new file mode 100644 index 0000000000..45397d80ee --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h @@ -0,0 +1,375 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__mesh_worker_vector_selector_h +#define __deal2__mesh_worker_vector_selector_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template 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 + void initialize(const NamedData&); + + /** + * 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 + void print (STREAM& s, const NamedData& 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 + 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 > >& values, + std::vector > > >& gradients, + std::vector > > >& hessians, + const FEValuesBase& fe, + const std::vector& 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 VectorData : + public VectorDataBase + { + public: + void initialize(const NamedData >&); + + virtual void fill( + std::vector > >& values, + std::vector > > >& gradients, + std::vector > > >& hessians, + const FEValuesBase& fe, + const std::vector& index, + unsigned int component, + unsigned int n_comp, + unsigned int start, + unsigned int size) const; + private: + SmartPointer > > 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 + inline void + VectorSelector::initialize(const NamedData& 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 + inline void + VectorSelector::print(STREAM& s, const NamedData& v) const + { + s << "values: "; + for (unsigned int i=0;i; +} diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h new file mode 100644 index 0000000000..b6cf863947 --- /dev/null +++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h @@ -0,0 +1,86 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + template + VectorDataBase::~VectorDataBase() + {} + + template + void + VectorDataBase::fill( + std::vector > >&, + std::vector > > >&, + std::vector > > >&, + const FEValuesBase&, + const std::vector&, + unsigned int, + unsigned int, + unsigned int, + unsigned int) const + {} +//----------------------------------------------------------------------// + + template + void + VectorData::initialize(const NamedData >& d) + { + data = &d; + VectorSelector::initialize(d); + } + + + template + void + VectorData::fill( + std::vector > >& values, + std::vector > > >& gradients, + std::vector > > >& hessians, + const FEValuesBase& fe, + const std::vector& index, + unsigned int component, + unsigned int n_comp, + unsigned int start, + unsigned int size) const + { + for (unsigned int i=0;in_values();++i) + { + const VECTOR& src = *(*data)(this->value_index(i)); + VectorSlice > > dst(values[i], component, n_comp); + fe.get_function_values(src, make_slice(index, start, size), dst, true); + } + + for (unsigned int i=0;in_gradients();++i) + { + const VECTOR& src = *(*data)(this->value_index(i)); + VectorSlice > > > dst(gradients[i], component, n_comp); + fe.get_function_gradients(src, make_slice(index, start, size), dst, true); + } + + for (unsigned int i=0;in_hessians();++i) + { + const VECTOR& src = *(*data)(this->value_index(i)); + VectorSlice > > > dst(hessians[i], component, n_comp); + fe.get_function_hessians(src, make_slice(index, start, size), dst, true); + } + } +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/fe_field_function.inst.in b/deal.II/deal.II/source/numerics/fe_field_function.inst.in index 4082dd26a4..54a65c46a7 100644 --- a/deal.II/deal.II/source/numerics/fe_field_function.inst.in +++ b/deal.II/deal.II/source/numerics/fe_field_function.inst.in @@ -1,8 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ -// Version: $Name$ // -// Copyright (C) 2007, 2008 by the deal.II authors +// Copyright (C) 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/source/numerics/mesh_worker.cc b/deal.II/deal.II/source/numerics/mesh_worker.cc new file mode 100644 index 0000000000..465eabd06a --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker.cc @@ -0,0 +1,116 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006, 2007, 2008, 2009 by Guido Kanschat +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +using namespace dealii; +using namespace MeshWorker; + +template +IntegrationWorker::IntegrationWorker () +{ + cell_flags = update_JxW_values; + bdry_flags = UpdateFlags(update_JxW_values | update_normal_vectors); + face_flags = bdry_flags; + ngbr_flags = update_default; +} + + +template +void +IntegrationWorker::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 +void +IntegrationWorker::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 +void +IntegrationWorker::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 +void +IntegrationWorker::initialize_gauss_quadrature( + unsigned int cp, + unsigned int bp, + unsigned int fp) +{ + cell_quadrature = QGauss(cp); + bdry_quadrature = QGauss(bp); + face_quadrature = QGauss(fp); + +} + + +template class IntegrationWorker; + diff --git a/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc b/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc new file mode 100644 index 0000000000..c57593172c --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc @@ -0,0 +1,130 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace MeshWorker +{ + namespace Assembler + { +//----------------------------------------------------------------------// + + template + inline void + Functional::initialize(unsigned int n) + { + results.resize(n); + } + + + template + template + inline void + Functional::assemble(const DoFInfo& info) + { + for (unsigned int i=0;i + template + inline void + Functional::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;i + inline number + Functional::operator() (unsigned int i) const + { + AssertIndexRange(i, results.size()); + return results[i]; + } + +//----------------------------------------------------------------------// + + template + inline void + CellsAndFaces::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 + template + inline void + CellsAndFaces::assemble(const DoFInfo& info) + { + for (unsigned int i=0;iuser_index()) += info.J[i]; + else + results.vector(0).block(i)(info.cell->user_index()) += info.J[i]; + } + } + + + template + template + inline void + CellsAndFaces::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;iuser_index()) += J; + if (info2.face != info1.face) + results.vector(1).block(i)(info2.face->user_index()) += J; + } + else + { + results.vector(0).block(i)(info1.cell->user_index()) += info1.J[i]; + results.vector(0).block(i)(info2.cell->user_index()) += info2.J[i]; + } + } + } + + +//----------------------------------------------------------------------// + + + + } +} + + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/mesh_worker_info.cc b/deal.II/deal.II/source/numerics/mesh_worker_info.cc new file mode 100644 index 0000000000..d7c6e7729b --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker_info.cc @@ -0,0 +1,42 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + template class LocalResults; + template class DoFInfo; + +#include "mesh_worker_info.inst" + + template void IntegrationInfo > + ::initialize >( + const FiniteElement&, const Mapping&, + const Quadrature::integral_dimension>&, const UpdateFlags); + template void IntegrationInfo > + ::initialize >( + const FiniteElement&, const Mapping&, + const Quadrature::integral_dimension>&, const UpdateFlags); + template void IntegrationInfo > + ::initialize >( + const FiniteElement&, const Mapping&, + const Quadrature::integral_dimension>&, const UpdateFlags); +} + + +DEAL_II_NAMESPACE_CLOSE + diff --git a/deal.II/deal.II/source/numerics/mesh_worker_info.inst.in b/deal.II/deal.II/source/numerics/mesh_worker_info.inst.in new file mode 100644 index 0000000000..8e123d5444 --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker_info.inst.in @@ -0,0 +1,34 @@ +//--------------------------------------------------------------------------- +// $Id: fe_field_function.inst.in 19046 2009-07-08 19:30:23Z bangerth $ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +for (FEV : FEVALUES_BASES) +{ + template class IntegrationInfo; +} + +for (FEV : FEVALUES_BASES) +{ +// for (TYPE : DERIVATIVE_TENSORS) +// { +// template void IntegrationInfo::fill_local_data( +// std::vector > >&, bool) const; +// template void IntegrationInfo::fill_local_data( +// std::vector > >&, VectorSelector&, bool) const; +// } +// for (VECTOR : SERIAL_VECTORS) +// { +// } + } + + + diff --git a/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc new file mode 100644 index 0000000000..9296026559 --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.cc @@ -0,0 +1,28 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + template class VectorDataBase; + +#include "mesh_worker_vector_selector.inst" +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in new file mode 100644 index 0000000000..17ba91a9e5 --- /dev/null +++ b/deal.II/deal.II/source/numerics/mesh_worker_vector_selector.inst.in @@ -0,0 +1,17 @@ +//--------------------------------------------------------------------------- +// $Id: fe_field_function.inst.in 19046 2009-07-08 19:30:23Z bangerth $ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +for (VECTOR : SERIAL_VECTORS) +{ + template class VectorData; +} diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 2e79616d83..2feb006baf 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -38,7 +38,8 @@ // Here come the new include files // for using the MeshWorker framework - +#include +#include #include #include @@ -178,280 +179,219 @@ void BoundaryValues::value_list(const std::vector > &points, } - // @sect3{Class: DGTransportEquation} + // @sect3{Integrating cell and face matrices} + // + // In order to use the MeshWorker + // framework, we need to define a + // worker object, which serves three + // purposes. First, it performs the + // local integration on each mesh + // cell. Second, it controls the loop + // over cells and faces and the data + // being generated. This part is + // easily implemented here, since the + // default structure can be + // used. Finally, through its base + // class from the + // MeshWorker::Assembler namespace, + // this object determins, how locally + // produced data is assembled into + // the global system. // - // Next we define the - // equation-dependent and - // DG-method-dependent class - // DGTransportEquation. Its - // member functions were already - // mentioned in the Introduction and - // will be explained - // below. Furthermore it includes - // objects of the previously defined - // Beta, RHS and - // BoundaryValues 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 -class DGTransportEquation +class DGIntegrator : public Subscriptor { public: - DGTransportEquation(); - - void assemble_cell_term(const FEValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const; - - void assemble_boundary_term(const FEFaceValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const; - - void assemble_face_term1(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix) const; - - void assemble_face_term2(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix, - FullMatrix &ui_ve_matrix, - FullMatrix &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::CellInfo CellInfo; + typedef typename MeshWorker::IntegrationWorker::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 beta_function; const RHS rhs_function; const BoundaryValues 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 -DGTransportEquation::DGTransportEquation () - : - beta_function (), - rhs_function (), - boundary_function () -{} - - - // @sect4{Function: assemble_cell_term} - // - // The assemble_cell_term - // function assembles the cell terms - // of the discretization. - // ui_vi_matrix is a cell matrix, - // i.e. for a DG method of degree 1, - // it is of size 4 times 4, and - // cell_vector is of size 4. - // When this function is invoked, - // fe_v is already reinit'ed with the - // current cell before and includes - // all shape values needed. -template -void DGTransportEquation::assemble_cell_term( - const FEValues &fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const +void DGIntegrator::cell(CellInfo& info) const { - // First we ask fe_v 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& fe_v = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + Vector& local_vector = info.R[0].block(0); + + // With these objects, we continue + // local integration like in step 12. const std::vector &JxW = fe_v.get_JxW_values (); - - // Then the flow field beta and the - // rhs_function are evaluated at - // the quadrature points, + std::vector > beta (fe_v.n_quadrature_points); std::vector 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; pointassemble_boundary_term - // function assembles the face terms - // at boundary faces. When this - // function is invoked, fe_v 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 -void DGTransportEquation::assemble_boundary_term( - const FEFaceValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const +void DGIntegrator::bdry(FaceInfo& info) const { - // Again, as in the previous - // function, we ask the - // FEValues object for the - // quadrature weights + const FEFaceValuesBase& fe_v = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + Vector& local_vector = info.R[0].block(0); + const std::vector &JxW = fe_v.get_JxW_values (); - // but here also for the normals. const std::vector > &normals = fe_v.get_normal_vectors (); - - // We evaluate the flow field - // and the boundary values at the - // quadrature points. + std::vector > beta (fe_v.n_quadrature_points); std::vector 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; point0) for (unsigned int i=0; iassemble_face_term1 - // 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, - // fe_v and fe_v_neighbor 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 - // ui_vi_matrix this function - // gets a new argument - // ue_vi_matrix, that stores - // contributions to the system matrix - // that are based on exterior values - // of $u$ and interior values of - // $v$. Here we note that ue is - // the short notation for u - // exterior and represents $u_h^-$, - // see the introduction. -template -void DGTransportEquation::assemble_face_term1( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &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 &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); - - // and we evaluate the flow field - // at the quadrature points. - std::vector > 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; point0) - for (unsigned int i=0; iassemble_face_term2 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 - // ui_ve_matrix and - // ue_ve_matrix 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 -void DGTransportEquation::assemble_face_term2( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix, - FullMatrix &ui_ve_matrix, - FullMatrix &ue_ve_matrix) const +void DGIntegrator::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& fe_v = info1.fe(); + + // For additional shape functions, + // we have to ask the neighbors + // FEFaceValuseBase. + const FEFaceValuesBase& 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& u1_v1_matrix = info1.M1[0].matrix; + FullMatrix& u2_v1_matrix = info1.M2[0].matrix; + FullMatrix& u1_v2_matrix = info2.M2[0].matrix; + FullMatrix& 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 &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); @@ -468,7 +408,7 @@ void DGTransportEquation::assemble_face_term2( // seen. for (unsigned int i=0; i::assemble_face_term2( // \kappa_+}$, for (unsigned int k=0; k::assemble_face_term2( // seen, too. for (unsigned int i=0; i::assemble_face_term2( // \kappa_-}$. for (unsigned int k=0; k &solution); void refine_grid (); void output_results (const unsigned int cycle) const; @@ -556,13 +495,6 @@ class DGMethod SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - // We define the quadrature - // formulae for the cell and the - // face terms of the - // discretization. - const QGauss quadrature; - const QGauss face_quadrature; - // And there are two solution // vectors, that store the // solutions to the problems @@ -570,30 +502,19 @@ class DGMethod // different assembling routines // assemble_system1 and // assemble_system2; - Vector solution1; - Vector solution2; - Vector right_hand_side; - - // Finally this class includes an - // object of the - // DGTransportEquations class - // described above. - const DGTransportEquation dg; + Vector solution; + Vector right_hand_side; }; template DGMethod::DGMethod () : - mapping (), // Change here for DG // methods of // different degrees. fe (1), - dof_handler (triangulation), - quadrature (4), - face_quadrature (4), - dg () + dof_handler (triangulation) {} @@ -635,723 +556,97 @@ void DGMethod::setup_system () system_matrix.reinit (sparsity_pattern); - solution1.reinit (dof_handler.n_dofs()); - solution2.reinit (dof_handler.n_dofs()); + solution.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); } - - // @sect4{Function: assemble_system1} - // - // We proceed with the - // assemble_system1 function that - // implements the DG discretization - // in its first version. This - // function repeatedly calls the - // assemble_cell_term, - // assemble_boundary_term and - // assemble_face_term1 functions - // of the DGTransportEquation - // object. The - // assemble_boundary_term covers - // the first case mentioned in the - // introduction. - // - // 1. face is at boundary - // - // This function takes a - // FEFaceValues object as - // argument. In contrast to that - // assemble_face_term1 - // takes two FEFaceValuesBase - // 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 FEFaceValues - // or of class FESubfaceValues - // (both derived from - // FEFaceValuesBase) according to - // the remaining cases mentioned - // in the introduction: - // - // 2. neighboring cell is finer - // (current cell: FESubfaceValues, - // neighboring cell: FEFaceValues); - // - // 3. neighboring cell is of the same - // refinement level (both, current - // and neighboring cell: - // FEFaceValues); - // - // 4. neighboring cell is coarser - // (current cell: FEFaceValues, - // neighboring cell: - // FESubfaceValues). - // - // 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 -void DGMethod::assemble_system1 () -{ - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - std::vector dofs_neighbor (dofs_per_cell); - - // First we create the - // update_flags for the - // FEValues and the - // FEFaceValues 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 FEValues - // object. Note, that since version - // 3.2.0 of deal.II the constructor - // of this class takes a - // Mapping object as first - // argument. Although the - // constructor without Mapping - // 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 MappingQ1 mapping) - // and makes it easier to change - // the mapping object later. - FEValues fe_v ( - mapping, fe, quadrature, update_flags); - - // Similarly we create the - // FEFaceValues and - // FESubfaceValues 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 fe_v_face ( - mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( - mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - FESubfaceValues 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 vi (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 ui_vi_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); - - Vector cell_vector (dofs_per_cell); - - // Furthermore we need some cell - // iterators. - typename DoFHandler::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 - // assemble_face_term1 - // function contributions to - // the cell matrices and the - // cell vector are only - // ADDED. Therefore on each - // cell we need to reset the - // ui_vi_matrix and - // cell_vector to zero, - // before assembling the cell terms. - ui_vi_matrix = 0; - cell_vector = 0; - - // Now we reinit the FEValues - // object for the current cell - fe_v.reinit (cell); - - // and call the function - // that assembles the cell - // terms. The first argument is - // the FEValues 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::faces_per_cell; ++face_no) - { - // First we set the face - // iterator - typename DoFHandler::face_iterator face=cell->face(face_no); - - // and clear the - // ue_vi_matrix 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 - // FEFaceValues - // 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::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 dim==1. - 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 - // cell. - const unsigned int neighbor2= - cell->neighbor_of_neighbor(face_no); - - - // We loop over - // subfaces - for (unsigned int subface_no=0; - subface_non_children(); ++subface_no) - { - // and set the - // cell - // iterator - // neighbor_child - // to the cell - // placed - // `behind' the - // current - // subface. - typename DoFHandler::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 - // ue_vi_matrix - // on each - // subface - // because on - // each subface - // the un - // belong to - // different - // neighboring - // cells. - ue_vi_matrix = 0; - - // As already - // mentioned - // above for - // the current - // case (case - // 2) we employ - // the - // FESubfaceValues - // 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 - // ue_vi_matrix - // to the - // system_matrix - for (unsigned int i=0; iif - // (face-@>has_children()) - } - 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 - // FEFaceValues - // 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 if - // (neighbor-@>level() - // == - // cell-@>level()) - } - 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 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 - // FEFaceValues - // 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 - // neighbor_child - // cell, - neighbor->get_dof_indices (dofs_neighbor); - - // and distribute the - // ue_vi_matrix. - for (unsigned int i=0; iface not at boundary: - } - // End of loop over all faces: - } - - // Finally we distribute the - // ui_vi_matrix - for (unsigned int i=0; iassemble_system2 function that - // implements the DG discretization - // in its second version. This - // function is very similar to the - // assemble_system1 - // function. Therefore, here we only - // discuss the differences between - // the two functions. This function - // repeatedly calls the - // assemble_face_term2 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 -void DGMethod::assemble_system2 () +void DGMethod::assemble_system () { - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - std::vector 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 - // fe_v_face_neighbor as case 4 - // does not occur. - FEValues fe_v ( - mapping, fe, quadrature, update_flags); - FEFaceValues fe_v_face ( - mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( - mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - - - FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix 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 ve (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 ui_ve_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix 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 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, Vector >, DGIntegrator > + integrator(dg); - Vector cell_vector (dofs_per_cell); - - // The following lines are roughly - // the same as in the previous - // function. - typename DoFHandler::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::faces_per_cell; ++face_no) - { - typename DoFHandler::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::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_non_children(); ++subface_no) - { - typename DoFHandler::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; ilevel() == 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 info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + + // Finally, the integration loop + // over all active cells + // (determined by the first + // argument, which is an active iterator). + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); } @@ -1453,7 +748,7 @@ void DGMethod::refine_grid () // are computed DerivativeApproximation::approximate_gradient (mapping, dof_handler, - solution2, + solution, gradient_indicator); // and they are cell-wise scaled by @@ -1509,7 +804,7 @@ void DGMethod::output_results (const unsigned int cycle) const DataOut data_out; data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution2, "u"); + data_out.add_data_vector (solution, "u"); data_out.build_patches (); @@ -1557,53 +852,9 @@ void DGMethod::run () << dof_handler.n_dofs() << std::endl; - // The constructor of the Timer - // class automatically starts - // the time measurement. - Timer assemble_timer; - // First assembling routine. - assemble_system1 (); - // The operator () accesses the - // current time without - // disturbing the time - // measurement. - std::cout << "Time of assemble_system1: " - << assemble_timer() - << std::endl; - solve (solution1); - - // As preparation for the - // second assembling routine we - // reinit the system matrix, the - // right hand side vector and - // the Timer object. - system_matrix = 0; - right_hand_side = 0; - assemble_timer.reset(); - - // We start the Timer, - assemble_timer.start(); - // call the second assembling routine - assemble_system2 (); - // and access the current time. - std::cout << "Time of assemble_system2: " - << assemble_timer() - << std::endl; - solve (solution2); - - // To make sure that both - // versions of the DG method - // yield the same - // discretization and hence the - // same solution we check the - // two solutions for equality. - solution1-=solution2; - const double difference=solution1.linfty_norm(); - if (difference>1e-13) - std::cout << "solution1 and solution2 differ!!" << std::endl; - else - std::cout << "solution1 and solution2 coincide." << std::endl; - + assemble_system (); + solve (solution); + // Finally we perform the // output. output_results (cycle); -- 2.39.5