From 4f61f93dc7aa7d10d9bc6b4a5f2186b290485b58 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 13 Dec 2015 11:27:09 -0600 Subject: [PATCH] Update the documentation of the DoFTools::make_*_sparsity_pattern() functions. --- include/deal.II/dofs/dof_tools.h | 432 +++++++++++++++++++------------ 1 file changed, 270 insertions(+), 162 deletions(-) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index 975d1eac13..560dd62bdd 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -203,7 +203,7 @@ namespace DoFTools /** * Two components couple only if their shape functions are both nonzero on * a given face. This flag is only used when computing integrals over - * faces of cells. + * faces of cells, e.g., in DoFTools::make_flux_sparsity_pattern(). */ nonzero }; @@ -342,80 +342,165 @@ namespace DoFTools */ /** - * Locate non-zero entries of the system matrix. + * Compute which entries of a matrix built on the given + * @p dof_handler may possibly be nonzero, and create a sparsity + * pattern object that represents these nonzero locations. * - * This function computes the possible positions of non-zero entries in the - * global system matrix. We assume that a certain finite element basis + * This function computes the possible positions of non-zero entries + * in the global system matrix by simulating which entries + * one would write to during the actual assembly of a matrix. For + * this, the function assumes that each finite element basis * function is non-zero on a cell only if its degree of freedom is - * associated with the interior, a face, an edge or a vertex of this cell. - * As a result, the matrix entry between two basis functions can be non-zero - * only if they correspond to degrees of freedom of at least one common - * cell. Therefore, @p make_sparsity_pattern just loops over all cells and - * enters all couplings local to that cell. As the generation of the - * sparsity pattern is irrespective of the equation which is solved later - * on, the resulting sparsity pattern is symmetric. - * - * Remember using SparsityPattern::compress() after generating the pattern. - * - * The actual type of the sparsity pattern may be SparsityPattern, - * DynamicSparsityPattern, BlockSparsityPattern, - * BlockDynamicSparsityPattern, or any other class that satisfies similar - * requirements. It is assumed that the size of the sparsity pattern matches - * the number of degrees of freedom and that enough unused nonzero entries - * are left to fill the sparsity pattern. The nonzero entries generated by - * this function are overlaid to possible previous content of the object, - * that is previously added entries are not deleted. - * - * Since this process is purely local, the sparsity pattern does not provide - * for entries introduced by the elimination of hanging nodes. They have to - * be taken care of by a call to ConstraintMatrix::condense() afterwards. - * - * Alternatively, the constraints on degrees of freedom can already be taken - * into account at the time of creating the sparsity pattern. For this, pass - * the ConstraintMatrix object as the third argument to the current - * function. No call to ConstraintMatrix::condense() is then necessary. This - * process is explained in step-27. - * - * In case the constraints are already taken care of in this function, it is - * possible to neglect off-diagonal entries in the sparsity pattern. When - * using ConstraintMatrix::distribute_local_to_global during assembling, no - * entries will ever be written into these matrix position, so that one can - * save some computing time in matrix-vector products by not even creating - * these elements. In that case, the variable keep_constrained_dofs - * needs to be set to false. - * - * If the @p subdomain_id parameter is given, the sparsity pattern is built - * only on cells that have a subdomain_id equal to the given argument. This - * is useful in parallel contexts where the matrix and sparsity pattern (for - * example a TrilinosWrappers::SparsityPattern) may be distributed and not - * every MPI process needs to build the entire sparsity pattern; in that - * case, it is sufficient if every process only builds that part of the - * sparsity pattern that corresponds to the subdomain_id for which it is - * responsible. This feature is used in step-32. + * associated with the interior, a face, an edge or a vertex of this + * cell. As a result, a matrix entry $A_{ij}$ that is computed from + * two basis functions $\varphi_i$ and $\varphi_j$ with (global) + * indices $i$ and $j$ (for example, using a bilinear form + * $A_{ij}=a(\varphi_i,\varphi_j)$) can be non-zero only if these + * shape functions correspond to degrees of freedom that are defined + * on at least one common cell. Therefore, this function just loops + * over all cells, figures out the global indices of all degrees of + * freedom, and presumes that all matrix entries that couple any of + * these indices will result in a nonzero matrix entry. These will + * then be added to the sparsity pattern. As this process of + * generating the sparsity pattern does not take into account the + * equation to be solved later on, the resulting sparsity pattern is + * symmetric. + * + * This algorithm makes no distinction between shape functions on + * each cell, i.e., it simply couples all degrees of freedom on a + * cell with all other degrees of freedom on a cell. This is often + * the case, and always a safe assumption. However, if you know + * something about the structure of your operator and that it does + * not couple certain shape functions with certain test functions, + * then you can get a sparser sparsity pattern by calling a variant + * of the current function described below that allows to specify + * which vector components couple with which other vector + * components. + * + * The method described above lives on the assumption that coupling + * between degrees of freedom only happens if shape functions + * overlap on at least one cell. This is the case with most usual + * finite element formulations involving conforming + * elements. However, for formulations such as the Discontinuous + * Galerkin finite element method, the bilinear form contains terms + * on interfaces between cells that couple shape functions that live + * on one cell with shape functions that live on a neighboring + * cell. The current function would not see these couplings, and + * would consequently not allocate entries in the sparsity + * pattern. You would then get into trouble during matrix assembly + * because you try to write into matrix entries for which no space + * has been allocated in the sparsity pattern. This can be avoided + * by calling the DoFTools::make_flux_sparsity_pattern() function + * instead, which takes into account coupling between degrees of + * freedom on neighboring cells. + * + * There are other situations where bilinear forms contain non-local + * terms, for example in treating integral equations. These require + * different methods for building the sparsity patterns that depend + * on the exact formulation of the problem. You will have to do this + * yourself then. + * + * @param[in] dof_handler The DoFHandler or hp::DoFHandler object + * that describes which degrees of freedom live on which cells. + * + * @param[out] sparsity_pattern The sparsity pattern to be filled with + * entries. + * + * @param[in] constraints The process for generating entries + * described above is purely local to each cell. Consequently, the + * sparsity pattern does not provide for matrix entries that will + * only be written into during the elimination of hanging nodes or + * other constraints. They have to be taken care of by a + * subsequent call to ConstraintMatrix::condense(). + * Alternatively, the constraints on degrees of freedom can + * already be taken into account at the time of creating the + * sparsity pattern. For this, pass the ConstraintMatrix object as + * the third argument to the current function. No call to + * ConstraintMatrix::condense() is then necessary. This process is + * explained in step-6, step-27, and other tutorial programs. + * + * @param[in] keep_constrained_dofs In case the constraints are + * already taken care of in this function by passing in a + * ConstraintMatrix object, it is possible to abandon some + * off-diagonal entries in the sparsity pattern if these entries + * will also not be written into during the actual assembly of the + * matrix this sparsity pattern later serves. Specifically, when + * using an assembly method that uses + * ConstraintMatrix::distribute_local_to_global(), no entries will + * ever be written into those matrix rows or columns that + * correspond to constrained degrees of freedom. In such cases, + * you can set the argument @p keep_constrained_dofs to @p false + * to avoid allocating these entries in the sparsity pattern. + * + * @param[in] subdomain_id If specified, the sparsity pattern is + * built only on cells that have a subdomain_id equal to the given + * argument. This is useful in parallel contexts where the matrix + * and sparsity pattern (for example a + * TrilinosWrappers::SparsityPattern) may be distributed and not + * every MPI process needs to build the entire sparsity pattern; + * in that case, it is sufficient if every process only builds + * that part of the sparsity pattern that corresponds to the + * subdomain_id for which it is responsible. This feature is used + * in step-32. (This argument is not usually needed for objects of + * type parallel::distributed::Triangulation because the current + * function only loops over locally owned cells anyway; thus, this + * argument typically only makes sense if you want to use the + * subdomain_id for anything other than indicating which processor + * owns a cell, for example which geometric component of the + * domain a cell belongs to.) + * + * @note The actual type of the sparsity pattern may be + * SparsityPattern, DynamicSparsityPattern, BlockSparsityPattern, + * BlockDynamicSparsityPattern, or any other class that satisfies + * similar requirements. It is assumed that the size of the + * sparsity pattern matches the number of degrees of freedom and + * that enough unused nonzero entries are left to fill the + * sparsity pattern if the sparsity pattern is of "static" kind + * (see @ref Sparsity for more information on what this + * means). The nonzero entries generated by this function are + * added to possible previous content of the object, i.e., + * previously added entries are not removed. + * + * @note If the sparsity pattern is represented by an object of type + * SparsityPattern (as opposed to, for example, + * DynamicSparsityPattern), you need to remember using + * SparsityPattern::compress() after generating the pattern. * * @ingroup constraints */ template void - make_sparsity_pattern (const DoFHandlerType &dof, + make_sparsity_pattern (const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const ConstraintMatrix &constraints = ConstraintMatrix(), const bool keep_constrained_dofs = true, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id); /** - * Locate non-zero entries for vector valued finite elements. This function - * does mostly the same as the previous make_sparsity_pattern(), but it is - * specialized for vector finite elements and allows to specify which - * variables couple in which equation. For example, if wanted to solve the - * Stokes equations, + * Compute which entries of a matrix built on the given + * @p dof_handler may possibly be nonzero, and create a sparsity + * pattern object that represents these nonzero locations. + * + * This function is a simple variation on the previous + * make_sparsity_pattern() function (see there for a description of + * all of the common arguments), but it provides functionality for + * vector finite elements that allows to be more specific about + * which variables couple in which equation. + * + * For example, if you wanted to solve the Stokes equations, * * @f{align*} -\Delta \mathbf u + \nabla p &= 0,\\ \text{div}\ u &= 0 @f} * - * in two space dimensions, using stable Q2/Q1 mixed elements (using the - * FESystem class), then you don't want all degrees of freedom to couple in - * each equation. You rather may want to give the following pattern of - * couplings: + * in two space dimensions, using stable Q2/Q1 mixed elements (using + * the FESystem class), then you don't want all degrees of freedom + * to couple in each equation. More specifically, in the first + * equation, only $u_x$ and $p$ appear; in the second equation, only + * $u_y$ and $p$ appear; and in the third equation, only $u_x$ and + * $u_y$ appear. (Note that this discussion only talks about vector + * components of the solution variable and the different equation, + * and has nothing to do with degrees of freedom, or in fact with + * any kind of discretization.) We can describe this by the + * following pattern of "couplings": * * @f[ * \left[ @@ -427,11 +512,18 @@ namespace DoFTools * \right] * @f] * - * where "1" indicates that two variables (i.e. components of the FESystem) - * couple in the respective equation, and a "0" means no coupling, in which - * case it is not necessary to allocate space in the matrix structure. - * Obviously, the mask refers to components of the composed FESystem, rather - * than to the degrees of freedom contained in there. + * where "1" indicates that two variables (i.e., vector components + * of the FESystem) couple in the respective equation, and a "0" + * means no coupling. These zeros imply that upon discretization via + * a standard finite element formulation, we will not write entries + * into the matrix that, for example, couple pressure test functions + * with pressure shape functions (and similar for the other zeros + * above). It is then a waste to allocate memory for these entries + * in the matrix and the sparsity pattern, and you can avoid this by + * creating a mask such as the one above that describes this to the + * (current) function that computes the sparsity pattern. As stated + * above, the mask shown above refers to components of the composed + * FESystem, rather than to degrees of freedom or shape functions. * * This function is designed to accept a coupling pattern, like the one * shown above, through the @p couplings parameter, which contains values of @@ -440,53 +532,18 @@ namespace DoFTools * coupling pattern. If the couplings are symmetric, then so will be the * resulting sparsity pattern. * - * The actual type of the sparsity pattern may be SparsityPattern, - * DynamicSparsityPattern, BlockSparsityPattern, - * BlockDynamicSparsityPattern, BlockDynamicSetSparsityPattern, or any other - * class that satisfies similar requirements. - * * There is a complication if some or all of the shape functions of the * finite element in use are non-zero in more than one component (in deal.II - * speak: they are non-primitive). In this case, the coupling element + * speak: they are @ref GlossPrimitive "non-primitive finite elements"). In + * this case, the coupling element * corresponding to the first non-zero component is taken and additional * ones for this component are ignored. * - * @todo Not implemented for hp::DoFHandler. - * - * As mentioned before, the creation of the sparsity pattern is a purely - * local process and the sparsity pattern does not provide for entries - * introduced by the elimination of hanging nodes. They have to be taken - * care of by a call to ConstraintMatrix::condense() afterwards. - * - * Alternatively, the constraints on degrees of freedom can already be taken - * into account at the time of creating the sparsity pattern. For this, pass - * the ConstraintMatrix object as the third argument to the current - * function. No call to ConstraintMatrix::condense() is then necessary. This - * process is explained in - * @ref step_27 "step-27". - * - * In case the constraints are already taken care of in this function, it is - * possible to neglect off-diagonal entries in the sparsity pattern. When - * using ConstraintMatrix::distribute_local_to_global during assembling, no - * entries will ever be written into these matrix position, so that one can - * save some computing time in matrix-vector products by not even creating - * these elements. In that case, the variable keep_constrained_dofs - * needs to be set to false. - * - * If the @p subdomain_id parameter is given, the sparsity pattern is built - * only on cells that have a subdomain_id equal to the given argument. This - * is useful in parallel contexts where the matrix and sparsity pattern (for - * example a TrilinosWrappers::SparsityPattern) may be distributed and not - * every MPI process needs to build the entire sparsity pattern; in that - * case, it is sufficient if every process only builds that part of the - * sparsity pattern that corresponds to the subdomain_id for which it is - * responsible. This feature is used in step-32. - * * @ingroup constraints */ template void - make_sparsity_pattern (const DoFHandlerType &dof, + make_sparsity_pattern (const DoFHandlerType &dof_handler, const Table<2, Coupling> &coupling, SparsityPatternType &sparsity_pattern, const ConstraintMatrix &constraints = ConstraintMatrix(), @@ -520,52 +577,51 @@ namespace DoFTools SparsityPatternType &sparsity); /** - * Create the sparsity pattern for boundary matrices. See the general - * documentation of this class for more information. - * - * The actual type of the sparsity pattern may be SparsityPattern, - * DynamicSparsityPattern, BlockSparsityPattern, - * BlockDynamicSparsityPattern, BlockDynamicSetSparsityPattern, or any other - * class that satisfies similar requirements. It is assumed that the size of - * the sparsity pattern is already correct. - */ - template - void - make_boundary_sparsity_pattern (const DoFHandlerType &dof, - const std::vector &dof_to_boundary_mapping, - SparsityPatternType &sparsity_pattern); - - /** - * Write the sparsity structure of the matrix composed of the basis - * functions on the boundary into the matrix structure. In contrast to the - * previous function, only those parts of the boundary are considered of - * which the boundary indicator is listed in the set of numbers passed to - * this function. + * Compute which entries of a matrix built on the given @p + * dof_handler may possibly be nonzero, and create a sparsity + * pattern object that represents these nonzero locations. This + * function is a variation of the make_sparsity_pattern() functions + * above in that it assumes that the bilinear form you want to use + * to generate the matrix also contains terms that integrate over + * the faces between cells (i.e., it contains "fluxes" + * between cells, explaining the name of the function). + * + * This function is useful for Discontinuous Galerkin methods where + * the standard make_sparsity_pattern() function would only create + * nonzero entries for all degrees of freedom on one cell coupling + * to all other degrees of freedom on the same cell; however, in DG + * methods, all or some degrees of freedom on each cell also couple + * to the degrees of freedom on other cells connected to the current + * one by a common face. The current function also creates the + * nonzero entries in the matrix resulting from these additional + * couplings. In other words, this function computes a strict + * super-set of nonzero entries compared to the work done by + * make_sparsity_pattern(). + * + * @param[in] dof_handler The DoFHandler or hp::DoFHandler object + * that describes which degrees of freedom live on which cells. + * + * @param[out] sparsity_pattern The sparsity pattern to be filled with + * entries. + * + * @note The actual type of the sparsity pattern may be + * SparsityPattern, DynamicSparsityPattern, BlockSparsityPattern, + * BlockDynamicSparsityPattern, or any other class that satisfies + * similar requirements. It is assumed that the size of the + * sparsity pattern matches the number of degrees of freedom and + * that enough unused nonzero entries are left to fill the + * sparsity pattern if the sparsity pattern is of "static" kind + * (see @ref Sparsity for more information on what this + * means). The nonzero entries generated by this function are + * added to possible previous content of the object, i.e., + * previously added entries are not removed. + * + * @note If the sparsity pattern is represented by an object of type + * SparsityPattern (as opposed to, for example, + * DynamicSparsityPattern), you need to remember using + * SparsityPattern::compress() after generating the pattern. * - * In fact, rather than a @p set of boundary indicators, a @p map needs to - * be passed, since most of the functions handling with boundary indicators - * take a mapping of boundary indicators and the respective boundary - * functions. The boundary function, however, is ignored in this function. - * If you have no functions at hand, but only the boundary indicators, set - * the function pointers to null pointers. - * - * For the type of the sparsity pattern, the same holds as said above. - */ - template - void - make_boundary_sparsity_pattern - (const DoFHandlerType &dof, - const typename FunctionMap::type &boundary_ids, - const std::vector &dof_to_boundary_mapping, - SparsityPatternType &sparsity); - - /** - * Generate sparsity pattern for fluxes, i.e. formulations of the discrete - * problem with discontinuous elements which couple across faces of cells. - * This is a replacement of the function @p make_sparsity_pattern for - * discontinuous methods. Since the fluxes include couplings between - * neighboring elements, the normal couplings and these extra matrix entries - * are considered. + * @ingroup constraints */ template void @@ -573,9 +629,11 @@ namespace DoFTools SparsityPatternType &sparsity_pattern); /** - * This function does the same as the other with the same name, but it gets - * a ConstraintMatrix additionally. This is for the case where you have - * fluxes but constraints as well. + * This function does essentially the same as the other + * make_flux_sparsity_pattern() function but allows the + * specification of a number of additional arguments. These carry + * the same meaning as discussed in the first + * make_sparsity_pattern() function above. * * @ingroup constraints */ @@ -588,22 +646,71 @@ namespace DoFTools const types::subdomain_id subdomain_id = numbers::invalid_unsigned_int); /** - * This function does the same as the other with the same name, but it gets - * two additional coefficient matrices. A matrix entry will only be - * generated for two basis functions, if there is a non-zero entry linking - * their associated components in the coefficient matrix. + * This function does essentially the same as the other + * make_flux_sparsity_pattern() function but allows the + * specification of coupling matrices that state which components of + * the solution variable couple in each of the equations you are + * discretizing. This works in complete analogy as discussed in the + * second make_sparsity_pattern() function above. * - * There is one matrix for couplings in a cell and one for the couplings - * occurring in fluxes. + * In fact, this function takes two such masks, one describing which + * variables couple with each other in the cell integrals that make + * up your bilinear form, and which variables coupld with each other + * in the face integrals. If you passed masks consisting of only 1s + * to both of these, then you would get the same sparsity pattern as + * if you had called the first of the make_sparsity_pattern() + * functions above. By setting some of the entries of these masks to + * zeros, you can get a sparser sparsity pattern. * - * @todo Not implemented for hp::DoFHandler. + * @ingroup constraints */ template void make_flux_sparsity_pattern (const DoFHandlerType &dof, SparsityPatternType &sparsity, - const Table<2,Coupling> &int_mask, - const Table<2,Coupling> &flux_mask); + const Table<2,Coupling> &cell_integrals_mask, + const Table<2,Coupling> &face_integrals_mask); + + /** + * Create the sparsity pattern for boundary matrices. See the general + * documentation of this class for more information. + * + * The function does essentially what the other + * make_sparsity_pattern() functions do, but assumes that the + * bilinear form that is used to build the matrix does not consist + * of domain integrals, but only of integrals over the boundary of + * the domain. + */ + template + void + make_boundary_sparsity_pattern (const DoFHandlerType &dof, + const std::vector &dof_to_boundary_mapping, + SparsityPatternType &sparsity_pattern); + + /** + * This function is a variation of the previous + * make_boundary_sparsity_pattern() function in which we assume that + * the boundary integrals that will give rise to the matrix extends + * only over those parts of the boundary whose boundary indicators + * are listed in the @p boundary_ids argument to this function. + * + * This function could have been written by passing a @p set of + * boundary_id numbers. However, most of the functions throughout + * deal.II dealing with boundary indicators take a mapping of + * boundary indicators and the corresponding boundary function, + * i.e., a FunctionMap argument. Correspondingly, this function does + * the same, though the actual boundary function is ignored here. + * (Consequently, if you don't have any such boundary functions, + * just create a map with the boundary indicators you want and set + * the function pointers to null pointers). + */ + template + void + make_boundary_sparsity_pattern + (const DoFHandlerType &dof, + const typename FunctionMap::type &boundary_ids, + const std::vector &dof_to_boundary_mapping, + SparsityPatternType &sparsity); /** * @} @@ -1468,9 +1575,10 @@ namespace DoFTools /** - * For each DoF, return in the output array to which subdomain (as given by - * the cell->subdomain_id() function) it belongs. The output array - * is supposed to have the right size already when calling this function. + * For each degree of freedom, return in the output array to which + * subdomain (as given by the cell->subdomain_id() + * function) it belongs. The output array is supposed to have the + * right size already when calling this function. * * Note that degrees of freedom associated with faces, edges, and vertices * may be associated with multiple subdomains if they are sitting on -- 2.39.5