From bf9eac97cc49b165376dfc1d1363432d2dd784d6 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 28 Apr 2011 03:34:42 +0000 Subject: [PATCH] Make DoFTools a namespace, rather than a class with all public, static member functions. git-svn-id: https://svn.dealii.org/trunk@23656 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/dofs/dof_tools.h | 4265 +++++++------- deal.II/source/dofs/dof_tools.cc | 6872 +++++++++++----------- 2 files changed, 5561 insertions(+), 5576 deletions(-) diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 9c2fa7407b..c76523c422 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -171,2013 +171,2008 @@ template class Mapping; * * * @ingroup dofs - * @author Wolfgang Bangerth, Guido Kanschat and others, 1998 - 2008 + * @author Wolfgang Bangerth, Guido Kanschat and others */ -class DoFTools +namespace DoFTools { - public: - /** - * The flags used in tables by certain - * make_*_pattern functions to - * describe whether two components of the - * solution couple in the bilinear forms - * corresponding to cell or face - * terms. An example of using these flags - * is shown in the introduction of - * step-46. - * - * In the descriptions of the individual - * elements below, remember that these - * flags are used as elements of tables - * of size FiniteElement::n_components - * times FiniteElement::n_components - * where each element indicates whether - * two components do or do not couple. - */ - enum Coupling - { - /** - * Two components do not - * couple. - */ - none, - /** - * Two components do couple. - */ - always, - /** - * 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. - */ - nonzero - }; - - /** - * @name Auxiliary functions - * @{ - */ - /** - * Maximal number of degrees of - * freedom on a cell. - */ - template - static unsigned int - max_dofs_per_cell (const DoFHandler &dh); - - template - static unsigned int - max_dofs_per_cell (const hp::DoFHandler &dh); - - - /** - * Maximal number of degrees of - * freedom on a face. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - max_dofs_per_face (const DoFHandler &dh); - - /** - * Maximal number of degrees of - * freedom on a face. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - max_dofs_per_face (const hp::DoFHandler &dh); - - /** - * Maximal number of degrees of - * freedom on a vertex. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - max_dofs_per_vertex (const DoFHandler &dh); - - /** - * Maximal number of degrees of - * freedom on a vertex. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - max_dofs_per_vertex (const hp::DoFHandler &dh); - - /** - * Number of vector components in the - * finite element object used by this - * DoFHandler. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - n_components (const DoFHandler &dh); - - /** - * Number of vector components in the - * finite element object used by this - * DoFHandler. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static unsigned int - n_components (const hp::DoFHandler &dh); - - /** - * Find out whether the FiniteElement - * used by this DoFHandler is primitive - * or not. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static bool - fe_is_primitive (const DoFHandler &dh); - - /** - * Find out whether the FiniteElement - * used by this DoFHandler is primitive - * or not. - * - * This function exists for both non-hp - * and hp DoFHandlers, to allow for a - * uniform interface to query this - * property. - */ - template - static bool - fe_is_primitive (const hp::DoFHandler &dh); - - /** - * @} - */ - - /** - * @name Sparsity Pattern Generation - * @{ - */ - - /** - * Locate non-zero entries of the - * system matrix. - * - * This function computes the - * possible positions of non-zero - * entries in the global system - * matrix. We assume that a - * certain 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, - * CompressedSparsityPattern, - * BlockSparsityPattern, - * BlockCompressedSparsityPattern, - * BlockCompressedSetSparsityPattern, - * BlockCompressedSimpleSparsityPattern, - * 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. - * - * @ingroup constraints - */ - template - static - void - make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); - - /** - * Locate non-zero entries for - * vector valued finite elements. - * This function does mostly the - * same as the previous - * @p 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, - * @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: - * @f[ - * \begin{array}{ccc} - * 1 & 0 & 1 \\ - * 0 & 1 & 1 \\ - * 1 & 1 & 0 - * \end{array} - * @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. - * - * This function is designed to - * accept a coupling pattern, like the one - * shown above, through the - * @p couplings parameter, which - * contains values of type #Coupling. It - * builds the matrix structure - * just like the previous - * function, but does not create - * matrix elements if not - * specified by the 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, - * CompressedSparsityPattern, - * BlockSparsityPattern, - * BlockCompressedSparsityPattern, - * BlockCompressedSetSparsityPattern, - * 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 - * correspoding to the first - * non-zero component is taken - * and additional ones for this - * component are ignored. - * - * 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 - static - void - make_sparsity_pattern (const DH &dof, - const Table<2, Coupling> &coupling, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); - - /** - * @deprecated This is the old - * form of the previous - * function. It generates a table - * of DoFTools::Coupling values - * (where a true - * value in the mask is - * translated into a - * Coupling::always value in the - * table) and calls the function - * above. - */ - template - static - void - make_sparsity_pattern (const DH &dof, - const std::vector > &mask, - SparsityPattern &sparsity_pattern); - - /** - * Construct a sparsity pattern that - * allows coupling degrees of freedom on - * two different but related meshes. - * - * The idea is that if the two given - * DoFHandler objects correspond to two - * different meshes (and potentially to - * different finite elements used on - * these cells), but that if the two - * triangulations they are based on are - * derived from the same coarse mesh - * through hierarchical refinement, then - * one may set up a problem where one - * would like to test shape functions - * from one mesh against the shape - * functions from another mesh. In - * particular, this means that shape - * functions from a cell on the first - * mesh are tested against those on the - * second cell that are located on the - * corresponding cell; this - * correspondence is something that the - * IntergridMap class can determine. - * - * This function then constructs a - * sparsity pattern for which the degrees - * of freedom that represent the rows - * come from the first given DoFHandler, - * whereas the ones that correspond to - * columns come from the second - * DoFHandler. - */ - template - static - void - make_sparsity_pattern (const DH &dof_row, - const DH &dof_col, - SparsityPattern &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, - * CompressedSparsityPattern, - * BlockSparsityPattern, - * BlockCompressedSparsityPattern, - * BlockCompressedSetSparsityPattern, or - * any other class that satisfies similar - * requirements. It is assumed that the - * size of the sparsity pattern is - * already correct. - */ - template - static void - make_boundary_sparsity_pattern (const DH &dof, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &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. - * - * 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 - static void - make_boundary_sparsity_pattern (const DH &dof, - const typename FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &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. - */ - template - static void - make_flux_sparsity_pattern (const DH &dof_handler, - SparsityPattern &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. - * Not implemented for - * hp::DoFHandler. - * - * @ingroup constraints - */ - template - static void - make_flux_sparsity_pattern (const DH &dof_handler, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs = true, - const types::subdomain_id_t 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. - * - * There is one matrix for - * couplings in a cell and one - * for the couplings occuring in - * fluxes. - * - * Not implemented for - * hp::DoFHandler. - */ - template - static void - make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const Table<2,Coupling> &int_mask, - const Table<2,Coupling> &flux_mask); - - //@} - /** - * @name Hanging Nodes - * @{ - */ - - /** - * Compute the constraints resulting from - * the presence of hanging nodes. Hanging - * nodes are best explained using a small - * picture: - * - * @image html hanging_nodes.png - * - * In order to make a finite element - * function globally continuous, we have - * to make sure that the dark red nodes - * have values that are compatible with - * the adjacent yellow nodes, so that the - * function has no jump when coming from - * the small cells to the large one at - * the top right. We therefore have to - * add conditions that constrain those - * "hanging nodes". - * - * The object into - * which these are inserted is - * later used to condense the - * global system matrix and right - * hand side, and to extend the - * solution vectors from the true - * degrees of freedom also to the - * constraint nodes. This - * function is explained in - * detail in the @ref step_6 - * "step-6" tutorial program and - * is used in almost all - * following programs as well. - * - * This function does not clear - * the constraint matrix object - * before use, in order to allow - * adding constraints from - * different sources to the same - * object. You therefore need to - * make sure it contains only - * constraints you still want; - * otherwise call the - * ConstraintMatrix::clear() - * function. Likewise, this - * function does not close the - * object since you may want to - * enter other constraints later - * on yourself. - * - * In the hp-case, i.e. when the - * argument is of type - * hp::DoFHandler, we consider - * constraints due to different - * finite elements used on two - * sides of a face between cells - * as hanging nodes as well. In - * other words, for hp finite - * elements, this function - * computes all constraints due - * to differing mesh sizes (h) or - * polynomial degrees (p) between - * adjacent cells. - * - * The template argument (and by - * consequence the type of the - * first argument to this - * function) can be either a - * ::DoFHandler, hp::DoFHandler, - * or MGDoFHandler. - * - * @ingroup constraints - */ - template - static void - make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints); - //@} - - /** - * Take a vector of values which live on - * cells (e.g. an error per cell) and - * distribute it to the dofs in such a - * way that a finite element field - * results, which can then be further - * processed, e.g. for output. You should - * note that the resulting field will not - * be continuous at hanging nodes. This - * can, however, easily be arranged by - * calling the appropriate @p distribute - * function of a ConstraintMatrix - * object created for this - * DoFHandler object, after the - * vector has been fully assembled. - * - * It is assumed that the number - * of elements in @p cell_data - * equals the number of active - * cells and that the number of - * elements in @p dof_data equals - * dof_handler.n_dofs(). - * - * Note that the input vector may - * be a vector of any data type - * as long as it is convertible - * to @p double. The output - * vector, being a data vector on - * a DoF handler, always consists of - * elements of type @p double. - * - * In case the finite element - * used by this DoFHandler - * consists of more than one - * component, you need to specify - * which component in the output - * vector should be used to store - * the finite element field in; - * the default is zero (no other - * value is allowed if the finite - * element consists only of one - * component). All other - * components of the vector - * remain untouched, i.e. their - * contents are not changed. - * - * This function cannot be used - * if the finite element in use - * has shape functions that are - * non-zero in more than one - * vector component (in deal.II - * speak: they are - * non-primitive). - */ - template - static void - distribute_cell_to_dof_vector (const DH &dof_handler, - const Vector &cell_data, - Vector &dof_data, - const unsigned int component = 0); - - /** - * Extract the indices of the - * degrees of freedom belonging - * to certain vector components - * or blocks (if the last - * argument is true) of - * a vector-valued finite - * element. The bit vector @p - * select defines, which - * components or blocks of an - * FESystem are to be extracted - * from the DoFHandler @p - * dof. The entries in @p - * selected_dofs corresponding to - * degrees of freedom belonging - * to these components are then - * flagged @p true, while all - * others are set to @p false. - * - * The size of @p select must - * equal the number of components - * or blocks in the FiniteElement - * used by @p dof, depending on - * the argument - * blocks. The size of - * @p selected_dofs must equal - * DoFHandler::n_dofs(). Previous - * contents of this array are - * overwritten. - * - * If the finite element under - * consideration is not - * primitive, that is some or all - * of its shape functions are - * non-zero in more than one - * vector component (which holds, - * for example, for FE_Nedelec or - * FE_RaviartThomas elements), then - * shape functions cannot be - * associated with a single - * vector component. In this - * case, if one shape - * vector component of this - * element is flagged in - * @p component_select, then - * this is equivalent to - * selecting all vector - * components corresponding to - * this non-primitive base - * element. - */ - template - static void - extract_dofs (const DoFHandler &dof_handler, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); - - /** - * The same function as above, - * but for a hp::DoFHandler. - */ - template - static void - extract_dofs (const hp::DoFHandler &dof_handler, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); - - /** - * Do the same thing as - * extract_dofs() for one level - * of a multi-grid DoF numbering. - */ - template - static void - extract_level_dofs (const unsigned int level, - const MGDoFHandler &dof, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); - - /** - * Extract all degrees of freedom - * which are at the boundary and - * belong to specified components - * of the solution. The function - * returns its results in the - * last non-default-valued - * parameter which contains - * @p true if a degree of - * freedom is at the boundary and - * belongs to one of the selected - * components, and @p false - * otherwise. - * - * By specifying the - * @p boundary_indicator - * variable, you can select which - * boundary indicators the faces - * have to have on which the - * degrees of freedom are located - * that shall be extracted. If it - * is an empty list, then all - * boundary indicators are - * accepted. - * - * The size of @p component_select - * shall equal the number of - * components in the finite - * element used by @p dof. The - * size of @p selected_dofs shall - * equal - * dof_handler.n_dofs(). Previous - * contents of this array or - * overwritten. - * - * Using the usual convention, if - * a shape function is non-zero - * in more than one component - * (i.e. it is non-primitive), - * then the element in the - * component mask is used that - * corresponds to the first - * non-zero components. Elements - * in the mask corresponding to - * later components are ignored. - */ - template - static void - extract_boundary_dofs (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); - - /** - * This function is similar to - * the extract_boundary_dofs() - * function but it extracts those - * degrees of freedom whose shape - * functions are nonzero on at - * least part of the selected - * boundary. For continuous - * elements, this is exactly the - * set of shape functions whose - * degrees of freedom are defined - * on boundary faces. On the - * other hand, if the finite - * element in used is a - * discontinuous element, all - * degrees of freedom are defined - * in the inside of cells and - * consequently none would be - * boundary degrees of - * freedom. Several of those - * would have shape functions - * that are nonzero on the - * boundary, however. This - * function therefore extracts - * all those for which the - * FiniteElement::has_support_on_face - * function says that it is - * nonzero on any face on one of - * the selected boundary parts. - */ - template - static void - extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); - - /** - * @name Hanging Nodes - * @{ - */ - - /** - * Select all dofs that will be - * constrained by interface - * constraints, i.e. all hanging - * nodes. - * - * The size of @p selected_dofs - * shall equal - * dof_handler.n_dofs(). Previous - * contents of this array or - * overwritten. - */ - template - static void - extract_hanging_node_dofs (const DoFHandler &dof_handler, - std::vector &selected_dofs); - //@} - - /** - * Flag all those degrees of - * freedom which are on cells - * with the given subdomain - * id. Note that DoFs on faces - * can belong to cells with - * differing subdomain ids, so - * the sets of flagged degrees of - * freedom are not mutually - * exclusive for different - * subdomain ids. - * - * If you want to get a unique - * association of degree of freedom with - * subdomains, use the - * @p get_subdomain_association - * function. - */ - template - static void - extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id_t subdomain_id, - std::vector &selected_dofs); - - - /** - * Extract the set of global DoF - * indices that are owned by the - * current processor. For regular - * DoFHandler objects, this set - * is the complete set with all - * DoF indices. In either case, - * it equals what - * DoFHandler::locally_owned_dofs() - * returns. - */ - template - static void - extract_locally_owned_dofs (const DH & dof_handler, - IndexSet & dof_set); - - - /** - * Extract the set of global DoF - * indices that are active on the - * current DoFHandler. For - * regular DoFHandlers, these are - * all DoF indices, but for - * DoFHandler objects built on - * parallel::distributed::Triangulation - * this set is a superset of - * DoFHandler::locally_owned_dofs() - * and contains all DoF indices - * that live on all locally owned - * cells (including on the - * interface to ghost - * cells). However, it does not - * contain the DoF indices that - * are exclusively defined on - * ghost or artificial cells (see - * @ref GlossArtificialCell "the - * glossary"). - * - * The degrees of freedom identified by - * this function equal those obtained - * from the - * dof_indices_with_subdomain_association() - * function when called with the locally - * owned subdomain id. - */ - template - static void - extract_locally_active_dofs (const DH & dof_handler, + /** + * The flags used in tables by certain + * make_*_pattern functions to + * describe whether two components of the + * solution couple in the bilinear forms + * corresponding to cell or face + * terms. An example of using these flags + * is shown in the introduction of + * step-46. + * + * In the descriptions of the individual + * elements below, remember that these + * flags are used as elements of tables + * of size FiniteElement::n_components + * times FiniteElement::n_components + * where each element indicates whether + * two components do or do not couple. + */ + enum Coupling + { + /** + * Two components do not + * couple. + */ + none, + /** + * Two components do couple. + */ + always, + /** + * 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. + */ + nonzero + }; + + /** + * @name Auxiliary functions + * @{ + */ + /** + * Maximal number of degrees of + * freedom on a cell. + */ + template + unsigned int + max_dofs_per_cell (const DoFHandler &dh); + + template + unsigned int + max_dofs_per_cell (const hp::DoFHandler &dh); + + + /** + * Maximal number of degrees of + * freedom on a face. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + max_dofs_per_face (const DoFHandler &dh); + + /** + * Maximal number of degrees of + * freedom on a face. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + max_dofs_per_face (const hp::DoFHandler &dh); + + /** + * Maximal number of degrees of + * freedom on a vertex. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + max_dofs_per_vertex (const DoFHandler &dh); + + /** + * Maximal number of degrees of + * freedom on a vertex. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + max_dofs_per_vertex (const hp::DoFHandler &dh); + + /** + * Number of vector components in the + * finite element object used by this + * DoFHandler. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + n_components (const DoFHandler &dh); + + /** + * Number of vector components in the + * finite element object used by this + * DoFHandler. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + unsigned int + n_components (const hp::DoFHandler &dh); + + /** + * Find out whether the FiniteElement + * used by this DoFHandler is primitive + * or not. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + bool + fe_is_primitive (const DoFHandler &dh); + + /** + * Find out whether the FiniteElement + * used by this DoFHandler is primitive + * or not. + * + * This function exists for both non-hp + * and hp DoFHandlers, to allow for a + * uniform interface to query this + * property. + */ + template + bool + fe_is_primitive (const hp::DoFHandler &dh); + + /** + * @} + */ + + /** + * @name Sparsity Pattern Generation + * @{ + */ + + /** + * Locate non-zero entries of the + * system matrix. + * + * This function computes the + * possible positions of non-zero + * entries in the global system + * matrix. We assume that a + * certain 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, + * CompressedSparsityPattern, + * BlockSparsityPattern, + * BlockCompressedSparsityPattern, + * BlockCompressedSetSparsityPattern, + * BlockCompressedSimpleSparsityPattern, + * 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. + * + * @ingroup constraints + */ + template + void + make_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, + const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); + + /** + * Locate non-zero entries for + * vector valued finite elements. + * This function does mostly the + * same as the previous + * @p 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, + * @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: + * @f[ + * \begin{array}{ccc} + * 1 & 0 & 1 \\ + * 0 & 1 & 1 \\ + * 1 & 1 & 0 + * \end{array} + * @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. + * + * This function is designed to + * accept a coupling pattern, like the one + * shown above, through the + * @p couplings parameter, which + * contains values of type #Coupling. It + * builds the matrix structure + * just like the previous + * function, but does not create + * matrix elements if not + * specified by the 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, + * CompressedSparsityPattern, + * BlockSparsityPattern, + * BlockCompressedSparsityPattern, + * BlockCompressedSetSparsityPattern, + * 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 + * correspoding to the first + * non-zero component is taken + * and additional ones for this + * component are ignored. + * + * 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 DH &dof, + const Table<2, Coupling> &coupling, + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, + const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); + + /** + * @deprecated This is the old + * form of the previous + * function. It generates a table + * of DoFTools::Coupling values + * (where a true + * value in the mask is + * translated into a + * Coupling::always value in the + * table) and calls the function + * above. + */ + template + void + make_sparsity_pattern (const DH &dof, + const std::vector > &mask, + SparsityPattern &sparsity_pattern); + + /** + * Construct a sparsity pattern that + * allows coupling degrees of freedom on + * two different but related meshes. + * + * The idea is that if the two given + * DoFHandler objects correspond to two + * different meshes (and potentially to + * different finite elements used on + * these cells), but that if the two + * triangulations they are based on are + * derived from the same coarse mesh + * through hierarchical refinement, then + * one may set up a problem where one + * would like to test shape functions + * from one mesh against the shape + * functions from another mesh. In + * particular, this means that shape + * functions from a cell on the first + * mesh are tested against those on the + * second cell that are located on the + * corresponding cell; this + * correspondence is something that the + * IntergridMap class can determine. + * + * This function then constructs a + * sparsity pattern for which the degrees + * of freedom that represent the rows + * come from the first given DoFHandler, + * whereas the ones that correspond to + * columns come from the second + * DoFHandler. + */ + template + void + make_sparsity_pattern (const DH &dof_row, + const DH &dof_col, + SparsityPattern &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, + * CompressedSparsityPattern, + * BlockSparsityPattern, + * BlockCompressedSparsityPattern, + * BlockCompressedSetSparsityPattern, 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 DH &dof, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &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. + * + * 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 DH &dof, + const typename FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &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. + */ + template + void + make_flux_sparsity_pattern (const DH &dof_handler, + SparsityPattern &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. + * Not implemented for + * hp::DoFHandler. + * + * @ingroup constraints + */ + template + void + make_flux_sparsity_pattern (const DH &dof_handler, + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs = true, + const types::subdomain_id_t 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. + * + * There is one matrix for + * couplings in a cell and one + * for the couplings occuring in + * fluxes. + * + * Not implemented for + * hp::DoFHandler. + */ + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const Table<2,Coupling> &int_mask, + const Table<2,Coupling> &flux_mask); + + //@} + /** + * @name Hanging Nodes + * @{ + */ + + /** + * Compute the constraints resulting from + * the presence of hanging nodes. Hanging + * nodes are best explained using a small + * picture: + * + * @image html hanging_nodes.png + * + * In order to make a finite element + * function globally continuous, we have + * to make sure that the dark red nodes + * have values that are compatible with + * the adjacent yellow nodes, so that the + * function has no jump when coming from + * the small cells to the large one at + * the top right. We therefore have to + * add conditions that constrain those + * "hanging nodes". + * + * The object into + * which these are inserted is + * later used to condense the + * global system matrix and right + * hand side, and to extend the + * solution vectors from the true + * degrees of freedom also to the + * constraint nodes. This + * function is explained in + * detail in the @ref step_6 + * "step-6" tutorial program and + * is used in almost all + * following programs as well. + * + * This function does not clear + * the constraint matrix object + * before use, in order to allow + * adding constraints from + * different sources to the same + * object. You therefore need to + * make sure it contains only + * constraints you still want; + * otherwise call the + * ConstraintMatrix::clear() + * function. Likewise, this + * function does not close the + * object since you may want to + * enter other constraints later + * on yourself. + * + * In the hp-case, i.e. when the + * argument is of type + * hp::DoFHandler, we consider + * constraints due to different + * finite elements used on two + * sides of a face between cells + * as hanging nodes as well. In + * other words, for hp finite + * elements, this function + * computes all constraints due + * to differing mesh sizes (h) or + * polynomial degrees (p) between + * adjacent cells. + * + * The template argument (and by + * consequence the type of the + * first argument to this + * function) can be either a + * ::DoFHandler, hp::DoFHandler, + * or MGDoFHandler. + * + * @ingroup constraints + */ + template + void + make_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints); + //@} + + /** + * Take a vector of values which live on + * cells (e.g. an error per cell) and + * distribute it to the dofs in such a + * way that a finite element field + * results, which can then be further + * processed, e.g. for output. You should + * note that the resulting field will not + * be continuous at hanging nodes. This + * can, however, easily be arranged by + * calling the appropriate @p distribute + * function of a ConstraintMatrix + * object created for this + * DoFHandler object, after the + * vector has been fully assembled. + * + * It is assumed that the number + * of elements in @p cell_data + * equals the number of active + * cells and that the number of + * elements in @p dof_data equals + * dof_handler.n_dofs(). + * + * Note that the input vector may + * be a vector of any data type + * as long as it is convertible + * to @p double. The output + * vector, being a data vector on + * a DoF handler, always consists of + * elements of type @p double. + * + * In case the finite element + * used by this DoFHandler + * consists of more than one + * component, you need to specify + * which component in the output + * vector should be used to store + * the finite element field in; + * the default is zero (no other + * value is allowed if the finite + * element consists only of one + * component). All other + * components of the vector + * remain untouched, i.e. their + * contents are not changed. + * + * This function cannot be used + * if the finite element in use + * has shape functions that are + * non-zero in more than one + * vector component (in deal.II + * speak: they are + * non-primitive). + */ + template + void + distribute_cell_to_dof_vector (const DH &dof_handler, + const Vector &cell_data, + Vector &dof_data, + const unsigned int component = 0); + + /** + * Extract the indices of the + * degrees of freedom belonging + * to certain vector components + * or blocks (if the last + * argument is true) of + * a vector-valued finite + * element. The bit vector @p + * select defines, which + * components or blocks of an + * FESystem are to be extracted + * from the DoFHandler @p + * dof. The entries in @p + * selected_dofs corresponding to + * degrees of freedom belonging + * to these components are then + * flagged @p true, while all + * others are set to @p false. + * + * The size of @p select must + * equal the number of components + * or blocks in the FiniteElement + * used by @p dof, depending on + * the argument + * blocks. The size of + * @p selected_dofs must equal + * DoFHandler::n_dofs(). Previous + * contents of this array are + * overwritten. + * + * If the finite element under + * consideration is not + * primitive, that is some or all + * of its shape functions are + * non-zero in more than one + * vector component (which holds, + * for example, for FE_Nedelec or + * FE_RaviartThomas elements), then + * shape functions cannot be + * associated with a single + * vector component. In this + * case, if one shape + * vector component of this + * element is flagged in + * @p component_select, then + * this is equivalent to + * selecting all vector + * components corresponding to + * this non-primitive base + * element. + */ + template + void + extract_dofs (const DoFHandler &dof_handler, + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); + + /** + * The same function as above, + * but for a hp::DoFHandler. + */ + template + void + extract_dofs (const hp::DoFHandler &dof_handler, + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); + + /** + * Do the same thing as + * extract_dofs() for one level + * of a multi-grid DoF numbering. + */ + template + void + extract_level_dofs (const unsigned int level, + const MGDoFHandler &dof, + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); + + /** + * Extract all degrees of freedom + * which are at the boundary and + * belong to specified components + * of the solution. The function + * returns its results in the + * last non-default-valued + * parameter which contains + * @p true if a degree of + * freedom is at the boundary and + * belongs to one of the selected + * components, and @p false + * otherwise. + * + * By specifying the + * @p boundary_indicator + * variable, you can select which + * boundary indicators the faces + * have to have on which the + * degrees of freedom are located + * that shall be extracted. If it + * is an empty list, then all + * boundary indicators are + * accepted. + * + * The size of @p component_select + * shall equal the number of + * components in the finite + * element used by @p dof. The + * size of @p selected_dofs shall + * equal + * dof_handler.n_dofs(). Previous + * contents of this array or + * overwritten. + * + * Using the usual convention, if + * a shape function is non-zero + * in more than one component + * (i.e. it is non-primitive), + * then the element in the + * component mask is used that + * corresponds to the first + * non-zero components. Elements + * in the mask corresponding to + * later components are ignored. + */ + template + void + extract_boundary_dofs (const DH &dof_handler, + const std::vector &component_select, + std::vector &selected_dofs, + const std::set &boundary_indicators = std::set()); + + /** + * This function is similar to + * the extract_boundary_dofs() + * function but it extracts those + * degrees of freedom whose shape + * functions are nonzero on at + * least part of the selected + * boundary. For continuous + * elements, this is exactly the + * set of shape functions whose + * degrees of freedom are defined + * on boundary faces. On the + * other hand, if the finite + * element in used is a + * discontinuous element, all + * degrees of freedom are defined + * in the inside of cells and + * consequently none would be + * boundary degrees of + * freedom. Several of those + * would have shape functions + * that are nonzero on the + * boundary, however. This + * function therefore extracts + * all those for which the + * FiniteElement::has_support_on_face + * function says that it is + * nonzero on any face on one of + * the selected boundary parts. + */ + template + void + extract_dofs_with_support_on_boundary (const DH &dof_handler, + const std::vector &component_select, + std::vector &selected_dofs, + const std::set &boundary_indicators = std::set()); + + /** + * @name Hanging Nodes + * @{ + */ + + /** + * Select all dofs that will be + * constrained by interface + * constraints, i.e. all hanging + * nodes. + * + * The size of @p selected_dofs + * shall equal + * dof_handler.n_dofs(). Previous + * contents of this array or + * overwritten. + */ + template + void + extract_hanging_node_dofs (const DoFHandler &dof_handler, + std::vector &selected_dofs); + //@} + + /** + * Flag all those degrees of + * freedom which are on cells + * with the given subdomain + * id. Note that DoFs on faces + * can belong to cells with + * differing subdomain ids, so + * the sets of flagged degrees of + * freedom are not mutually + * exclusive for different + * subdomain ids. + * + * If you want to get a unique + * association of degree of freedom with + * subdomains, use the + * @p get_subdomain_association + * function. + */ + template + void + extract_subdomain_dofs (const DH &dof_handler, + const types::subdomain_id_t subdomain_id, + std::vector &selected_dofs); + + + /** + * Extract the set of global DoF + * indices that are owned by the + * current processor. For regular + * DoFHandler objects, this set + * is the complete set with all + * DoF indices. In either case, + * it equals what + * DoFHandler::locally_owned_dofs() + * returns. + */ + template + void + extract_locally_owned_dofs (const DH & dof_handler, + IndexSet & dof_set); + + + /** + * Extract the set of global DoF + * indices that are active on the + * current DoFHandler. For + * regular DoFHandlers, these are + * all DoF indices, but for + * DoFHandler objects built on + * parallel::distributed::Triangulation + * this set is a superset of + * DoFHandler::locally_owned_dofs() + * and contains all DoF indices + * that live on all locally owned + * cells (including on the + * interface to ghost + * cells). However, it does not + * contain the DoF indices that + * are exclusively defined on + * ghost or artificial cells (see + * @ref GlossArtificialCell "the + * glossary"). + * + * The degrees of freedom identified by + * this function equal those obtained + * from the + * dof_indices_with_subdomain_association() + * function when called with the locally + * owned subdomain id. + */ + template + void + extract_locally_active_dofs (const DH & dof_handler, + IndexSet & dof_set); + + /** + * Extract the set of global DoF + * indices that are active on the + * current DoFHandler. For + * regular DoFHandlers, these are + * all DoF indices, but for + * DoFHandler objects built on + * parallel::distributed::Triangulation + * this set is the union of + * DoFHandler::locally_owned_dofs() + * and the DoF indices on all + * ghost cells. In essence, it is + * the DoF indices on all cells + * that are not artificial (see + * @ref GlossArtificialCell "the glossary"). + */ + template + void + extract_locally_relevant_dofs (const DH & dof_handler, IndexSet & dof_set); - /** - * Extract the set of global DoF - * indices that are active on the - * current DoFHandler. For - * regular DoFHandlers, these are - * all DoF indices, but for - * DoFHandler objects built on - * parallel::distributed::Triangulation - * this set is the union of - * DoFHandler::locally_owned_dofs() - * and the DoF indices on all - * ghost cells. In essence, it is - * the DoF indices on all cells - * that are not artificial (see - * @ref GlossArtificialCell "the glossary"). - */ - template - static void - extract_locally_relevant_dofs (const DH & dof_handler, - IndexSet & dof_set); - - /** - * Extract a vector that represents the - * constant modes of the DoFHandler for - * the components chosen by - * component_select. The - * constant modes on a discretization are - * the null space of a Laplace operator - * on the selected components with - * Neumann boundary conditions - * applied. The null space is a necessary - * ingredient for obtaining a good AMG - * preconditioner when using the class - * TrilinosWrappers::PreconditionAMG. - * Since the ML AMG package only works on - * algebraic properties of the respective - * matrix, it has no chance to detect - * whether the matrix comes from a scalar - * or a vector valued problem. However, a - * near null space supplies exactly the - * needed information about these - * components. The null space will - * consist of as many vectors as there - * are true arguments in - * component_select, each of - * which will be one in one vector component and - * zero in all others. We store this - * object in a vector of vectors, where - * the outer vector is of the size of the - * number of selected components, and - * each inner vector has as many - * components as there are (locally owned) degrees of - * freedom in the selected - * components. Note that any matrix - * associated with this null space must - * have been constructed using the same - * component_select argument, - * since the numbering of DoFs is done - * relative to the selected dofs, not to - * all dofs. - * - * The main reason for this - * program is the use of the - * null space with the - * AMG preconditioner. - */ - template - static void - extract_constant_modes (const DH &dof_handler, - const std::vector &component_select, - std::vector > &constant_modes); - - /** - * For each active cell of a DoFHandler - * or hp::DoFHandler, extract the active - * finite element index and fill the - * vector given as second argument. This - * vector is assumed to have as many - * entries as there are active cells. - * - * For non-hp DoFHandler objects given as - * first argument, the returned vector - * will consist of only zeros, indicating - * that all cells use the same finite - * element. For a hp::DoFHandler, the - * values may be different, though. - */ - template - static void - get_active_fe_indices (const DH &dof_handler, - std::vector &active_fe_indices); - - /** - * 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. - * - * Note that degrees of freedom - * associated with faces, edges, and - * vertices may be associated with - * multiple subdomains if they are - * sitting on partition boundaries. In - * these cases, we put them into one of - * the associated partitions in an - * undefined way. This may sometimes lead - * to different numbers of degrees of - * freedom in partitions, even if the - * number of cells is perfectly - * equidistributed. While this is - * regrettable, it is not a problem in - * practice since the number of degrees - * of freedom on partition boundaries is - * asymptotically vanishing as we refine - * the mesh as long as the number of - * partitions is kept constant. - * - * This function returns the association - * of each DoF with one subdomain. If you - * are looking for the association of - * each @em cell with a subdomain, either - * query the - * cell->subdomain_id() - * function, or use the - * GridTools::get_subdomain_association - * function. - * - * Note that this function is of - * questionable use for DoFHandler objects built on - * parallel::distributed::Triangulation - * since in that case ownership of - * individual degrees of freedom by MPI - * processes is controlled by the DoF - * handler object, not based on some - * geometric algorithm in conjunction - * with subdomain id. In particular, the - * degrees of freedom identified by the - * functions in this namespace as - * associated with a subdomain are not - * the same the - * DoFHandler class - * identifies as those it owns. - */ - template - static void - get_subdomain_association (const DH &dof_handler, - std::vector &subdomain); - - /** - * Count how many degrees of freedom are - * uniquely associated with the given - * @p subdomain index. - * - * Note that there may be rare cases - * where cells with the given @p - * subdomain index exist, but none of its - * degrees of freedom are actually - * associated with it. In that case, the - * returned value will be zero. - * - * This function will generate an - * exception if there are no cells with - * the given @p subdomain index. - * - * This function returns the number of - * DoFs associated with one subdomain. If - * you are looking for the association of - * @em cells with this subdomain, use the - * GridTools::count_cells_with_subdomain_association - * function. - * - * Note that this function is of - * questionable use for DoFHandler objects built on - * parallel::distributed::Triangulation - * since in that case ownership of - * individual degrees of freedom by MPI - * processes is controlled by the DoF - * handler object, not based on some - * geometric algorithm in conjunction - * with subdomain id. In particular, the - * degrees of freedom identified by the - * functions in this namespace as - * associated with a subdomain are not - * the same the - * DoFHandler class - * identifies as those it owns. - */ - template - static unsigned int - count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain); - - /** - * Count how many degrees of freedom are - * uniquely associated with the given - * @p subdomain index. - * - * This function does what the previous - * one does except that it splits the - * result among the vector components of - * the finite element in use by the - * DoFHandler object. The last argument - * (which must have a length equal to the - * number of vector components) will - * therefore store how many degrees of - * freedom of each vector component are - * associated with the given subdomain. - * - * Note that this function is of - * questionable use for DoFHandler objects built on - * parallel::distributed::Triangulation - * since in that case ownership of - * individual degrees of freedom by MPI - * processes is controlled by the DoF - * handler object, not based on some - * geometric algorithm in conjunction - * with subdomain id. In particular, the - * degrees of freedom identified by the - * functions in this namespace as - * associated with a subdomain are not - * the same the - * DoFHandler class - * identifies as those it owns. - */ - template - static void - count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain, - std::vector &n_dofs_on_subdomain); - - /** - * Return a set of indices that denotes - * the degrees of freedom that live on - * the given subdomain, i.e. that are on - * cells owned by the current - * processor. Note that this includes the - * ones that this subdomain "owns" - * (i.e. the ones for which - * get_subdomain_association() returns a - * value equal to the subdomain given - * here and that are selected by the - * extract_locally_owned() function) but - * also all of those that sit on the - * boundary between the given subdomain - * and other subdomain. In essence, - * degrees of freedom that sit on - * boundaries between subdomain will be - * in the index sets returned by this - * function for more than one subdomain. - * - * Note that this function is of - * questionable use for DoFHandler objects built on - * parallel::distributed::Triangulation - * since in that case ownership of - * individual degrees of freedom by MPI - * processes is controlled by the DoF - * handler object, not based on some - * geometric algorithm in conjunction - * with subdomain id. In particular, the - * degrees of freedom identified by the - * functions in this namespace as - * associated with a subdomain are not - * the same the - * DoFHandler class - * identifies as those it owns. - */ - template - static - IndexSet - dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain); - - /** - * Count how many degrees of - * freedom out of the total - * number belong to each - * component. If the number of - * components the finite element - * has is one (i.e. you only have - * one scalar variable), then the - * number in this component - * obviously equals the total - * number of degrees of - * freedom. Otherwise, the sum of - * the DoFs in all the components - * needs to equal the total - * number. - * - * However, the last statement - * does not hold true if the - * finite element is not - * primitive, i.e. some or all of - * its shape functions are - * non-zero in more than one - * vector component. This - * applies, for example, to the - * Nedelec or Raviart-Thomas - * elements. In this case, a - * degree of freedom is counted - * in each component in which it - * is non-zero, so that the sum - * mentioned above is greater - * than the total number of - * degrees of freedom. - * - * This behavior can be switched - * off by the optional parameter - * vector_valued_once. If - * this is true, the - * number of components of a - * nonprimitive vector valued - * element is collected only in - * the first component. All other - * components will have a count - * of zero. - * - * The additional optional - * argument @p target_component - * allows for a re-sorting and - * grouping of components. To - * this end, it contains for each - * component the component number - * it shall be counted as. Having - * the same number entered - * several times sums up several - * components as the same. One of - * the applications of this - * argument is when you want to - * form block matrices and - * vectors, but want to pack - * several components into the - * same block (for example, when - * you have @p dim velocities - * and one pressure, to put all - * velocities into one block, and - * the pressure into another). - * - * The result is returned in @p - * dofs_per_component. Note that - * the size of @p - * dofs_per_component needs to be - * enough to hold all the indices - * specified in @p - * target_component. If this is - * not the case, an assertion is - * thrown. The indices not - * targetted by target_components - * are left untouched. - */ - template - static void - count_dofs_per_component (const DoFHandler& dof_handler, - std::vector& dofs_per_component, - const bool vector_valued_once = false, - std::vector target_component - = std::vector()); - - /** - * Count the degrees of freedom - * in each block. This function - * is similar to - * count_dofs_per_component(), - * with the difference that the - * counting is done by - * blocks. See @ref GlossBlock - * "blocks" in the glossary for - * details. Again the vectors are - * assumed to have the correct - * size before calling this - * function. If this is not the - * case, an assertion is thrown. - * - * This function is used in the - * step-22, - * step-31, and - * step-32 tutorial - * programs. - * - * @pre The dofs_per_block - * variable has as many - * components as the finite - * element used by the - * dof_handler argument has - * blocks, or alternatively as - * many blocks as are enumerated - * in the target_blocks argument - * if given. - */ - template - static void - count_dofs_per_block (const DoFHandler& dof_handler, - std::vector& dofs_per_block, - std::vector target_blocks - = std::vector()); - - /** - * @deprecated See the previous - * function with the same name - * for a description. This - * function exists for - * compatibility with older - * versions only. - */ - template - static void - count_dofs_per_component (const DoFHandler& dof_handler, - std::vector& dofs_per_component, - std::vector target_component); - - /** - * This function can be used when - * different variables shall be - * discretized on different - * grids, where one grid is - * coarser than the other. This - * idea might seem nonsensical at - * first, but has reasonable - * applications in inverse - * (parameter estimation) - * problems, where there might - * not be enough information to - * recover the parameter on the - * same grid as the state - * variable; furthermore, the - * smoothness properties of state - * variable and parameter might - * not be too much related, so - * using different grids might be - * an alternative to using - * stronger regularization of the - * problem. - * - * The basic idea of this - * function is explained in the - * following. Let us, for - * convenience, denote by - * ``parameter grid'' the coarser - * of the two grids, and by - * ``state grid'' the finer of - * the two. We furthermore assume - * that the finer grid can be - * obtained by refinement of the - * coarser one, i.e. the fine - * grid is at least as much - * refined as the coarse grid at - * each point of the - * domain. Then, each shape - * function on the coarse grid - * can be represented as a linear - * combination of shape functions - * on the fine grid (assuming - * identical ansatz - * spaces). Thus, if we - * discretize as usual, using - * shape functions on the fine - * grid, we can consider the - * restriction that the parameter - * variable shall in fact be - * discretized by shape functions - * on the coarse grid as a - * constraint. These constraints - * are linear and happen to have - * the form managed by the - * ``ConstraintMatrix'' class. - * - * The construction of these - * constraints is done as - * follows: for each of the - * degrees of freedom (i.e. shape - * functions) on the coarse grid, - * we compute its representation - * on the fine grid, i.e. how the - * linear combination of shape - * functions on the fine grid - * looks like that resembles the - * shape function on the coarse - * grid. From this information, - * we can then compute the - * constraints which have to hold - * if a solution of a linear - * equation on the fine grid - * shall be representable on the - * coarse grid. The exact - * algorithm how these - * constraints can be computed is - * rather complicated and is best - * understood by reading the - * source code, which contains - * many comments. - * - * Before explaining the use of - * this function, we would like - * to state that the total number - * of degrees of freedom used for - * the discretization is not - * reduced by the use of this - * function, i.e. even though we - * discretize one variable on a - * coarser grid, the total number - * of degrees of freedom is that - * of the fine grid. This seems - * to be counter-productive, - * since it does not give us a - * benefit from using a coarser - * grid. The reason why it may be - * useful to choose this approach - * nonetheless is three-fold: - * first, as stated above, there - * might not be enough - * information to recover a - * parameter on a fine grid, - * i.e. we chose to discretize it - * on the coarse grid not to save - * DoFs, but for other - * reasons. Second, the - * ``ConstraintMatrix'' includes - * the constraints into the - * linear system of equations, by - * which constrained nodes become - * dummy nodes; we may therefore - * exclude them from the linear - * algebra, for example by - * sorting them to the back of - * the DoF numbers and simply - * calling the solver for the - * upper left block of the matrix - * which works on the - * non-constrained nodes only, - * thus actually realizing the - * savings in numerical effort - * from the reduced number of - * actual degrees of freedom. The - * third reason is that for some - * or other reason we have chosen - * to use two different grids, it - * may be actually quite - * difficult to write a function - * that assembles the system - * matrix for finite element - * spaces on different grids; - * using the approach of - * constraints as with this - * function allows to use - * standard techniques when - * discretizing on only one grid - * (the finer one) without having - * to take care of the fact that - * one or several of the variable - * actually belong to different - * grids. - * - * The use of this function is as - * follows: it accepts as - * parameters two DoF Handlers, - * the first of which refers to - * the coarse grid and the second - * of which is the fine grid. On - * both, a finite element is - * represented by the DoF handler - * objects, which will usually - * have several components, which - * may belong to different finite - * elements. The second and - * fourth parameter of this - * function therefore state which - * variable on the coarse grid - * shall be used to restrict the - * stated component on the fine - * grid. Of course, the finite - * elements used for the - * respective components on the - * two grids need to be the - * same. An example may clarify - * this: consider the parameter - * estimation mentioned briefly - * above; there, on the fine grid - * the whole discretization is - * done, thus the variables are - * ``u'', ``q'', and the Lagrange - * multiplier ``lambda'', which - * are discretized using - * continuous linear, piecewise - * constant discontinuous, and - * continuous linear elements, - * respectively. Only the - * parameter ``q'' shall be - * represented on the coarse - * grid, thus the DoFHandler - * object on the coarse grid - * represents only one variable, - * discretized using piecewise - * constant discontinuous - * elements. Then, the parameter - * denoting the component on the - * coarse grid would be zero (the - * only possible choice, since - * the variable on the coarse - * grid is scalar), and one on - * the fine grid (corresponding - * to the variable ``q''; zero - * would be ``u'', two would be - * ``lambda''). Furthermore, an - * object of type IntergridMap - * is needed; this could in - * principle be generated by the - * function itself from the two - * DoFHandler objects, but since - * it is probably available - * anyway in programs that use - * this function, we shall use it - * instead of re-generating - * it. Finally, the computed - * constraints are entered into a - * variable of type - * ConstraintMatrix; the - * constraints are added, - * i.e. previous contents which - * may have, for example, be - * obtained from hanging nodes, - * are not deleted, so that you - * only need one object of this - * type. - */ - template - static void - compute_intergrid_constraints (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - ConstraintMatrix &constraints); - - - /** - * This function generates a - * matrix such that when a vector - * of data with as many elements - * as there are degrees of - * freedom of this component on - * the coarse grid is multiplied - * to this matrix, we obtain a - * vector with as many elements - * are there are global degrees - * of freedom on the fine - * grid. All the elements of the - * other components of the finite - * element fields on the fine - * grid are not touched. - * - * The output of this function is - * a compressed format that can - * be given to the @p reinit - * functions of the - * SparsityPattern ad - * SparseMatrix classes. - */ - template - static void - compute_intergrid_transfer_representation (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - std::vector > &transfer_representation); - - /** - * Create a mapping from degree - * of freedom indices to the - * index of that degree of - * freedom on the boundary. After - * this operation, mapping[dof] - * gives the index of the - * degree of freedom with global - * number @p dof in the list of - * degrees of freedom on the - * boundary. If the degree of - * freedom requested is not on - * the boundary, the value of - * mapping[dof] is - * @p invalid_dof_index. This - * function is mainly used when - * setting up matrices and - * vectors on the boundary from - * the trial functions, which - * have global numbers, while the - * matrices and vectors use - * numbers of the trial functions - * local to the boundary. - * - * Prior content of @p mapping - * is deleted. - */ - template - static void - map_dof_to_boundary_indices (const DH &dof_handler, - std::vector &mapping); - - /** - * Same as the previous function, - * except that only those parts - * of the boundary are considered - * for which the boundary - * indicator is listed in the - * second argument. - * - * See the general doc of this - * class for more information. - */ - template - static void - map_dof_to_boundary_indices (const DH &dof_handler, - const std::set &boundary_indicators, - std::vector &mapping); - - /** - * Return a list of support - * points for all the degrees of - * freedom handled by this DoF - * handler object. This function, - * of course, only works if the - * finite element object used by - * the DoF handler object - * actually provides support - * points, i.e. no edge elements - * or the like. Otherwise, an - * exception is thrown. - * - * The given array must have a - * length of as many elements as - * there are degrees of freedom. - */ - template - static void - map_dofs_to_support_points (const Mapping &mapping, - const DoFHandler &dof_handler, - std::vector > &support_points); - - /** - * This is the opposite function - * to the one above. It generates - * a map where the keys are the - * support points of the degrees - * of freedom, while the values - * are the DoF indices. - * - * Since there is no natural - * order in the space of points - * (except for the 1d case), you - * have to provide a map with an - * explicitly specified - * comparator object. This - * function is therefore - * templatized on the comparator - * object. Previous content of - * the map object is deleted in - * this function. - * - * Just as with the function - * above, it is assumed that the - * finite element in use here - * actually supports the notion - * of support points of all its - * components. - */ - template - static void - map_support_points_to_dofs (const Mapping &mapping, - const DH &dof_handler, - std::map, unsigned int, Comp> &point_to_index_map); - - /** - * Map a coupling table from the - * user friendly organization by - * components to the organization - * by blocks. Specializations of - * this function for DoFHandler - * and hp::DoFHandler are - * required due to the different - * results of their finite - * element access. - * - * The return vector will be - * initialized to the correct - * length inside this function. - */ - template - static void - convert_couplings_to_blocks (const hp::DoFHandler& dof_handler, - const Table<2, Coupling>& table_by_component, - std::vector >& tables_by_block); - - /** - * Make a constraint matrix for the - * constraints that result from zero - * boundary values. - * - * This function constrains all - * degrees of freedom on the - * boundary. Optionally, you can - * add a component mask, which - * restricts this functionality - * to a subset of an FESystem. - * - * For non-@ref GlossPrimitive "primitive" - * shape functions, any degree of freedom - * is affected that belongs to a - * shape function where at least - * one of its nonzero components - * is affected. - * - * This function is used - * in step-36, for - * example. - * - * @ingroup constraints - */ - template class DH> - static void - make_zero_boundary_constraints (const DH &dof, - ConstraintMatrix &zero_boundary_constraints, - const std::vector &component_mask_=std::vector()); - - /** - * Map a coupling table from the - * user friendly organization by - * components to the organization - * by blocks. Specializations of - * this function for DoFHandler - * and hp::DoFHandler are - * required due to the different - * results of their finite - * element access. - * - * The return vector will be - * initialized to the correct - * length inside this function. - */ - template - static void - convert_couplings_to_blocks (const DoFHandler& dof_handler, - const Table<2, Coupling>& table_by_component, - std::vector >& tables_by_block); - - /** - * Given a finite element and a table how - * the vector components of it couple - * with each other, compute and return a - * table that describes how the - * individual shape functions couple with - * each other. - */ - template - static - Table<2,Coupling> - dof_couplings_from_component_couplings (const FiniteElement &fe, - const Table<2,Coupling> &component_couplings); - - /** - * Same function as above for a - * collection of finite elements, - * returning a collection of tables. - * - * The function currently treats - * DoFTools::Couplings::nonzero the same - * as DoFTools::Couplings::always . - */ - template - static - std::vector > - dof_couplings_from_component_couplings (const hp::FECollection &fe, - const Table<2,Coupling> &component_couplings); + /** + * Extract a vector that represents the + * constant modes of the DoFHandler for + * the components chosen by + * component_select. The + * constant modes on a discretization are + * the null space of a Laplace operator + * on the selected components with + * Neumann boundary conditions + * applied. The null space is a necessary + * ingredient for obtaining a good AMG + * preconditioner when using the class + * TrilinosWrappers::PreconditionAMG. + * Since the ML AMG package only works on + * algebraic properties of the respective + * matrix, it has no chance to detect + * whether the matrix comes from a scalar + * or a vector valued problem. However, a + * near null space supplies exactly the + * needed information about these + * components. The null space will + * consist of as many vectors as there + * are true arguments in + * component_select, each of + * which will be one in one vector component and + * zero in all others. We store this + * object in a vector of vectors, where + * the outer vector is of the size of the + * number of selected components, and + * each inner vector has as many + * components as there are (locally owned) degrees of + * freedom in the selected + * components. Note that any matrix + * associated with this null space must + * have been constructed using the same + * component_select argument, + * since the numbering of DoFs is done + * relative to the selected dofs, not to + * all dofs. + * + * The main reason for this + * program is the use of the + * null space with the + * AMG preconditioner. + */ + template + void + extract_constant_modes (const DH &dof_handler, + const std::vector &component_select, + std::vector > &constant_modes); + + /** + * For each active cell of a DoFHandler + * or hp::DoFHandler, extract the active + * finite element index and fill the + * vector given as second argument. This + * vector is assumed to have as many + * entries as there are active cells. + * + * For non-hp DoFHandler objects given as + * first argument, the returned vector + * will consist of only zeros, indicating + * that all cells use the same finite + * element. For a hp::DoFHandler, the + * values may be different, though. + */ + template + void + get_active_fe_indices (const DH &dof_handler, + std::vector &active_fe_indices); + + /** + * 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. + * + * Note that degrees of freedom + * associated with faces, edges, and + * vertices may be associated with + * multiple subdomains if they are + * sitting on partition boundaries. In + * these cases, we put them into one of + * the associated partitions in an + * undefined way. This may sometimes lead + * to different numbers of degrees of + * freedom in partitions, even if the + * number of cells is perfectly + * equidistributed. While this is + * regrettable, it is not a problem in + * practice since the number of degrees + * of freedom on partition boundaries is + * asymptotically vanishing as we refine + * the mesh as long as the number of + * partitions is kept constant. + * + * This function returns the association + * of each DoF with one subdomain. If you + * are looking for the association of + * each @em cell with a subdomain, either + * query the + * cell->subdomain_id() + * function, or use the + * GridTools::get_subdomain_association + * function. + * + * Note that this function is of + * questionable use for DoFHandler objects built on + * parallel::distributed::Triangulation + * since in that case ownership of + * individual degrees of freedom by MPI + * processes is controlled by the DoF + * handler object, not based on some + * geometric algorithm in conjunction + * with subdomain id. In particular, the + * degrees of freedom identified by the + * functions in this namespace as + * associated with a subdomain are not + * the same the + * DoFHandler class + * identifies as those it owns. + */ + template + void + get_subdomain_association (const DH &dof_handler, + std::vector &subdomain); + + /** + * Count how many degrees of freedom are + * uniquely associated with the given + * @p subdomain index. + * + * Note that there may be rare cases + * where cells with the given @p + * subdomain index exist, but none of its + * degrees of freedom are actually + * associated with it. In that case, the + * returned value will be zero. + * + * This function will generate an + * exception if there are no cells with + * the given @p subdomain index. + * + * This function returns the number of + * DoFs associated with one subdomain. If + * you are looking for the association of + * @em cells with this subdomain, use the + * GridTools::count_cells_with_subdomain_association + * function. + * + * Note that this function is of + * questionable use for DoFHandler objects built on + * parallel::distributed::Triangulation + * since in that case ownership of + * individual degrees of freedom by MPI + * processes is controlled by the DoF + * handler object, not based on some + * geometric algorithm in conjunction + * with subdomain id. In particular, the + * degrees of freedom identified by the + * functions in this namespace as + * associated with a subdomain are not + * the same the + * DoFHandler class + * identifies as those it owns. + */ + template + unsigned int + count_dofs_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain); + + /** + * Count how many degrees of freedom are + * uniquely associated with the given + * @p subdomain index. + * + * This function does what the previous + * one does except that it splits the + * result among the vector components of + * the finite element in use by the + * DoFHandler object. The last argument + * (which must have a length equal to the + * number of vector components) will + * therefore store how many degrees of + * freedom of each vector component are + * associated with the given subdomain. + * + * Note that this function is of + * questionable use for DoFHandler objects built on + * parallel::distributed::Triangulation + * since in that case ownership of + * individual degrees of freedom by MPI + * processes is controlled by the DoF + * handler object, not based on some + * geometric algorithm in conjunction + * with subdomain id. In particular, the + * degrees of freedom identified by the + * functions in this namespace as + * associated with a subdomain are not + * the same the + * DoFHandler class + * identifies as those it owns. + */ + template + void + count_dofs_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain, + std::vector &n_dofs_on_subdomain); + + /** + * Return a set of indices that denotes + * the degrees of freedom that live on + * the given subdomain, i.e. that are on + * cells owned by the current + * processor. Note that this includes the + * ones that this subdomain "owns" + * (i.e. the ones for which + * get_subdomain_association() returns a + * value equal to the subdomain given + * here and that are selected by the + * extract_locally_owned() function) but + * also all of those that sit on the + * boundary between the given subdomain + * and other subdomain. In essence, + * degrees of freedom that sit on + * boundaries between subdomain will be + * in the index sets returned by this + * function for more than one subdomain. + * + * Note that this function is of + * questionable use for DoFHandler objects built on + * parallel::distributed::Triangulation + * since in that case ownership of + * individual degrees of freedom by MPI + * processes is controlled by the DoF + * handler object, not based on some + * geometric algorithm in conjunction + * with subdomain id. In particular, the + * degrees of freedom identified by the + * functions in this namespace as + * associated with a subdomain are not + * the same the + * DoFHandler class + * identifies as those it owns. + */ + template + IndexSet + dof_indices_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain); + + /** + * Count how many degrees of + * freedom out of the total + * number belong to each + * component. If the number of + * components the finite element + * has is one (i.e. you only have + * one scalar variable), then the + * number in this component + * obviously equals the total + * number of degrees of + * freedom. Otherwise, the sum of + * the DoFs in all the components + * needs to equal the total + * number. + * + * However, the last statement + * does not hold true if the + * finite element is not + * primitive, i.e. some or all of + * its shape functions are + * non-zero in more than one + * vector component. This + * applies, for example, to the + * Nedelec or Raviart-Thomas + * elements. In this case, a + * degree of freedom is counted + * in each component in which it + * is non-zero, so that the sum + * mentioned above is greater + * than the total number of + * degrees of freedom. + * + * This behavior can be switched + * off by the optional parameter + * vector_valued_once. If + * this is true, the + * number of components of a + * nonprimitive vector valued + * element is collected only in + * the first component. All other + * components will have a count + * of zero. + * + * The additional optional + * argument @p target_component + * allows for a re-sorting and + * grouping of components. To + * this end, it contains for each + * component the component number + * it shall be counted as. Having + * the same number entered + * several times sums up several + * components as the same. One of + * the applications of this + * argument is when you want to + * form block matrices and + * vectors, but want to pack + * several components into the + * same block (for example, when + * you have @p dim velocities + * and one pressure, to put all + * velocities into one block, and + * the pressure into another). + * + * The result is returned in @p + * dofs_per_component. Note that + * the size of @p + * dofs_per_component needs to be + * enough to hold all the indices + * specified in @p + * target_component. If this is + * not the case, an assertion is + * thrown. The indices not + * targetted by target_components + * are left untouched. + */ + template + void + count_dofs_per_component (const DoFHandler& dof_handler, + std::vector& dofs_per_component, + const bool vector_valued_once = false, + std::vector target_component + = std::vector()); + + /** + * Count the degrees of freedom + * in each block. This function + * is similar to + * count_dofs_per_component(), + * with the difference that the + * counting is done by + * blocks. See @ref GlossBlock + * "blocks" in the glossary for + * details. Again the vectors are + * assumed to have the correct + * size before calling this + * function. If this is not the + * case, an assertion is thrown. + * + * This function is used in the + * step-22, + * step-31, and + * step-32 tutorial + * programs. + * + * @pre The dofs_per_block + * variable has as many + * components as the finite + * element used by the + * dof_handler argument has + * blocks, or alternatively as + * many blocks as are enumerated + * in the target_blocks argument + * if given. + */ + template + void + count_dofs_per_block (const DoFHandler& dof_handler, + std::vector& dofs_per_block, + std::vector target_blocks + = std::vector()); + + /** + * @deprecated See the previous + * function with the same name + * for a description. This + * function exists for + * compatibility with older + * versions only. + */ + template + void + count_dofs_per_component (const DoFHandler& dof_handler, + std::vector& dofs_per_component, + std::vector target_component); + + /** + * This function can be used when + * different variables shall be + * discretized on different + * grids, where one grid is + * coarser than the other. This + * idea might seem nonsensical at + * first, but has reasonable + * applications in inverse + * (parameter estimation) + * problems, where there might + * not be enough information to + * recover the parameter on the + * same grid as the state + * variable; furthermore, the + * smoothness properties of state + * variable and parameter might + * not be too much related, so + * using different grids might be + * an alternative to using + * stronger regularization of the + * problem. + * + * The basic idea of this + * function is explained in the + * following. Let us, for + * convenience, denote by + * ``parameter grid'' the coarser + * of the two grids, and by + * ``state grid'' the finer of + * the two. We furthermore assume + * that the finer grid can be + * obtained by refinement of the + * coarser one, i.e. the fine + * grid is at least as much + * refined as the coarse grid at + * each point of the + * domain. Then, each shape + * function on the coarse grid + * can be represented as a linear + * combination of shape functions + * on the fine grid (assuming + * identical ansatz + * spaces). Thus, if we + * discretize as usual, using + * shape functions on the fine + * grid, we can consider the + * restriction that the parameter + * variable shall in fact be + * discretized by shape functions + * on the coarse grid as a + * constraint. These constraints + * are linear and happen to have + * the form managed by the + * ``ConstraintMatrix'' class. + * + * The construction of these + * constraints is done as + * follows: for each of the + * degrees of freedom (i.e. shape + * functions) on the coarse grid, + * we compute its representation + * on the fine grid, i.e. how the + * linear combination of shape + * functions on the fine grid + * looks like that resembles the + * shape function on the coarse + * grid. From this information, + * we can then compute the + * constraints which have to hold + * if a solution of a linear + * equation on the fine grid + * shall be representable on the + * coarse grid. The exact + * algorithm how these + * constraints can be computed is + * rather complicated and is best + * understood by reading the + * source code, which contains + * many comments. + * + * Before explaining the use of + * this function, we would like + * to state that the total number + * of degrees of freedom used for + * the discretization is not + * reduced by the use of this + * function, i.e. even though we + * discretize one variable on a + * coarser grid, the total number + * of degrees of freedom is that + * of the fine grid. This seems + * to be counter-productive, + * since it does not give us a + * benefit from using a coarser + * grid. The reason why it may be + * useful to choose this approach + * nonetheless is three-fold: + * first, as stated above, there + * might not be enough + * information to recover a + * parameter on a fine grid, + * i.e. we chose to discretize it + * on the coarse grid not to save + * DoFs, but for other + * reasons. Second, the + * ``ConstraintMatrix'' includes + * the constraints into the + * linear system of equations, by + * which constrained nodes become + * dummy nodes; we may therefore + * exclude them from the linear + * algebra, for example by + * sorting them to the back of + * the DoF numbers and simply + * calling the solver for the + * upper left block of the matrix + * which works on the + * non-constrained nodes only, + * thus actually realizing the + * savings in numerical effort + * from the reduced number of + * actual degrees of freedom. The + * third reason is that for some + * or other reason we have chosen + * to use two different grids, it + * may be actually quite + * difficult to write a function + * that assembles the system + * matrix for finite element + * spaces on different grids; + * using the approach of + * constraints as with this + * function allows to use + * standard techniques when + * discretizing on only one grid + * (the finer one) without having + * to take care of the fact that + * one or several of the variable + * actually belong to different + * grids. + * + * The use of this function is as + * follows: it accepts as + * parameters two DoF Handlers, + * the first of which refers to + * the coarse grid and the second + * of which is the fine grid. On + * both, a finite element is + * represented by the DoF handler + * objects, which will usually + * have several components, which + * may belong to different finite + * elements. The second and + * fourth parameter of this + * function therefore state which + * variable on the coarse grid + * shall be used to restrict the + * stated component on the fine + * grid. Of course, the finite + * elements used for the + * respective components on the + * two grids need to be the + * same. An example may clarify + * this: consider the parameter + * estimation mentioned briefly + * above; there, on the fine grid + * the whole discretization is + * done, thus the variables are + * ``u'', ``q'', and the Lagrange + * multiplier ``lambda'', which + * are discretized using + * continuous linear, piecewise + * constant discontinuous, and + * continuous linear elements, + * respectively. Only the + * parameter ``q'' shall be + * represented on the coarse + * grid, thus the DoFHandler + * object on the coarse grid + * represents only one variable, + * discretized using piecewise + * constant discontinuous + * elements. Then, the parameter + * denoting the component on the + * coarse grid would be zero (the + * only possible choice, since + * the variable on the coarse + * grid is scalar), and one on + * the fine grid (corresponding + * to the variable ``q''; zero + * would be ``u'', two would be + * ``lambda''). Furthermore, an + * object of type IntergridMap + * is needed; this could in + * principle be generated by the + * function itself from the two + * DoFHandler objects, but since + * it is probably available + * anyway in programs that use + * this function, we shall use it + * instead of re-generating + * it. Finally, the computed + * constraints are entered into a + * variable of type + * ConstraintMatrix; the + * constraints are added, + * i.e. previous contents which + * may have, for example, be + * obtained from hanging nodes, + * are not deleted, so that you + * only need one object of this + * type. + */ + template + void + compute_intergrid_constraints (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + ConstraintMatrix &constraints); + + + /** + * This function generates a + * matrix such that when a vector + * of data with as many elements + * as there are degrees of + * freedom of this component on + * the coarse grid is multiplied + * to this matrix, we obtain a + * vector with as many elements + * are there are global degrees + * of freedom on the fine + * grid. All the elements of the + * other components of the finite + * element fields on the fine + * grid are not touched. + * + * The output of this function is + * a compressed format that can + * be given to the @p reinit + * functions of the + * SparsityPattern ad + * SparseMatrix classes. + */ + template + void + compute_intergrid_transfer_representation (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + std::vector > &transfer_representation); + + /** + * Create a mapping from degree + * of freedom indices to the + * index of that degree of + * freedom on the boundary. After + * this operation, mapping[dof] + * gives the index of the + * degree of freedom with global + * number @p dof in the list of + * degrees of freedom on the + * boundary. If the degree of + * freedom requested is not on + * the boundary, the value of + * mapping[dof] is + * @p invalid_dof_index. This + * function is mainly used when + * setting up matrices and + * vectors on the boundary from + * the trial functions, which + * have global numbers, while the + * matrices and vectors use + * numbers of the trial functions + * local to the boundary. + * + * Prior content of @p mapping + * is deleted. + */ + template + void + map_dof_to_boundary_indices (const DH &dof_handler, + std::vector &mapping); + + /** + * Same as the previous function, + * except that only those parts + * of the boundary are considered + * for which the boundary + * indicator is listed in the + * second argument. + * + * See the general doc of this + * class for more information. + */ + template + void + map_dof_to_boundary_indices (const DH &dof_handler, + const std::set &boundary_indicators, + std::vector &mapping); + + /** + * Return a list of support + * points for all the degrees of + * freedom handled by this DoF + * handler object. This function, + * of course, only works if the + * finite element object used by + * the DoF handler object + * actually provides support + * points, i.e. no edge elements + * or the like. Otherwise, an + * exception is thrown. + * + * The given array must have a + * length of as many elements as + * there are degrees of freedom. + */ + template + void + map_dofs_to_support_points (const Mapping &mapping, + const DoFHandler &dof_handler, + std::vector > &support_points); + + /** + * This is the opposite function + * to the one above. It generates + * a map where the keys are the + * support points of the degrees + * of freedom, while the values + * are the DoF indices. + * + * Since there is no natural + * order in the space of points + * (except for the 1d case), you + * have to provide a map with an + * explicitly specified + * comparator object. This + * function is therefore + * templatized on the comparator + * object. Previous content of + * the map object is deleted in + * this function. + * + * Just as with the function + * above, it is assumed that the + * finite element in use here + * actually supports the notion + * of support points of all its + * components. + */ + template + void + map_support_points_to_dofs (const Mapping &mapping, + const DH &dof_handler, + std::map, unsigned int, Comp> &point_to_index_map); + + /** + * Map a coupling table from the + * user friendly organization by + * components to the organization + * by blocks. Specializations of + * this function for DoFHandler + * and hp::DoFHandler are + * required due to the different + * results of their finite + * element access. + * + * The return vector will be + * initialized to the correct + * length inside this function. + */ + template + void + convert_couplings_to_blocks (const hp::DoFHandler& dof_handler, + const Table<2, Coupling>& table_by_component, + std::vector >& tables_by_block); + + /** + * Make a constraint matrix for the + * constraints that result from zero + * boundary values. + * + * This function constrains all + * degrees of freedom on the + * boundary. Optionally, you can + * add a component mask, which + * restricts this functionality + * to a subset of an FESystem. + * + * For non-@ref GlossPrimitive "primitive" + * shape functions, any degree of freedom + * is affected that belongs to a + * shape function where at least + * one of its nonzero components + * is affected. + * + * This function is used + * in step-36, for + * example. + * + * @ingroup constraints + */ + template class DH> + void + make_zero_boundary_constraints (const DH &dof, + ConstraintMatrix &zero_boundary_constraints, + const std::vector &component_mask_=std::vector()); + + /** + * Map a coupling table from the + * user friendly organization by + * components to the organization + * by blocks. Specializations of + * this function for DoFHandler + * and hp::DoFHandler are + * required due to the different + * results of their finite + * element access. + * + * The return vector will be + * initialized to the correct + * length inside this function. + */ + template + void + convert_couplings_to_blocks (const DoFHandler& dof_handler, + const Table<2, Coupling>& table_by_component, + std::vector >& tables_by_block); + + /** + * Given a finite element and a table how + * the vector components of it couple + * with each other, compute and return a + * table that describes how the + * individual shape functions couple with + * each other. + */ + template + Table<2,Coupling> + dof_couplings_from_component_couplings (const FiniteElement &fe, + const Table<2,Coupling> &component_couplings); + + /** + * Same function as above for a + * collection of finite elements, + * returning a collection of tables. + * + * The function currently treats + * DoFTools::Couplings::nonzero the same + * as DoFTools::Couplings::always . + */ + template + std::vector > + dof_couplings_from_component_couplings (const hp::FECollection &fe, + const Table<2,Coupling> &component_couplings); - /** - * Exception - */ - DeclException0 (ExcFEHasNoSupportPoints); - /** - * Exception - */ - DeclException0 (ExcFENotPrimitive); - /** - * Exception - */ - DeclException2 (ExcWrongSize, - int, int, - << "The dimension " << arg1 << " of the vector is wrong. " - << "It should be " << arg2); - /** - * Exception - */ - DeclException2 (ExcInvalidComponent, - int, int, - << "The component you gave (" << arg1 << ") " - << "is invalid with respect to the number " - << "of components in the finite element " - << "(" << arg2 << ")"); - /** - * Exception - */ - DeclException0 (ExcFiniteElementsDontMatch); - /** - * Exception - */ - DeclException0 (ExcGridNotCoarser); - /** - * Exception - */ - DeclException0 (ExcGridsDontMatch); - /** - * Exception - */ - DeclException0 (ExcNoFESelected); - /** - * Exception - */ - DeclException0 (ExcInvalidBoundaryIndicator); -}; + /** + * Exception + */ + DeclException0 (ExcFEHasNoSupportPoints); + /** + * Exception + */ + DeclException0 (ExcFENotPrimitive); + /** + * Exception + */ + DeclException2 (ExcWrongSize, + int, int, + << "The dimension " << arg1 << " of the vector is wrong. " + << "It should be " << arg2); + /** + * Exception + */ + DeclException2 (ExcInvalidComponent, + int, int, + << "The component you gave (" << arg1 << ") " + << "is invalid with respect to the number " + << "of components in the finite element " + << "(" << arg2 << ")"); + /** + * Exception + */ + DeclException0 (ExcFiniteElementsDontMatch); + /** + * Exception + */ + DeclException0 (ExcGridNotCoarser); + /** + * Exception + */ + DeclException0 (ExcGridsDontMatch); + /** + * Exception + */ + DeclException0 (ExcNoFESelected); + /** + * Exception + */ + DeclException0 (ExcInvalidBoundaryIndicator); +} /* ------------------------- inline functions -------------- */ +#ifndef DOXYGEN +namespace DoFTools +{ /** * Operator computing the maximum coupling out of two. * * @relates DoFTools */ -inline -DoFTools::Coupling operator |= (DoFTools::Coupling& c1, - const DoFTools::Coupling c2) -{ - if (c2 == DoFTools::always) - c1 = DoFTools::always; - else if (c1 != DoFTools::always && c2 == DoFTools::nonzero) - return c1 = DoFTools::nonzero; - return c1; -} + inline + Coupling operator |= (Coupling& c1, + const Coupling c2) + { + if (c2 == Coupling::always) + c1 = Coupling::always; + else if (c1 != Coupling::always && c2 == Coupling::nonzero) + return c1 = Coupling::nonzero; + return c1; + } /** @@ -2185,151 +2180,161 @@ DoFTools::Coupling operator |= (DoFTools::Coupling& c1, * * @relates DoFTools */ -inline -DoFTools::Coupling operator | (const DoFTools::Coupling c1, - const DoFTools::Coupling c2) -{ - if (c1 == DoFTools::always || c2 == DoFTools::always) - return DoFTools::always; - if (c1 == DoFTools::nonzero || c2 == DoFTools::nonzero) - return DoFTools::nonzero; - return DoFTools::none; -} + inline + Coupling operator | (const Coupling c1, + const Coupling c2) + { + if (c1 == Coupling::always || c2 == Coupling::always) + return Coupling::always; + if (c1 == Coupling::nonzero || c2 == Coupling::nonzero) + return Coupling::nonzero; + return Coupling::none; + } // ---------------------- inline and template functions -------------------- -template -inline unsigned int -DoFTools::max_dofs_per_cell (const DoFHandler &dh) -{ - return dh.get_fe().dofs_per_cell; -} - - -template -inline unsigned int -DoFTools::max_dofs_per_face (const DoFHandler &dh) -{ - return dh.get_fe().dofs_per_face; -} - - -template -inline unsigned int -DoFTools::max_dofs_per_vertex (const DoFHandler &dh) -{ - return dh.get_fe().dofs_per_vertex; -} - - -template -inline unsigned int -DoFTools::n_components (const DoFHandler &dh) -{ - return dh.get_fe().n_components(); -} - - - -template -inline -bool -DoFTools::fe_is_primitive (const DoFHandler &dh) -{ - return dh.get_fe().is_primitive(); -} - - -template -inline unsigned int -DoFTools::max_dofs_per_cell (const hp::DoFHandler &dh) -{ - return dh.get_fe().max_dofs_per_cell (); -} - - -template -inline unsigned int -DoFTools::max_dofs_per_face (const hp::DoFHandler &dh) -{ - return dh.get_fe().max_dofs_per_face (); -} - - -template -inline unsigned int -DoFTools::max_dofs_per_vertex (const hp::DoFHandler &dh) -{ - return dh.get_fe().max_dofs_per_vertex (); -} - - -template -inline unsigned int -DoFTools::n_components (const hp::DoFHandler &dh) -{ - return dh.get_fe()[0].n_components(); -} - - -template -inline -bool -DoFTools::fe_is_primitive (const hp::DoFHandler &dh) -{ - return dh.get_fe()[0].is_primitive(); -} - - -template -inline -void -DoFTools::make_sparsity_pattern (const DH &dof, - const std::vector > &mask, - SparsityPattern &sparsity_pattern) -{ - const unsigned int ncomp = dof.get_fe().n_components(); - - Assert (mask.size() == ncomp, - ExcDimensionMismatch(mask.size(), ncomp)); - for (unsigned int i=0; i couplings(ncomp, ncomp); - for (unsigned int i=0;i -void -DoFTools::map_support_points_to_dofs ( - const Mapping &mapping, - const DH &dof_handler, - std::map, unsigned int, Comp> &point_to_index_map) -{ - // let the checking of arguments be - // done by the function first - // called - std::vector > support_points (dof_handler.n_dofs()); - map_dofs_to_support_points (mapping, dof_handler, support_points); - // now copy over the results of the - // previous function into the - // output arg - point_to_index_map.clear (); - for (unsigned int i=0; i + inline + unsigned int + max_dofs_per_cell (const DoFHandler &dh) + { + return dh.get_fe().dofs_per_cell; + } + + + template + inline + unsigned int + max_dofs_per_face (const DoFHandler &dh) + { + return dh.get_fe().dofs_per_face; + } + + + template + inline + unsigned int + max_dofs_per_vertex (const DoFHandler &dh) + { + return dh.get_fe().dofs_per_vertex; + } + + + template + inline + unsigned int + n_components (const DoFHandler &dh) + { + return dh.get_fe().n_components(); + } + + + + template + inline + bool + fe_is_primitive (const DoFHandler &dh) + { + return dh.get_fe().is_primitive(); + } + + + template + inline + unsigned int + max_dofs_per_cell (const hp::DoFHandler &dh) + { + return dh.get_fe().max_dofs_per_cell (); + } + + + template + inline + unsigned int + max_dofs_per_face (const hp::DoFHandler &dh) + { + return dh.get_fe().max_dofs_per_face (); + } + + + template + inline + unsigned int + max_dofs_per_vertex (const hp::DoFHandler &dh) + { + return dh.get_fe().max_dofs_per_vertex (); + } + + + template + inline + unsigned int + n_components (const hp::DoFHandler &dh) + { + return dh.get_fe()[0].n_components(); + } + + + template + inline + bool + fe_is_primitive (const hp::DoFHandler &dh) + { + return dh.get_fe()[0].is_primitive(); + } + + + template + inline + void + make_sparsity_pattern (const DH &dof, + const std::vector > &mask, + SparsityPattern &sparsity_pattern) + { + const unsigned int ncomp = dof.get_fe().n_components(); + + Assert (mask.size() == ncomp, + ExcDimensionMismatch(mask.size(), ncomp)); + for (unsigned int i=0; i couplings(ncomp, ncomp); + for (unsigned int i=0;i + void + map_support_points_to_dofs ( + const Mapping &mapping, + const DH &dof_handler, + std::map, unsigned int, Comp> &point_to_index_map) + { + // let the checking of arguments be + // done by the function first + // called + std::vector > support_points (dof_handler.n_dofs()); + map_dofs_to_support_points (mapping, dof_handler, support_points); + // now copy over the results of the + // previous function into the + // output arg + point_to_index_map.clear (); + for (unsigned int i=0; i -void -DoFTools::make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) +namespace DoFTools { - const unsigned int n_dofs = dof.n_dofs(); - - Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); - Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); - - std::vector dofs_on_this_cell; - dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - - // In case we work with a distributed - // sparsity pattern of Trilinos type, we - // only have to do the work if the - // current cell is owned by the calling - // processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if (((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) - && - !cell->is_artificial() - && - !cell->is_ghost()) - { - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - dofs_on_this_cell.resize (dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); - - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs); - } -} + + template + void + make_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id_t subdomain_id) + { + const unsigned int n_dofs = dof.n_dofs(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (sparsity.n_cols() == n_dofs, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + + std::vector dofs_on_this_cell; + dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // In case we work with a distributed + // sparsity pattern of Trilinos type, we + // only have to do the work if the + // current cell is owned by the calling + // processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if (((subdomain_id == types::invalid_subdomain_id) + || + (subdomain_id == cell->subdomain_id())) + && + !cell->is_artificial() + && + !cell->is_ghost()) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + dofs_on_this_cell.resize (dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this + // cell. if no constraints pattern was + // given, then the following call acts + // as if simply no constraints existed + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs); + } + } -template -void -DoFTools::make_sparsity_pattern (const DH &dof, - const Table<2,Coupling> &couplings, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) -{ - const unsigned int n_dofs = dof.n_dofs(); - - Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); - Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); - Assert (couplings.n_rows() == dof.get_fe().n_components(), - ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components())); - Assert (couplings.n_cols() == dof.get_fe().n_components(), - ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); - - const hp::FECollection fe_collection (dof.get_fe()); - - // first, for each finite element, build a - // mask for each dof, not like the one - // given which represents components. make - // sure we do the right thing also with - // respect to non-primitive shape - // functions, which takes some additional - // thought - std::vector > dof_mask(fe_collection.size()); - - // check whether the table of couplings - // contains only true arguments, i.e., we - // do not exclude any index. that is the - // easy case, since we don't have to set - // up the tables - bool need_dof_mask = false; - for (unsigned int i=0; i + void + make_sparsity_pattern (const DH &dof, + const Table<2,Coupling> &couplings, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id_t subdomain_id) + { + const unsigned int n_dofs = dof.n_dofs(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (sparsity.n_cols() == n_dofs, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + Assert (couplings.n_rows() == dof.get_fe().n_components(), + ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components())); + Assert (couplings.n_cols() == dof.get_fe().n_components(), + ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); + + const hp::FECollection fe_collection (dof.get_fe()); + + // first, for each finite element, build a + // mask for each dof, not like the one + // given which represents components. make + // sure we do the right thing also with + // respect to non-primitive shape + // functions, which takes some additional + // thought + std::vector > dof_mask(fe_collection.size()); + + // check whether the table of couplings + // contains only true arguments, i.e., we + // do not exclude any index. that is the + // easy case, since we don't have to set + // up the tables + bool need_dof_mask = false; + for (unsigned int i=0; i dofs_on_this_cell(fe_collection.max_dofs_per_cell()); - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - - // In case we work with a distributed - // sparsity pattern of Trilinos type, we - // only have to do the work if the - // current cell is owned by the calling - // processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if (((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) - && - !cell->is_artificial() - && - !cell->is_ghost()) - { - const unsigned int fe_index = cell->active_fe_index(); - const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell; - dofs_on_this_cell.resize (dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); + std::vector dofs_on_this_cell(fe_collection.max_dofs_per_cell()); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // In case we work with a distributed + // sparsity pattern of Trilinos type, we + // only have to do the work if the + // current cell is owned by the calling + // processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if (((subdomain_id == types::invalid_subdomain_id) + || + (subdomain_id == cell->subdomain_id())) + && + !cell->is_artificial() + && + !cell->is_ghost()) + { + const unsigned int fe_index = cell->active_fe_index(); + const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell; + dofs_on_this_cell.resize (dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs, - dof_mask[fe_index]); - } -} + // make sparsity pattern for this + // cell. if no constraints pattern was + // given, then the following call acts + // as if simply no constraints existed + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs, + dof_mask[fe_index]); + } + } -template -void -DoFTools::make_sparsity_pattern ( - const DH &dof_row, - const DH &dof_col, - SparsityPattern &sparsity) -{ - const unsigned int n_dofs_row = dof_row.n_dofs(); - const unsigned int n_dofs_col = dof_col.n_dofs(); - Assert (sparsity.n_rows() == n_dofs_row, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row)); - Assert (sparsity.n_cols() == n_dofs_col, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col)); + template + void + make_sparsity_pattern ( + const DH &dof_row, + const DH &dof_col, + SparsityPattern &sparsity) + { + const unsigned int n_dofs_row = dof_row.n_dofs(); + const unsigned int n_dofs_col = dof_col.n_dofs(); + Assert (sparsity.n_rows() == n_dofs_row, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row)); + Assert (sparsity.n_cols() == n_dofs_col, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col)); - const std::list > - cell_list - = GridTools::get_finest_common_cells (dof_row, dof_col); + const std::list > + cell_list + = GridTools::get_finest_common_cells (dof_row, dof_col); - typename std::list > - ::const_iterator - cell_iter = cell_list.begin(); - for (; cell_iter!=cell_list.end(); ++cell_iter) - { - const typename DH::cell_iterator cell_row = cell_iter->first; - const typename DH::cell_iterator cell_col = cell_iter->second; - - if (!cell_row->has_children() && !cell_col->has_children()) - { - const unsigned int dofs_per_cell_row = - cell_row->get_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; ihas_children()) - { - const std::vector - child_cells = GridTools::get_active_child_cells (cell_row); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row_child->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; i - child_cells = GridTools::get_active_child_cells (cell_col); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col_child->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col_child->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; i > + ::const_iterator + cell_iter = cell_list.begin(); + for (; cell_iter!=cell_list.end(); ++cell_iter) + { + const typename DH::cell_iterator cell_row = cell_iter->first; + const typename DH::cell_iterator cell_col = cell_iter->second; + if (!cell_row->has_children() && !cell_col->has_children()) + { + const unsigned int dofs_per_cell_row = + cell_row->get_fe().dofs_per_cell; + const unsigned int dofs_per_cell_col = + cell_col->get_fe().dofs_per_cell; + std::vector + local_dof_indices_row(dofs_per_cell_row); + std::vector + local_dof_indices_col(dofs_per_cell_col); + cell_row->get_dof_indices (local_dof_indices_row); + cell_col->get_dof_indices (local_dof_indices_col); + for (unsigned int i=0; ihas_children()) + { + const std::vector + child_cells = GridTools::get_active_child_cells (cell_row); + for (unsigned int i=0; iget_fe().dofs_per_cell; + const unsigned int dofs_per_cell_col = + cell_col->get_fe().dofs_per_cell; + std::vector + local_dof_indices_row(dofs_per_cell_row); + std::vector + local_dof_indices_col(dofs_per_cell_col); + cell_row_child->get_dof_indices (local_dof_indices_row); + cell_col->get_dof_indices (local_dof_indices_col); + for (unsigned int i=0; i + child_cells = GridTools::get_active_child_cells (cell_col); + for (unsigned int i=0; iget_fe().dofs_per_cell; + const unsigned int dofs_per_cell_col = + cell_col_child->get_fe().dofs_per_cell; + std::vector + local_dof_indices_row(dofs_per_cell_row); + std::vector + local_dof_indices_col(dofs_per_cell_col); + cell_row->get_dof_indices (local_dof_indices_row); + cell_col_child->get_dof_indices (local_dof_indices_col); + for (unsigned int i=0; i -void -DoFTools::make_boundary_sparsity_pattern ( - const DH &dof, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity) -{ - if (DH::dimension == 1) - { - // there are only 2 boundary - // indicators in 1d, so it is no - // performance problem to call the - // other function - typename DH::FunctionMap boundary_indicators; - boundary_indicators[0] = 0; - boundary_indicators[1] = 0; - make_boundary_sparsity_pattern (dof, - boundary_indicators, - dof_to_boundary_mapping, - sparsity); - return; - } - const unsigned int n_dofs = dof.n_dofs(); - AssertDimension (dof_to_boundary_mapping.size(), n_dofs); - AssertDimension (sparsity.n_rows(), dof.n_boundary_dofs()); - AssertDimension (sparsity.n_cols(), dof.n_boundary_dofs()); + template + void + make_boundary_sparsity_pattern ( + const DH &dof, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &sparsity) + { + if (DH::dimension == 1) + { + // there are only 2 boundary + // indicators in 1d, so it is no + // performance problem to call the + // other function + typename DH::FunctionMap boundary_indicators; + boundary_indicators[0] = 0; + boundary_indicators[1] = 0; + make_boundary_sparsity_pattern (dof, + boundary_indicators, + dof_to_boundary_mapping, + sparsity); + return; + } + + const unsigned int n_dofs = dof.n_dofs(); + + AssertDimension (dof_to_boundary_mapping.size(), n_dofs); + AssertDimension (sparsity.n_rows(), dof.n_boundary_dofs()); + AssertDimension (sparsity.n_cols(), dof.n_boundary_dofs()); #ifdef DEBUG - if (sparsity.n_rows() != 0) - { - unsigned int max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); - }; + if (sparsity.n_rows() != 0) + { + unsigned int max_element = 0; + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != DH::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + AssertDimension (max_element, sparsity.n_rows()-1); + }; #endif - std::vector dofs_on_this_face; - dofs_on_this_face.reserve (max_dofs_per_face(dof)); - - // loop over all faces to check - // whether they are at a - // boundary. note that we need not - // take special care of single - // lines (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every boundary line is - // also part of a boundary face. - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->at_boundary(f)) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_this_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_this_face, - cell->active_fe_index()); - - // make sparsity pattern for this cell - for (unsigned int i=0; i dofs_on_this_face; + dofs_on_this_face.reserve (max_dofs_per_face(dof)); + + // loop over all faces to check + // whether they are at a + // boundary. note that we need not + // take special care of single + // lines (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every boundary line is + // also part of a boundary face. + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_this_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_this_face, + cell->active_fe_index()); + + // make sparsity pattern for this cell + for (unsigned int i=0; i -void DoFTools::make_boundary_sparsity_pattern ( - const DH &dof, - const typename FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity) -{ - if (DH::dimension == 1) - { - // first check left, then right - // boundary point - for (unsigned int direction=0; direction<2; ++direction) - { - // if this boundary is not - // requested, then go on with next one - if (boundary_indicators.find(direction) == - boundary_indicators.end()) - continue; - - // find active cell at that - // boundary: first go to - // left/right, then to children - typename DH::cell_iterator cell = dof.begin(0); - while (!cell->at_boundary(direction)) - cell = cell->neighbor(direction); - while (!cell->active()) - cell = cell->child(direction); - - const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; - std::vector boundary_dof_boundary_indices (dofs_per_vertex); - - // next get boundary mapped dof - // indices of boundary dofs - for (unsigned int i=0; ivertex_dof_index(direction,i)]; - - for (unsigned int i=0; i + void make_boundary_sparsity_pattern ( + const DH &dof, + const typename FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &sparsity) + { + if (DH::dimension == 1) + { + // first check left, then right + // boundary point + for (unsigned int direction=0; direction<2; ++direction) + { + // if this boundary is not + // requested, then go on with next one + if (boundary_indicators.find(direction) == + boundary_indicators.end()) + continue; + + // find active cell at that + // boundary: first go to + // left/right, then to children + typename DH::cell_iterator cell = dof.begin(0); + while (!cell->at_boundary(direction)) + cell = cell->neighbor(direction); + while (!cell->active()) + cell = cell->child(direction); + + const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; + std::vector boundary_dof_boundary_indices (dofs_per_vertex); + + // next get boundary mapped dof + // indices of boundary dofs + for (unsigned int i=0; ivertex_dof_index(direction,i)]; + + for (unsigned int i=0; i::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); - }; + if (sparsity.n_rows() != 0) + { + unsigned int max_element = 0; + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != DH::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + AssertDimension (max_element, sparsity.n_rows()-1); + }; #endif - std::vector dofs_on_this_face; - dofs_on_this_face.reserve (max_dofs_per_face(dof)); - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (boundary_indicators.find(cell->face(f)->boundary_indicator()) != - boundary_indicators.end()) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_this_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_this_face, - cell->active_fe_index()); - - // make sparsity pattern for this cell - for (unsigned int i=0; i dofs_on_this_face; + dofs_on_this_face.reserve (max_dofs_per_face(dof)); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (boundary_indicators.find(cell->face(f)->boundary_indicator()) != + boundary_indicators.end()) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_this_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_this_face, + cell->active_fe_index()); + + // make sparsity pattern for this cell + for (unsigned int i=0; i -void -DoFTools::make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) -{ - const unsigned int n_dofs = dof.n_dofs(); - - AssertDimension (sparsity.n_rows(), n_dofs); - AssertDimension (sparsity.n_cols(), n_dofs); - - std::vector dofs_on_this_cell; - std::vector dofs_on_other_cell; - dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); - dofs_on_other_cell.reserve (max_dofs_per_cell(dof)); - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - - // TODO: in an old implementation, we used - // user flags before to tag faces that were - // already touched. this way, we could reduce - // the work a little bit. now, we instead add - // only data from one side. this should be OK, - // but we need to actually verify it. - - // In case we work with a distributed - // sparsity pattern of Trilinos type, we - // only have to do the work if the - // current cell is owned by the calling - // processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if ((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) - { - const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell; - dofs_on_this_cell.resize (n_dofs_on_this_cell); - cell->get_dof_indices (dofs_on_this_cell); - - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs); - - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; - ++face) + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id_t subdomain_id) + { + const unsigned int n_dofs = dof.n_dofs(); + + AssertDimension (sparsity.n_rows(), n_dofs); + AssertDimension (sparsity.n_cols(), n_dofs); + + std::vector dofs_on_this_cell; + std::vector dofs_on_other_cell; + dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); + dofs_on_other_cell.reserve (max_dofs_per_cell(dof)); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // TODO: in an old implementation, we used + // user flags before to tag faces that were + // already touched. this way, we could reduce + // the work a little bit. now, we instead add + // only data from one side. this should be OK, + // but we need to actually verify it. + + // In case we work with a distributed + // sparsity pattern of Trilinos type, we + // only have to do the work if the + // current cell is owned by the calling + // processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if ((subdomain_id == types::invalid_subdomain_id) + || + (subdomain_id == cell->subdomain_id())) { - typename DH::face_iterator cell_face = cell->face(face); - if (! cell->at_boundary(face) ) - { - typename DH::cell_iterator neighbor = cell->neighbor(face); + const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell; + dofs_on_this_cell.resize (n_dofs_on_this_cell); + cell->get_dof_indices (dofs_on_this_cell); - if (cell_face->has_children()) + // make sparsity pattern for this + // cell. if no constraints pattern was + // given, then the following call acts + // as if simply no constraints existed + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs); + + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + typename DH::face_iterator cell_face = cell->face(face); + if (! cell->at_boundary(face) ) { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->number_of_children(); - ++sub_nr) + typename DH::cell_iterator neighbor = cell->neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->number_of_children(); + ++sub_nr) + { + const typename DH::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + const unsigned int n_dofs_on_neighbor + = sub_neighbor->get_fe().dofs_per_cell; + dofs_on_other_cell.resize (n_dofs_on_neighbor); + sub_neighbor->get_dof_indices (dofs_on_other_cell); + + constraints.add_entries_local_to_global + (dofs_on_this_cell, dofs_on_other_cell, + sparsity, keep_constrained_dofs); + constraints.add_entries_local_to_global + (dofs_on_other_cell, dofs_on_this_cell, + sparsity, keep_constrained_dofs); + } + } + else { - const typename DH::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); + // Refinement edges are taken care of + // by coarser cells + + // TODO: in the distributed case, we miss out + // the constraints when the neighbor cell is + // coarser, but only the current cell is owned + // locally! + if (cell->neighbor_is_coarser(face)) + continue; const unsigned int n_dofs_on_neighbor - = sub_neighbor->get_fe().dofs_per_cell; + = neighbor->get_fe().dofs_per_cell; dofs_on_other_cell.resize (n_dofs_on_neighbor); - sub_neighbor->get_dof_indices (dofs_on_other_cell); + + neighbor->get_dof_indices (dofs_on_other_cell); constraints.add_entries_local_to_global (dofs_on_this_cell, dofs_on_other_cell, sparsity, keep_constrained_dofs); - constraints.add_entries_local_to_global - (dofs_on_other_cell, dofs_on_this_cell, - sparsity, keep_constrained_dofs); + + // only need to add these in case the neighbor + // cell is not locally owned - otherwise, we + // touch each face twice and hence put the + // indices the other way around + if (cell->neighbor(face)->subdomain_id() != + cell->subdomain_id()) + constraints.add_entries_local_to_global + (dofs_on_other_cell, dofs_on_this_cell, + sparsity, keep_constrained_dofs); } } - else - { - // Refinement edges are taken care of - // by coarser cells - - // TODO: in the distributed case, we miss out - // the constraints when the neighbor cell is - // coarser, but only the current cell is owned - // locally! - if (cell->neighbor_is_coarser(face)) - continue; - - const unsigned int n_dofs_on_neighbor - = neighbor->get_fe().dofs_per_cell; - dofs_on_other_cell.resize (n_dofs_on_neighbor); - - neighbor->get_dof_indices (dofs_on_other_cell); - - constraints.add_entries_local_to_global - (dofs_on_this_cell, dofs_on_other_cell, - sparsity, keep_constrained_dofs); - - // only need to add these in case the neighbor - // cell is not locally owned - otherwise, we - // touch each face twice and hence put the - // indices the other way around - if (cell->neighbor(face)->subdomain_id() != - cell->subdomain_id()) - constraints.add_entries_local_to_global - (dofs_on_other_cell, dofs_on_this_cell, - sparsity, keep_constrained_dofs); - } } } - } -} + } -template -void -DoFTools::make_flux_sparsity_pattern ( - const DH &dof, - SparsityPattern &sparsity) -{ - ConstraintMatrix constraints; - make_flux_sparsity_pattern (dof, sparsity, constraints); -} + template + void + make_flux_sparsity_pattern ( + const DH &dof, + SparsityPattern &sparsity) + { + ConstraintMatrix constraints; + make_flux_sparsity_pattern (dof, sparsity, constraints); + } -template -Table<2,DoFTools::Coupling> -DoFTools::dof_couplings_from_component_couplings -(const FiniteElement &fe, - const Table<2,Coupling> &component_couplings) -{ - Assert(component_couplings.n_rows() == fe.n_components(), - ExcDimensionMismatch(component_couplings.n_rows(), - fe.n_components())); - Assert(component_couplings.n_cols() == fe.n_components(), - ExcDimensionMismatch(component_couplings.n_cols(), - fe.n_components())); + template + Table<2,Coupling> + dof_couplings_from_component_couplings + (const FiniteElement &fe, + const Table<2,Coupling> &component_couplings) + { + Assert(component_couplings.n_rows() == fe.n_components(), + ExcDimensionMismatch(component_couplings.n_rows(), + fe.n_components())); + Assert(component_couplings.n_cols() == fe.n_components(), + ExcDimensionMismatch(component_couplings.n_cols(), + fe.n_components())); - const unsigned int n_dofs = fe.dofs_per_cell; + const unsigned int n_dofs = fe.dofs_per_cell; - Table<2,DoFTools::Coupling> dof_couplings (n_dofs, n_dofs); + Table<2,Coupling> dof_couplings (n_dofs, n_dofs); - for (unsigned int i=0; i -std::vector > -DoFTools::dof_couplings_from_component_couplings -(const hp::FECollection &fe, - const Table<2,Coupling> &component_couplings) -{ - std::vector > return_value (fe.size()); - for (unsigned int i=0; i + std::vector > + dof_couplings_from_component_couplings + (const hp::FECollection &fe, + const Table<2,Coupling> &component_couplings) + { + std::vector > return_value (fe.size()); + for (unsigned int i=0; i - void - make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const Table<2,dealii::DoFTools::Coupling> &int_mask, - const Table<2,dealii::DoFTools::Coupling> &flux_mask) + namespace { - const FiniteElement &fe = dof.get_fe(); - - std::vector dofs_on_this_cell(fe.dofs_per_cell); - std::vector dofs_on_other_cell(fe.dofs_per_cell); - - const Table<2,dealii::DoFTools::Coupling> - int_dof_mask = dealii::DoFTools::dof_couplings_from_component_couplings(fe, int_mask), - flux_dof_mask = dealii::DoFTools::dof_couplings_from_component_couplings(fe, flux_mask); - - Table<2,bool> support_on_face(fe.dofs_per_cell, - GeometryInfo::faces_per_cell); - for (unsigned int i=0; i::faces_per_cell;++f) - support_on_face(i,f) = fe.has_support_on_face(i,f); - - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this cell - for (unsigned int i=0; i::faces_per_cell; - ++face) - { - const typename DH::face_iterator - cell_face = cell->face(face); - if (cell_face->user_flag_set ()) - continue; + + // implementation of the same function in + // namespace DoFTools for non-hp + // DoFHandlers + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const Table<2,Coupling> &int_mask, + const Table<2,Coupling> &flux_mask) + { + const FiniteElement &fe = dof.get_fe(); - if (cell->at_boundary (face) ) - { - for (unsigned int i=0; i dofs_on_this_cell(fe.dofs_per_cell); + std::vector dofs_on_other_cell(fe.dofs_per_cell); - if ((flux_dof_mask(i,j) == dealii::DoFTools::always) - || - (flux_dof_mask(i,j) == dealii::DoFTools::nonzero - && - i_non_zero_i - && - j_non_zero_i)) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } - } - } - else - { - typename DH::cell_iterator - neighbor = cell->neighbor(face); - // Refinement edges are taken care of - // by coarser cells - if (cell->neighbor_is_coarser(face)) - continue; + const Table<2,Coupling> + int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask), + flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask); - typename DH::face_iterator cell_face = cell->face(face); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face); + Table<2,bool> support_on_face(fe.dofs_per_cell, + GeometryInfo::faces_per_cell); + for (unsigned int i=0; i::faces_per_cell;++f) + support_on_face(i,f) = fe.has_support_on_face(i,f); - if (cell_face->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename DH::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this cell + for (unsigned int i=0; i::faces_per_cell; + ++face) + { + const typename DH::face_iterator + cell_face = cell->face(face); + if (cell_face->user_flag_set ()) + continue; - sub_neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iat_boundary (face) ) + { + for (unsigned int i=0; ineighbor(face); + // Refinement edges are taken care of + // by coarser cells + if (cell->neighbor_is_coarser(face)) + continue; + + typename DH::face_iterator cell_face = cell->face(face); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->n_children(); + ++sub_nr) + { + const typename DH::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + sub_neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); - } - } - else - { - neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); + } + } + else + { + neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); - } - } - } - } - } - - - // implementation of the same function in - // namespace DoFTools for non-hp - // DoFHandlers - template - void - make_flux_sparsity_pattern (const dealii::hp::DoFHandler &dof, - SparsityPattern &sparsity, - const Table<2,dealii::DoFTools::Coupling> &int_mask, - const Table<2,dealii::DoFTools::Coupling> &flux_mask) - { - // while the implementation above is - // quite optimized and caches a lot of - // data (see e.g. the int/flux_dof_mask - // tables), this is no longer practical - // for the hp version since we would - // have to have it for all combinations - // of elements in the - // hp::FECollection. consequently, the - // implementation here is simpler and - // probably less efficient but at least - // readable... - - const dealii::hp::FECollection &fe = dof.get_fe(); - - std::vector dofs_on_this_cell(dealii::DoFTools::max_dofs_per_cell(dof)); - std::vector dofs_on_other_cell(dealii::DoFTools::max_dofs_per_cell(dof)); - - const std::vector > - int_dof_mask - = dealii::DoFTools::dof_couplings_from_component_couplings(fe, int_mask); - - typename dealii::hp::DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this cell - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - if (int_dof_mask[cell->active_fe_index()](i,j) != dealii::DoFTools::none) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - - // Loop over all interior neighbors - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; - ++face) - { - const typename dealii::hp::DoFHandler::face_iterator - cell_face = cell->face(face); - if (cell_face->user_flag_set ()) - continue; - - if (cell->at_boundary (face) ) - { - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - cell->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - cell->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::nonzero)) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } - else - { - typename dealii::hp::DoFHandler::cell_iterator - neighbor = cell->neighbor(face); - // Refinement edges are taken care of - // by coarser cells - if (cell->neighbor_is_coarser(face)) - continue; - - typename dealii::hp::DoFHandler::face_iterator - cell_face = cell->face(face); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face); - - if (cell_face->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename dealii::hp::DoFHandler::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); - - dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell); - sub_neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - { - for (unsigned int j=0; jget_fe().dofs_per_cell; - ++j) - { - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - sub_neighbor->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - sub_neighbor->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::nonzero)) - { + } + if (flux_dof_mask(i,j) == Coupling::nonzero) + { + if (i_non_zero_i && j_non_zero_e) sparsity.add (dofs_on_this_cell[i], dofs_on_other_cell[j]); + if (i_non_zero_e && j_non_zero_i) sparsity.add (dofs_on_other_cell[i], dofs_on_this_cell[j]); + if (i_non_zero_i && j_non_zero_i) sparsity.add (dofs_on_this_cell[i], dofs_on_this_cell[j]); + if (i_non_zero_e && j_non_zero_e) sparsity.add (dofs_on_other_cell[i], dofs_on_other_cell[j]); - } - - if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == dealii::DoFTools::always) - || - (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == dealii::DoFTools::nonzero)) - { + } + + if (flux_dof_mask(j,i) == Coupling::always) + { + sparsity.add (dofs_on_this_cell[j], + dofs_on_other_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_this_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_other_cell[i]); + } + if (flux_dof_mask(j,i) == Coupling::nonzero) + { + if (j_non_zero_i && i_non_zero_e) sparsity.add (dofs_on_this_cell[j], dofs_on_other_cell[i]); + if (j_non_zero_e && i_non_zero_i) sparsity.add (dofs_on_other_cell[j], dofs_on_this_cell[i]); + if (j_non_zero_i && i_non_zero_i) sparsity.add (dofs_on_this_cell[j], dofs_on_this_cell[i]); + if (j_non_zero_e && i_non_zero_e) sparsity.add (dofs_on_other_cell[j], dofs_on_other_cell[i]); - } - } - } - sub_neighbor->face(neighbor_face)->set_user_flag (); - } - } - else - { - dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell); - neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - { - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - { - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - neighbor->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - neighbor->get_fe().system_to_component_index(j).first) - == dealii::DoFTools::nonzero)) - { - sparsity.add (dofs_on_this_cell[i], - dofs_on_other_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - - if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == dealii::DoFTools::always) - || - (flux_mask(neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == dealii::DoFTools::nonzero)) - { - sparsity.add (dofs_on_this_cell[j], - dofs_on_other_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_this_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - } - } - neighbor->face(neighbor_face)->set_user_flag (); - } - } - } - } + } + } + } + neighbor->face(neighbor_face)->set_user_flag (); + } + } + } + } + } + + + // implementation of the same function in + // namespace DoFTools for non-hp + // DoFHandlers + template + void + make_flux_sparsity_pattern (const dealii::hp::DoFHandler &dof, + SparsityPattern &sparsity, + const Table<2,Coupling> &int_mask, + const Table<2,Coupling> &flux_mask) + { + // while the implementation above is + // quite optimized and caches a lot of + // data (see e.g. the int/flux_dof_mask + // tables), this is no longer practical + // for the hp version since we would + // have to have it for all combinations + // of elements in the + // hp::FECollection. consequently, the + // implementation here is simpler and + // probably less efficient but at least + // readable... + + const dealii::hp::FECollection &fe = dof.get_fe(); + + std::vector dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof)); + std::vector dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof)); + + const std::vector > + int_dof_mask + = dof_couplings_from_component_couplings(fe, int_mask); + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this cell + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + if (int_dof_mask[cell->active_fe_index()](i,j) != Coupling::none) + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + + // Loop over all interior neighbors + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + const typename dealii::hp::DoFHandler::face_iterator + cell_face = cell->face(face); + if (cell_face->user_flag_set ()) + continue; + + if (cell->at_boundary (face) ) + { + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + cell->get_fe().system_to_component_index(j).first) + == Coupling::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + cell->get_fe().system_to_component_index(j).first) + == Coupling::nonzero)) + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + } + else + { + typename dealii::hp::DoFHandler::cell_iterator + neighbor = cell->neighbor(face); + // Refinement edges are taken care of + // by coarser cells + if (cell->neighbor_is_coarser(face)) + continue; + + typename dealii::hp::DoFHandler::face_iterator + cell_face = cell->face(face); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->n_children(); + ++sub_nr) + { + const typename dealii::hp::DoFHandler::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell); + sub_neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + for (unsigned int j=0; jget_fe().dofs_per_cell; + ++j) + { + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + sub_neighbor->get_fe().system_to_component_index(j).first) + == Coupling::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + sub_neighbor->get_fe().system_to_component_index(j).first) + == Coupling::nonzero)) + { + sparsity.add (dofs_on_this_cell[i], + dofs_on_other_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_other_cell[j]); + } + + if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == Coupling::always) + || + (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == Coupling::nonzero)) + { + sparsity.add (dofs_on_this_cell[j], + dofs_on_other_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_this_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_other_cell[i]); + } + } + } + sub_neighbor->face(neighbor_face)->set_user_flag (); + } + } + else + { + dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell); + neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + { + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == Coupling::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == Coupling::nonzero)) + { + sparsity.add (dofs_on_this_cell[i], + dofs_on_other_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_other_cell[j]); + } + + if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == Coupling::always) + || + (flux_mask(neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == Coupling::nonzero)) + { + sparsity.add (dofs_on_this_cell[j], + dofs_on_other_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_this_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_other_cell[i]); + } + } + } + neighbor->face(neighbor_face)->set_user_flag (); + } + } + } + } + } } } -} -template -void -DoFTools:: -make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const Table<2,Coupling> &int_mask, - const Table<2,Coupling> &flux_mask) -{ - // do the error checking and frame code - // here, and then pass on to more - // specialized functions in the internal - // namespace - const unsigned int n_dofs = dof.n_dofs(); - const unsigned int n_comp = dof.get_fe().n_components(); - - Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); - Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); - Assert (int_mask.n_rows() == n_comp, - ExcDimensionMismatch (int_mask.n_rows(), n_comp)); - Assert (int_mask.n_cols() == n_comp, - ExcDimensionMismatch (int_mask.n_cols(), n_comp)); - Assert (flux_mask.n_rows() == n_comp, - ExcDimensionMismatch (flux_mask.n_rows(), n_comp)); - Assert (flux_mask.n_cols() == n_comp, - ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); + + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const Table<2,Coupling> &int_mask, + const Table<2,Coupling> &flux_mask) + { + // do the error checking and frame code + // here, and then pass on to more + // specialized functions in the internal + // namespace + const unsigned int n_dofs = dof.n_dofs(); + const unsigned int n_comp = dof.get_fe().n_components(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (sparsity.n_cols() == n_dofs, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + Assert (int_mask.n_rows() == n_comp, + ExcDimensionMismatch (int_mask.n_rows(), n_comp)); + Assert (int_mask.n_cols() == n_comp, + ExcDimensionMismatch (int_mask.n_cols(), n_comp)); + Assert (flux_mask.n_rows() == n_comp, + ExcDimensionMismatch (flux_mask.n_rows(), n_comp)); + Assert (flux_mask.n_cols() == n_comp, + ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); - // Clear user flags because we will - // need them. But first we save - // them and make sure that we - // restore them later such that at - // the end of this function the - // Triangulation will be in the - // same state as it was at the - // beginning of this function. - std::vector user_flags; - dof.get_tria().save_user_flags(user_flags); - const_cast &>(dof.get_tria()).clear_user_flags (); - - internal::DoFTools::make_flux_sparsity_pattern (dof, sparsity, - int_mask, flux_mask); + // Clear user flags because we will + // need them. But first we save + // them and make sure that we + // restore them later such that at + // the end of this function the + // Triangulation will be in the + // same state as it was at the + // beginning of this function. + std::vector user_flags; + dof.get_tria().save_user_flags(user_flags); + const_cast &>(dof.get_tria()).clear_user_flags (); + + internal::make_flux_sparsity_pattern (dof, sparsity, + int_mask, flux_mask); - // finally restore the user flags - const_cast &>(dof.get_tria()).load_user_flags(user_flags); -} + // finally restore the user flags + const_cast &>(dof.get_tria()).load_user_flags(user_flags); + } @@ -1174,9 +1178,7 @@ make_flux_sparsity_pattern (const DH &dof, -namespace internal -{ - namespace DoFTools + namespace internal { namespace { @@ -1521,29 +1523,29 @@ namespace internal } - /** - * Make sure that the given @p - * face_interpolation_matrix - * pointer points to a valid - * matrix. If the pointer is zero - * beforehand, create an entry - * with the correct data. If it - * is nonzero, don't touch it. - */ + /** + * Make sure that the given @p + * face_interpolation_matrix + * pointer points to a valid + * matrix. If the pointer is zero + * beforehand, create an entry + * with the correct data. If it + * is nonzero, don't touch it. + */ template void ensure_existence_of_face_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - std_cxx1x::shared_ptr > &matrix) + const FiniteElement &fe2, + std_cxx1x::shared_ptr > &matrix) { - if (matrix == std_cxx1x::shared_ptr >()) - { - matrix = std_cxx1x::shared_ptr > - (new FullMatrix (fe2.dofs_per_face, - fe1.dofs_per_face)); - fe1.get_face_interpolation_matrix (fe2, - *matrix); - } + if (matrix == std_cxx1x::shared_ptr >()) + { + matrix = std_cxx1x::shared_ptr > + (new FullMatrix (fe2.dofs_per_face, + fe1.dofs_per_face)); + fe1.get_face_interpolation_matrix (fe2, + *matrix); + } } @@ -1555,19 +1557,19 @@ namespace internal template void ensure_existence_of_subface_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - const unsigned int subface, - std_cxx1x::shared_ptr > &matrix) + const FiniteElement &fe2, + const unsigned int subface, + std_cxx1x::shared_ptr > &matrix) { - if (matrix == std_cxx1x::shared_ptr >()) - { - matrix = std_cxx1x::shared_ptr > - (new FullMatrix (fe2.dofs_per_face, - fe1.dofs_per_face)); - fe1.get_subface_interpolation_matrix (fe2, - subface, - *matrix); - } + if (matrix == std_cxx1x::shared_ptr >()) + { + matrix = std_cxx1x::shared_ptr > + (new FullMatrix (fe2.dofs_per_face, + fe1.dofs_per_face)); + fe1.get_subface_interpolation_matrix (fe2, + subface, + *matrix); + } } @@ -1579,9 +1581,6 @@ namespace internal * master part as explained in * the @ref hp_paper "hp paper". */ -#ifdef DEAL_II_ANON_NAMESPACE_BUG - static -#endif void ensure_existence_of_split_face_matrix (const FullMatrix &face_interpolation_matrix, const std::vector &master_dof_mask, @@ -1592,10 +1591,10 @@ namespace internal static_cast(face_interpolation_matrix.n()), ExcInternalError()); - if (split_matrix == + if (split_matrix == std_cxx1x::shared_ptr,FullMatrix > >()) - { - split_matrix + { + split_matrix = std_cxx1x::shared_ptr,FullMatrix > > (new std::pair,FullMatrix >()); @@ -1638,46 +1637,46 @@ namespace internal } - // a template that can - // determine statically whether - // a given DoFHandler class - // supports different finite - // element elements + // a template that can + // determine statically whether + // a given DoFHandler class + // supports different finite + // element elements template struct DoFHandlerSupportsDifferentFEs { - static const bool value = true; + static const bool value = true; }; template struct DoFHandlerSupportsDifferentFEs< dealii::DoFHandler > { - static const bool value = false; + static const bool value = false; }; template struct DoFHandlerSupportsDifferentFEs< dealii::MGDoFHandler > { - static const bool value = false; + static const bool value = false; }; - /** - * A function that returns how - * many different finite - * elements a dof handler - * uses. This is one for non-hp - * DoFHandlers and - * dof_handler.get_fe().size() - * for the hp-versions. - */ + /** + * A function that returns how + * many different finite + * elements a dof handler + * uses. This is one for non-hp + * DoFHandlers and + * dof_handler.get_fe().size() + * for the hp-versions. + */ template unsigned int n_finite_elements (const dealii::hp::DoFHandler &dof_handler) { - return dof_handler.get_fe().size(); + return dof_handler.get_fe().size(); } @@ -1685,19 +1684,19 @@ namespace internal unsigned int n_finite_elements (const DH &) { - return 1; + return 1; } - /** - * For a given face belonging - * to an active cell that - * borders to a more refined - * cell, return the fe_index of - * the most dominating finite - * element used on any of the - * face's subfaces. - */ + /** + * For a given face belonging + * to an active cell that + * borders to a more refined + * cell, return the fe_index of + * the most dominating finite + * element used on any of the + * face's subfaces. + */ template unsigned int get_most_dominating_subface_fe_index (const face_iterator &face) @@ -1707,85 +1706,82 @@ namespace internal const unsigned int spacedim = face_iterator::AccessorType::space_dimension; - unsigned int dominating_subface_no = 0; - for (; dominating_subface_non_children(); - ++dominating_subface_no) - { - // each of the subfaces - // can have only a single - // fe_index associated - // with them, since there - // is no cell on the - // other side - Assert (face->child(dominating_subface_no) - ->n_active_fe_indices() - == 1, - ExcInternalError()); - - const FiniteElement & - this_subface_fe = (face->child(dominating_subface_no) - ->get_fe (face->child(dominating_subface_no) - ->nth_active_fe_index(0))); - - FiniteElementDomination::Domination - domination = FiniteElementDomination::either_element_can_dominate; - for (unsigned int sf=0; sfn_children(); ++sf) - if (sf != dominating_subface_no) - { - const FiniteElement & - that_subface_fe = (face->child(sf) - ->get_fe (face->child(sf) - ->nth_active_fe_index(0))); - - domination = domination & - this_subface_fe.compare_for_face_domination(that_subface_fe); - } - - // see if the element - // on this subface is - // able to dominate - // the ones on all - // other subfaces, - // and if so take it - if ((domination == FiniteElementDomination::this_element_dominates) - || - (domination == FiniteElementDomination::either_element_can_dominate)) - break; - } - - // check that we have - // found one such subface - Assert (dominating_subface_no < face->n_children(), - ExcNotImplemented()); - - // return the finite element - // index used on it. note - // that only a single fe can - // be active on such subfaces - return face->child (dominating_subface_no)->nth_active_fe_index(0); + unsigned int dominating_subface_no = 0; + for (; dominating_subface_non_children(); + ++dominating_subface_no) + { + // each of the subfaces + // can have only a single + // fe_index associated + // with them, since there + // is no cell on the + // other side + Assert (face->child(dominating_subface_no) + ->n_active_fe_indices() + == 1, + ExcInternalError()); + + const FiniteElement & + this_subface_fe = (face->child(dominating_subface_no) + ->get_fe (face->child(dominating_subface_no) + ->nth_active_fe_index(0))); + + FiniteElementDomination::Domination + domination = FiniteElementDomination::either_element_can_dominate; + for (unsigned int sf=0; sfn_children(); ++sf) + if (sf != dominating_subface_no) + { + const FiniteElement & + that_subface_fe = (face->child(sf) + ->get_fe (face->child(sf) + ->nth_active_fe_index(0))); + + domination = domination & + this_subface_fe.compare_for_face_domination(that_subface_fe); + } + + // see if the element + // on this subface is + // able to dominate + // the ones on all + // other subfaces, + // and if so take it + if ((domination == FiniteElementDomination::this_element_dominates) + || + (domination == FiniteElementDomination::either_element_can_dominate)) + break; + } + + // check that we have + // found one such subface + Assert (dominating_subface_no < face->n_children(), + ExcNotImplemented()); + + // return the finite element + // index used on it. note + // that only a single fe can + // be active on such subfaces + return face->child (dominating_subface_no)->nth_active_fe_index(0); } - /** + /** * Copy constraints into a constraint * matrix object. * - * This function removes zero - * constraints and those, which - * constrain a DoF which was - * already eliminated in one of - * the previous steps of the hp - * hanging node procedure. + * This function removes zero + * constraints and those, which + * constrain a DoF which was + * already eliminated in one of + * the previous steps of the hp + * hanging node procedure. * * It also suppresses very small * entries in the constraint matrix to * avoid making the sparsity pattern * fuller than necessary. - */ -#ifdef DEAL_II_ANON_NAMESPACE_BUG - static -#endif + */ void filter_constraints (const std::vector &master_dofs, const std::vector &slave_dofs, @@ -1815,25 +1811,25 @@ namespace internal for (unsigned int row=0; row!=n_slave_dofs; ++row) - if (constraints.is_constrained (slave_dofs[row]) == false) - { - bool constraint_already_satisfied = false; - - // Check if we have an identity - // constraint, which is already - // satisfied by unification of - // the corresponding global dof - // indices - for (unsigned int i=0; i= 1e-14*abs_sum)) constraints.add_entry (slave_dofs[row], master_dofs[i], face_constraints (row,i)); constraints.set_inhomogeneity (slave_dofs[row], 0.); - } - } + } + } } } - static void make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &, ConstraintMatrix &) @@ -1886,7 +1881,6 @@ namespace internal - static void make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1> &, ConstraintMatrix &, @@ -1897,7 +1891,6 @@ namespace internal } - static void make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &, ConstraintMatrix &) @@ -1908,7 +1901,6 @@ namespace internal - static void make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1> &, ConstraintMatrix &, @@ -1919,7 +1911,6 @@ namespace internal } - static void make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, ConstraintMatrix &/*constraints*/) @@ -1933,7 +1924,6 @@ namespace internal - static void make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, ConstraintMatrix &/*constraints*/, @@ -1947,7 +1937,6 @@ namespace internal } - static void make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &, ConstraintMatrix &) @@ -1958,7 +1947,6 @@ namespace internal - static void make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, ConstraintMatrix &, @@ -1971,7 +1959,6 @@ namespace internal // currently not used but may be in the future: -// static // void // make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1,2> &, // ConstraintMatrix &) @@ -1982,7 +1969,6 @@ namespace internal -// static // void // make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1,2> &, // ConstraintMatrix &, @@ -1993,7 +1979,6 @@ namespace internal // } -// static // void // make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1,2> &/*dof_handler*/, // ConstraintMatrix &/*constraints*/, @@ -2010,7 +1995,6 @@ namespace internal template - static void make_oldstyle_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints, @@ -2167,7 +2151,6 @@ namespace internal template - static void make_oldstyle_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints, @@ -2275,10 +2258,10 @@ namespace internal const FiniteElement &fe = cell->get_fe(); const unsigned int fe_index = cell->active_fe_index(); - const unsigned int n_dofs_on_mother = fe.dofs_per_face; - const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+ - 12*fe.dofs_per_line+ - 4*fe.dofs_per_quad); + const unsigned int n_dofs_on_mother = fe.dofs_per_face; + const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+ + 12*fe.dofs_per_line+ + 4*fe.dofs_per_quad); //TODO[TL]: think about this and the following in case of anisotropic refinement dofs_on_mother.resize (n_dofs_on_mother); @@ -2406,7 +2389,6 @@ namespace internal template - static void make_hp_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints) @@ -2424,11 +2406,11 @@ namespace internal const unsigned int spacedim = DH::space_dimension; - // a matrix to be used for - // constraints below. declared - // here and simply resized down - // below to avoid permanent - // re-allocation of memory + // a matrix to be used for + // constraints below. declared + // here and simply resized down + // below to avoid permanent + // re-allocation of memory FullMatrix constraint_matrix; // similarly have arrays that @@ -2440,21 +2422,21 @@ namespace internal std::vector slave_dofs; std::vector scratch_dofs; - // caches for the face and - // subface interpolation - // matrices between different - // (or the same) finite - // elements. we compute them - // only once, namely the first - // time they are needed, and - // then just reuse them + // caches for the face and + // subface interpolation + // matrices between different + // (or the same) finite + // elements. we compute them + // only once, namely the first + // time they are needed, and + // then just reuse them Table<2,std_cxx1x::shared_ptr > > - face_interpolation_matrices (n_finite_elements (dof_handler), - n_finite_elements (dof_handler)); + face_interpolation_matrices (n_finite_elements (dof_handler), + n_finite_elements (dof_handler)); Table<3,std_cxx1x::shared_ptr > > - subface_interpolation_matrices (n_finite_elements (dof_handler), - n_finite_elements (dof_handler), - GeometryInfo::max_children_per_face); + subface_interpolation_matrices (n_finite_elements (dof_handler), + n_finite_elements (dof_handler), + GeometryInfo::max_children_per_face); // similarly have a cache for // the matrices that are split @@ -2466,7 +2448,7 @@ namespace internal // interpolation matrix as // described in the @ref hp_paper "hp paper" Table<2,std_cxx1x::shared_ptr,FullMatrix > > > - split_face_interpolation_matrices (n_finite_elements (dof_handler), + split_face_interpolation_matrices (n_finite_elements (dof_handler), n_finite_elements (dof_handler)); // finally, for each pair of finite @@ -3226,1530 +3208,1527 @@ namespace internal } } } -} - -template -void -DoFTools::make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints) -{ - // Decide whether to use the - // new or old make_hanging_node_constraints - // function. If all the FiniteElement - // or all elements in a FECollection support - // the new face constraint matrix, the - // new code will be used. - // Otherwise, the old implementation is used - // for the moment. - if (dof_handler.get_fe().hp_constraints_are_implemented ()) - internal::DoFTools:: - make_hp_hanging_node_constraints (dof_handler, - constraints); - else - internal::DoFTools:: - make_oldstyle_hanging_node_constraints (dof_handler, - constraints, - dealii::internal::int2type()); -} - - -namespace internal -{ - // this internal function assigns to each dof - // the respective component of the vector - // system. if use_blocks is set, then the - // assignment is done by blocks, not by - // components, as specified by - // component_select. The additional argument - // component_select only is used for - // non-primitive FEs, where we need it since - // more components couple, and no unique - // component can be assigned. Then, we sort - // them to the first selected component of the - // vector system. template - inline void - extract_dofs_by_component (const DH &dof, - const std::vector &component_select, - const bool sort_by_blocks, - std::vector &dofs_by_component) + make_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints) { - const dealii::hp::FECollection - fe_collection (dof.get_fe()); - Assert (fe_collection.n_components() < 256, ExcNotImplemented()); - Assert (dofs_by_component.size() == dof.n_locally_owned_dofs(), - ExcDimensionMismatch(dofs_by_component.size(), - dof.n_locally_owned_dofs())); - - // next set up a table for the - // degrees of freedom on each of - // the cells whether it is - // something interesting or not - std::vector > local_component_association - (fe_collection.size()); - for (unsigned int f=0; f &fe = - fe_collection[f]; - local_component_association[f].resize(fe.dofs_per_cell); - if (sort_by_blocks == true) - { - for (unsigned int i=0; i indices; - for (typename DH::active_cell_iterator c=dof.begin_active(); - c!=dof.end(); ++ c) - if (!c->is_artificial() && !c->is_ghost()) - { - const unsigned int fe_index = c->active_fe_index(); - const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell; - indices.resize(dofs_per_cell); - c->get_dof_indices(indices); - for (unsigned int i=0; i()); } -} -template -void DoFTools::distribute_cell_to_dof_vector ( - const DH &dof_handler, - const Vector &cell_data, - Vector &dof_data, - const unsigned int component) -{ - const Triangulation &tria = dof_handler.get_tria(); - - Assert (cell_data.size()==tria.n_active_cells(), - ExcWrongSize (cell_data.size(), tria.n_active_cells())); - Assert (dof_data.size()==dof_handler.n_dofs(), - ExcWrongSize (dof_data.size(), dof_handler.n_dofs())); - Assert (component < n_components(dof_handler), - ExcInvalidComponent(component, n_components(dof_handler))); - Assert (fe_is_primitive(dof_handler) == true, - ExcFENotPrimitive()); - - // store a flag whether we should care - // about different components. this is - // just a simplification, we could ask - // for this at every single place - // equally well - const bool consider_components = (n_components(dof_handler) != 1); - - // zero out the components that we - // will touch - if (consider_components == false) - dof_data = 0; - else + namespace internal + { + // this internal function assigns to each dof + // the respective component of the vector + // system. if use_blocks is set, then the + // assignment is done by blocks, not by + // components, as specified by + // component_select. The additional argument + // component_select only is used for + // non-primitive FEs, where we need it since + // more components couple, and no unique + // component can be assigned. Then, we sort + // them to the first selected component of the + // vector system. + template + inline + void + extract_dofs_by_component (const DH &dof, + const std::vector &component_select, + const bool sort_by_blocks, + std::vector &dofs_by_component) { - std::vector component_dofs (dof_handler.n_locally_owned_dofs()); - std::vector component_mask (dof_handler.get_fe().n_components(), - false); - component_mask[component] = true; - internal::extract_dofs_by_component (dof_handler, component_mask, - false, component_dofs); - - for (unsigned int i=0; i(component)) - dof_data(i) = 0; - } - - // count how often we have added a value - // in the sum for each dof - std::vector touch_count (dof_handler.n_dofs(), 0); - - typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - std::vector dof_indices; - dof_indices.reserve (max_dofs_per_cell(dof_handler)); + const dealii::hp::FECollection + fe_collection (dof.get_fe()); + Assert (fe_collection.n_components() < 256, ExcNotImplemented()); + Assert (dofs_by_component.size() == dof.n_locally_owned_dofs(), + ExcDimensionMismatch(dofs_by_component.size(), + dof.n_locally_owned_dofs())); + + // next set up a table for the + // degrees of freedom on each of + // the cells whether it is + // something interesting or not + std::vector > local_component_association + (fe_collection.size()); + for (unsigned int f=0; f &fe = + fe_collection[f]; + local_component_association[f].resize(fe.dofs_per_cell); + if (sort_by_blocks == true) + { + for (unsigned int i=0; iget_fe().dofs_per_cell; - dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (dof_indices); - - for (unsigned int i=0; iget_fe().system_to_component_index(i).first == component)) + // then loop over all cells and do + // the work + std::vector indices; + for (typename DH::active_cell_iterator c=dof.begin_active(); + c!=dof.end(); ++ c) + if (!c->is_artificial() && !c->is_ghost()) { - // sum up contribution of the - // present_cell to this dof - dof_data(dof_indices[i]) += cell_data(present_cell); - // note that we added another - // summand - ++touch_count[dof_indices[i]]; + const unsigned int fe_index = c->active_fe_index(); + const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell; + indices.resize(dofs_per_cell); + c->get_dof_indices(indices); + for (unsigned int i=0; i + void distribute_cell_to_dof_vector ( + const DH &dof_handler, + const Vector &cell_data, + Vector &dof_data, + const unsigned int component) + { + const Triangulation &tria = dof_handler.get_tria(); + + Assert (cell_data.size()==tria.n_active_cells(), + ExcWrongSize (cell_data.size(), tria.n_active_cells())); + Assert (dof_data.size()==dof_handler.n_dofs(), + ExcWrongSize (dof_data.size(), dof_handler.n_dofs())); + Assert (component < n_components(dof_handler), + ExcInvalidComponent(component, n_components(dof_handler))); + Assert (fe_is_primitive(dof_handler) == true, + ExcFENotPrimitive()); + + // store a flag whether we should care + // about different components. this is + // just a simplification, we could ask + // for this at every single place + // equally well + const bool consider_components = (n_components(dof_handler) != 1); + + // zero out the components that we + // will touch + if (consider_components == false) + dof_data = 0; + else + { + std::vector component_dofs (dof_handler.n_locally_owned_dofs()); + std::vector component_mask (dof_handler.get_fe().n_components(), + false); + component_mask[component] = true; + internal::extract_dofs_by_component (dof_handler, component_mask, + false, component_dofs); + + for (unsigned int i=0; i(component)) + dof_data(i) = 0; + } -template -void -DoFTools::extract_dofs ( - const DoFHandler &dof, - const std::vector &component_select, - std::vector &selected_dofs, - const bool count_by_blocks) -{ - const FiniteElement &fe = dof.get_fe(); + // count how often we have added a value + // in the sum for each dof + std::vector touch_count (dof_handler.n_dofs(), 0); - if (count_by_blocks == true) - { - Assert(component_select.size() == fe.n_blocks(), - ExcDimensionMismatch(component_select.size(), fe.n_blocks())); - } - else - { - Assert(component_select.size() == n_components(dof), - ExcDimensionMismatch(component_select.size(), n_components(dof))); - } + typename DH::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + std::vector dof_indices; + dof_indices.reserve (max_dofs_per_cell(dof_handler)); - Assert(selected_dofs.size() == dof.n_locally_owned_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs())); + for (unsigned int present_cell = 0; cell!=endc; ++cell, ++present_cell) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (dof_indices); - // two special cases: no component - // is selected, and all components - // are selected; both rather - // stupid, but easy to catch - if (std::count (component_select.begin(), component_select.end(), true) - == 0) - { - std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false); - return; - } - else if (std::count (component_select.begin(), component_select.end(), true) - == static_cast(component_select.size())) - { - std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true); - return; - } + for (unsigned int i=0; iget_fe().system_to_component_index(i).first == component)) + { + // sum up contribution of the + // present_cell to this dof + dof_data(dof_indices[i]) += cell_data(present_cell); + // note that we added another + // summand + ++touch_count[dof_indices[i]]; + } + } + // compute the mean value on all the + // dofs by dividing with the number + // of summands. + for (unsigned int i=0; i dofs_by_component (dof.n_locally_owned_dofs()); - internal::extract_dofs_by_component (dof, component_select, count_by_blocks, - dofs_by_component); - for (unsigned int i=0; i + void + extract_dofs ( + const DoFHandler &dof, + const std::vector &component_select, + std::vector &selected_dofs, + const bool count_by_blocks) + { + const FiniteElement &fe = dof.get_fe(); + if (count_by_blocks == true) + { + Assert(component_select.size() == fe.n_blocks(), + ExcDimensionMismatch(component_select.size(), fe.n_blocks())); + } + else + { + Assert(component_select.size() == n_components(dof), + ExcDimensionMismatch(component_select.size(), n_components(dof))); + } + Assert(selected_dofs.size() == dof.n_locally_owned_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs())); -template -void -DoFTools::extract_dofs ( - const hp::DoFHandler &dof, - const std::vector &component_select, - std::vector &selected_dofs, - const bool count_by_blocks) -{ - const FiniteElement &fe = dof.begin_active()->get_fe(); + // two special cases: no component + // is selected, and all components + // are selected; both rather + // stupid, but easy to catch + if (std::count (component_select.begin(), component_select.end(), true) + == 0) + { + std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false); + return; + } + else if (std::count (component_select.begin(), component_select.end(), true) + == static_cast(component_select.size())) + { + std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true); + return; + } - if (count_by_blocks == true) - { - Assert(component_select.size() == fe.n_blocks(), - ExcDimensionMismatch(component_select.size(), fe.n_blocks())); - } - else - { - Assert(component_select.size() == n_components(dof), - ExcDimensionMismatch(component_select.size(), n_components(dof))); - } - Assert(selected_dofs.size() == dof.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs())); + // preset all values by false + std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false); - // two special cases: no component - // is selected, and all components - // are selected; both rather - // stupid, but easy to catch - if (std::count (component_select.begin(), component_select.end(), true) - == 0) - { - std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); - return; - }; - if (std::count (component_select.begin(), component_select.end(), true) - == static_cast(component_select.size())) - { - std::fill_n (selected_dofs.begin(), dof.n_dofs(), true); - return; - }; - - - // preset all values by false - std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); - - // if we count by blocks, we need to extract - // the association of blocks with local dofs, - // and then go through all the cells and set - // the properties according to this - // info. Otherwise, we let the function - // extract_dofs_by_component function do the - // job. - std::vector dofs_by_component (dof.n_dofs()); - internal::extract_dofs_by_component (dof, component_select, count_by_blocks, - dofs_by_component); - - for (unsigned int i=0; i dofs_by_component (dof.n_locally_owned_dofs()); + internal::extract_dofs_by_component (dof, component_select, count_by_blocks, + dofs_by_component); + for (unsigned int i=0; i -void -DoFTools::extract_level_dofs( - const unsigned int level, - const MGDoFHandler &dof, - const std::vector &component_select, - std::vector &selected_dofs, - const bool count_by_blocks) -{ - const FiniteElement& fe = dof.get_fe(); - if (count_by_blocks == true) - { - Assert(component_select.size() == fe.n_blocks(), - ExcDimensionMismatch(component_select.size(), fe.n_blocks())); - } - else - { - Assert(component_select.size() == fe.n_components(), - ExcDimensionMismatch(component_select.size(), fe.n_components())); - } + template + void + extract_dofs ( + const hp::DoFHandler &dof, + const std::vector &component_select, + std::vector &selected_dofs, + const bool count_by_blocks) + { + const FiniteElement &fe = dof.begin_active()->get_fe(); - Assert(selected_dofs.size() == dof.n_dofs(level), - ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); + if (count_by_blocks == true) + { + Assert(component_select.size() == fe.n_blocks(), + ExcDimensionMismatch(component_select.size(), fe.n_blocks())); + } + else + { + Assert(component_select.size() == n_components(dof), + ExcDimensionMismatch(component_select.size(), n_components(dof))); + } - // two special cases: no component - // is selected, and all components - // are selected, both rather - // stupid, but easy to catch - if (std::count (component_select.begin(), component_select.end(), true) - == 0) - { - std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false); - return; - }; - if (std::count (component_select.begin(), component_select.end(), true) - == static_cast(component_select.size())) - { - std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true); - return; - }; + Assert(selected_dofs.size() == dof.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs())); - // preset all values by false - std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false); + // two special cases: no component + // is selected, and all components + // are selected; both rather + // stupid, but easy to catch + if (std::count (component_select.begin(), component_select.end(), true) + == 0) + { + std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); + return; + }; + if (std::count (component_select.begin(), component_select.end(), true) + == static_cast(component_select.size())) + { + std::fill_n (selected_dofs.begin(), dof.n_dofs(), true); + return; + }; - // next set up a table for the - // degrees of freedom on each of - // the cells whether it is - // something interesting or not - std::vector local_selected_dofs (fe.dofs_per_cell, false); - for (unsigned int i=0; i indices(fe.dofs_per_cell); - typename MGDoFHandler::cell_iterator c; - for (c = dof.begin(level) ; c != dof.end(level) ; ++ c) - { - c->get_mg_dof_indices(indices); - for (unsigned int i=0; i dofs_by_component (dof.n_dofs()); + internal::extract_dofs_by_component (dof, component_select, count_by_blocks, + dofs_by_component); + for (unsigned int i=0; i -void -DoFTools::extract_boundary_dofs (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators) -{ - Assert (component_select.size() == n_components(dof_handler), - ExcWrongSize (component_select.size(), - n_components(dof_handler))); - Assert (boundary_indicators.find (255) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); - const unsigned int dim=DH::dimension; - - // let's see whether we have to - // check for certain boundary - // indicators or whether we can - // accept all - const bool check_boundary_indicator = (boundary_indicators.size() != 0); - - // also see whether we have to - // check whether a certain vector - // component is selected, or all - const bool check_vector_component - = (component_select != std::vector(component_select.size(), - true)); - - // clear and reset array by default - // values - selected_dofs.clear (); - selected_dofs.resize (dof_handler.n_dofs(), false); - std::vector face_dof_indices; - face_dof_indices.reserve (max_dofs_per_face(dof_handler)); - - // now loop over all cells and - // check whether their faces are at - // the boundary. note that we need - // not take special care of single - // lines being at the boundary - // (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every isolated boundary - // line is also part of a boundary - // face which we will be visiting - // sooner or later - for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); - cell!=dof_handler.end(); ++cell) - for (unsigned int face=0; - face::faces_per_cell; ++face) - if (cell->at_boundary(face)) - if (! check_boundary_indicator || - (boundary_indicators.find (cell->face(face)->boundary_indicator()) - != boundary_indicators.end())) - { - const FiniteElement &fe = cell->get_fe(); - - const unsigned int dofs_per_face = fe.dofs_per_face; - face_dof_indices.resize (dofs_per_face); - cell->face(face)->get_dof_indices (face_dof_indices, - cell->active_fe_index()); - - for (unsigned int i=0; i + void + extract_level_dofs( + const unsigned int level, + const MGDoFHandler &dof, + const std::vector &component_select, + std::vector &selected_dofs, + const bool count_by_blocks) + { + const FiniteElement& fe = dof.get_fe(); -template -void -DoFTools::extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators) -{ - Assert (component_select.size() == n_components(dof_handler), - ExcWrongSize (component_select.size(), - n_components(dof_handler))); - Assert (boundary_indicators.find (255) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); - - // let's see whether we have to - // check for certain boundary - // indicators or whether we can - // accept all - const bool check_boundary_indicator = (boundary_indicators.size() != 0); - - // also see whether we have to - // check whether a certain vector - // component is selected, or all - const bool check_vector_component - = (component_select != std::vector(component_select.size(), - true)); - - // clear and reset array by default - // values - selected_dofs.clear (); - selected_dofs.resize (dof_handler.n_dofs(), false); - std::vector cell_dof_indices; - cell_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - - // now loop over all cells and - // check whether their faces are at - // the boundary. note that we need - // not take special care of single - // lines being at the boundary - // (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every isolated boundary - // line is also part of a boundary - // face which we will be visiting - // sooner or later - for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); - cell!=dof_handler.end(); ++cell) - for (unsigned int face=0; - face::faces_per_cell; ++face) - if (cell->at_boundary(face)) - if (! check_boundary_indicator || - (boundary_indicators.find (cell->face(face)->boundary_indicator()) - != boundary_indicators.end())) - { - const FiniteElement &fe = cell->get_fe(); + if (count_by_blocks == true) + { + Assert(component_select.size() == fe.n_blocks(), + ExcDimensionMismatch(component_select.size(), fe.n_blocks())); + } + else + { + Assert(component_select.size() == fe.n_components(), + ExcDimensionMismatch(component_select.size(), fe.n_components())); + } - const unsigned int dofs_per_cell = fe.dofs_per_cell; - cell_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (cell_dof_indices); + Assert(selected_dofs.size() == dof.n_dofs(level), + ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); - for (unsigned int i=0; i(component_select.size())) + { + std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true); + return; + }; - selected_dofs[cell_dof_indices[i]] - = (component_select[first_nonzero_comp] - == true); - } - } + // preset all values by false + std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false); + + // next set up a table for the + // degrees of freedom on each of + // the cells whether it is + // something interesting or not + std::vector local_selected_dofs (fe.dofs_per_cell, false); + for (unsigned int i=0; i indices(fe.dofs_per_cell); + typename MGDoFHandler::cell_iterator c; + for (c = dof.begin(level) ; c != dof.end(level) ; ++ c) + { + c->get_mg_dof_indices(indices); + for (unsigned int i=0; i + void + extract_boundary_dofs (const DH &dof_handler, + const std::vector &component_select, + std::vector &selected_dofs, + const std::set &boundary_indicators) { - template - void extract_hanging_node_dofs (const dealii::DoFHandler<1,spacedim> &dof_handler, - std::vector &selected_dofs) - { - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); - // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - - // there are no hanging nodes in 1d - } + Assert (component_select.size() == n_components(dof_handler), + ExcWrongSize (component_select.size(), + n_components(dof_handler))); + Assert (boundary_indicators.find (255) == boundary_indicators.end(), + ExcInvalidBoundaryIndicator()); + const unsigned int dim=DH::dimension; + + // let's see whether we have to + // check for certain boundary + // indicators or whether we can + // accept all + const bool check_boundary_indicator = (boundary_indicators.size() != 0); + + // also see whether we have to + // check whether a certain vector + // component is selected, or all + const bool check_vector_component + = (component_select != std::vector(component_select.size(), + true)); + + // clear and reset array by default + // values + selected_dofs.clear (); + selected_dofs.resize (dof_handler.n_dofs(), false); + std::vector face_dof_indices; + face_dof_indices.reserve (max_dofs_per_face(dof_handler)); + + // now loop over all cells and + // check whether their faces are at + // the boundary. note that we need + // not take special care of single + // lines being at the boundary + // (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every isolated boundary + // line is also part of a boundary + // face which we will be visiting + // sooner or later + for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + for (unsigned int face=0; + face::faces_per_cell; ++face) + if (cell->at_boundary(face)) + if (! check_boundary_indicator || + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != boundary_indicators.end())) + { + const FiniteElement &fe = cell->get_fe(); + const unsigned int dofs_per_face = fe.dofs_per_face; + face_dof_indices.resize (dofs_per_face); + cell->face(face)->get_dof_indices (face_dof_indices, + cell->active_fe_index()); - template - void extract_hanging_node_dofs (const dealii::DoFHandler<2,spacedim> &dof_handler, - std::vector &selected_dofs) - { - const unsigned int dim = 2; + for (unsigned int i=0; i &fe = dof_handler.get_fe(); - // this function is similar to the - // make_sparsity_pattern function, - // see there for more information - typename dealii::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->has_children()) + template + void + extract_dofs_with_support_on_boundary (const DH &dof_handler, + const std::vector &component_select, + std::vector &selected_dofs, + const std::set &boundary_indicators) + { + Assert (component_select.size() == n_components(dof_handler), + ExcWrongSize (component_select.size(), + n_components(dof_handler))); + Assert (boundary_indicators.find (255) == boundary_indicators.end(), + ExcInvalidBoundaryIndicator()); + + // let's see whether we have to + // check for certain boundary + // indicators or whether we can + // accept all + const bool check_boundary_indicator = (boundary_indicators.size() != 0); + + // also see whether we have to + // check whether a certain vector + // component is selected, or all + const bool check_vector_component + = (component_select != std::vector(component_select.size(), + true)); + + // clear and reset array by default + // values + selected_dofs.clear (); + selected_dofs.resize (dof_handler.n_dofs(), false); + std::vector cell_dof_indices; + cell_dof_indices.reserve (max_dofs_per_cell(dof_handler)); + + // now loop over all cells and + // check whether their faces are at + // the boundary. note that we need + // not take special care of single + // lines being at the boundary + // (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every isolated boundary + // line is also part of a boundary + // face which we will be visiting + // sooner or later + for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + for (unsigned int face=0; + face::faces_per_cell; ++face) + if (cell->at_boundary(face)) + if (! check_boundary_indicator || + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != boundary_indicators.end())) { - const typename dealii::DoFHandler::line_iterator - line = cell->face(face); + const FiniteElement &fe = cell->get_fe(); - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + cell_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (cell_dof_indices); - for (unsigned int child=0; child<2; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - selected_dofs[line->child(child)->dof_index(dof)] = true; + for (unsigned int i=0; i - void extract_hanging_node_dofs (const dealii::DoFHandler<3,spacedim> &dof_handler, - std::vector &selected_dofs) - { - const unsigned int dim = 3; - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); - // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - const FiniteElement &fe = dof_handler.get_fe(); + namespace internal + { + namespace + { + template + void extract_hanging_node_dofs (const dealii::DoFHandler<1,spacedim> &dof_handler, + std::vector &selected_dofs) + { + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + // preset all values by false + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - // this function is similar to the - // make_sparsity_pattern function, - // see there for more information + // there are no hanging nodes in 1d + } - typename dealii::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->has_children()) - { - const typename dealii::DoFHandler::face_iterator - face = cell->face(f); - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true; + template + void extract_hanging_node_dofs (const dealii::DoFHandler<2,spacedim> &dof_handler, + std::vector &selected_dofs) + { + const unsigned int dim = 2; + + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + // preset all values by false + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + + const FiniteElement &fe = dof_handler.get_fe(); + + // this function is similar to the + // make_sparsity_pattern function, + // see there for more information + typename dealii::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + { + const typename dealii::DoFHandler::line_iterator + line = cell->face(face); - // dof numbers on the centers of - // the lines bounding this face - for (unsigned int line=0; line<4; ++line) for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true; - - // next the dofs on the lines interior - // to the face; the order of these - // lines is laid down in the - // FiniteElement class documentation - for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(1)->line(2)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(2)->line(3)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(3)->line(0)->dof_index(dof)] = true; - - // dofs on the bordering lines - for (unsigned int line=0; line<4; ++line) + selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true; + for (unsigned int child=0; child<2; ++child) for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true; - - // finally, for the dofs interior - // to the four child faces - for (unsigned int child=0; child<4; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) - selected_dofs[face->child(child)->dof_index(dof)] = true; - } - } - } -} - + selected_dofs[line->child(child)->dof_index(dof)] = true; + } + } -template -void -DoFTools:: -extract_hanging_node_dofs (const DoFHandler &dof_handler, - std::vector &selected_dofs) -{ - internal::DoFTools::extract_hanging_node_dofs (dof_handler, - selected_dofs); -} + template + void extract_hanging_node_dofs (const dealii::DoFHandler<3,spacedim> &dof_handler, + std::vector &selected_dofs) + { + const unsigned int dim = 3; + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + // preset all values by false + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + const FiniteElement &fe = dof_handler.get_fe(); -template -void -DoFTools::extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id_t subdomain_id, - std::vector &selected_dofs) -{ - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); - - // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - - std::vector local_dof_indices; - local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - - // this function is similar to the - // make_sparsity_pattern function, - // see there for more information - typename DH::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == subdomain_id) - { - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->has_children()) + { + const typename dealii::DoFHandler::face_iterator + face = cell->face(f); + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true; -template -void -DoFTools::extract_locally_owned_dofs (const DH & dof_handler, - IndexSet & dof_set) -{ - // collect all the locally owned dofs - dof_set = dof_handler.locally_owned_dofs(); - dof_set.compress (); -} + // dof numbers on the centers of + // the lines bounding this face + for (unsigned int line=0; line<4; ++line) + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true; + // next the dofs on the lines interior + // to the face; the order of these + // lines is laid down in the + // FiniteElement class documentation + for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(1)->line(2)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(2)->line(3)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(3)->line(0)->dof_index(dof)] = true; + // dofs on the bordering lines + for (unsigned int line=0; line<4; ++line) + for (unsigned int child=0; child<2; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true; -template -void -DoFTools::extract_locally_active_dofs (const DH & dof_handler, - IndexSet & dof_set) -{ - // collect all the locally owned dofs - dof_set = dof_handler.locally_owned_dofs(); - - // add the DoF on the adjacent ghost cells - // to the IndexSet, cache them in a - // set. need to check each dof manually - // because we can't be sure that the dof - // range of locally_owned_dofs is really - // contiguous. - std::vector dof_indices; - std::set global_dof_indices; - - typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - if (!cell->is_ghost() && !cell->is_artificial()) - { - dof_indices.resize(cell->get_fe().dofs_per_cell); - cell->get_dof_indices(dof_indices); - - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); + // finally, for the dofs interior + // to the four child faces + for (unsigned int child=0; child<4; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) + selected_dofs[face->child(child)->dof_index(dof)] = true; + } } + } + } - dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); - dof_set.compress(); -} + + template + void + + extract_hanging_node_dofs (const DoFHandler &dof_handler, + std::vector &selected_dofs) + { + internal::extract_hanging_node_dofs (dof_handler, + selected_dofs); + } -template -void -DoFTools::extract_locally_relevant_dofs (const DH & dof_handler, - IndexSet & dof_set) -{ - // collect all the locally owned dofs - dof_set = dof_handler.locally_owned_dofs(); - - // add the DoF on the adjacent ghost cells - // to the IndexSet, cache them in a - // set. need to check each dof manually - // because we can't be sure that the dof - // range of locally_owned_dofs is really - // contiguous. - std::vector dof_indices; - std::set global_dof_indices; - - typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->is_ghost()) - { - dof_indices.resize(cell->get_fe().dofs_per_cell); - cell->get_dof_indices(dof_indices); - - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); - } + template + void + extract_subdomain_dofs (const DH &dof_handler, + const types::subdomain_id_t subdomain_id, + std::vector &selected_dofs) + { + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + + // preset all values by false + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + + std::vector local_dof_indices; + local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); + + // this function is similar to the + // make_sparsity_pattern function, + // see there for more information + typename DH::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == subdomain_id) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i + void + extract_locally_owned_dofs (const DH & dof_handler, + IndexSet & dof_set) + { + // collect all the locally owned dofs + dof_set = dof_handler.locally_owned_dofs(); + dof_set.compress (); + } -template -void -DoFTools::extract_constant_modes (const DH &dof_handler, - const std::vector &component_select, - std::vector > &constant_modes) -{ - const unsigned int n_components = dof_handler.get_fe().n_components(); - Assert (n_components == component_select.size(), - ExcDimensionMismatch(n_components, - component_select.size())); - std::vector localized_component (n_components, - numbers::invalid_unsigned_int); - unsigned int n_components_selected = 0; - for (unsigned int i=0; i dofs_by_component (dof_handler.n_locally_owned_dofs()); - internal::extract_dofs_by_component (dof_handler, component_select, false, - dofs_by_component); - unsigned int n_selected_dofs = 0; - for (unsigned int i=0; i(n_selected_dofs, - false)); - std::vector component_list (n_components, 0); - for (unsigned int d=0; d + void + extract_locally_active_dofs (const DH & dof_handler, + IndexSet & dof_set) + { + // collect all the locally owned dofs + dof_set = dof_handler.locally_owned_dofs(); + + // add the DoF on the adjacent ghost cells + // to the IndexSet, cache them in a + // set. need to check each dof manually + // because we can't be sure that the dof + // range of locally_owned_dofs is really + // contiguous. + std::vector dof_indices; + std::set global_dof_indices; + + typename DH::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (!cell->is_ghost() && !cell->is_artificial()) + { + dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(dof_indices); + + for (std::vector::iterator it=dof_indices.begin(); + it!=dof_indices.end(); + ++it) + if (!dof_set.is_element(*it)) + global_dof_indices.insert(*it); + } + dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); -template -void -DoFTools::get_active_fe_indices (const DH &dof_handler, - std::vector &active_fe_indices) -{ - Assert (active_fe_indices.size() == dof_handler.get_tria().n_active_cells(), - ExcWrongSize (active_fe_indices.size(), - dof_handler.get_tria().n_active_cells())); - - typename DH::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (unsigned int index=0; cell!=endc; ++cell, ++index) - active_fe_indices[index] = cell->active_fe_index(); -} + dof_set.compress(); + } -template -void -DoFTools::get_subdomain_association (const DH &dof_handler, - std::vector &subdomain_association) -{ - // if the Triangulation is distributed, the - // only thing we can usefully ask is for - // its locally owned subdomain - Assert ((dynamic_cast*> - (&dof_handler.get_tria()) == 0), - ExcMessage ("For parallel::distributed::Triangulation objects and " - "associated DoF handler objects, asking for any subdomain other " - "than the locally owned one does not make sense.")); - - Assert(subdomain_association.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(subdomain_association.size(), - dof_handler.n_dofs())); - - // preset all values by an invalid value - std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(), - types::invalid_subdomain_id); - - std::vector local_dof_indices; - local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - - // pseudo-randomly assign variables - // which lie on the interface - // between subdomains to each of - // the two or more - bool coin_flip = true; - - // loop over all cells and record - // which subdomain a DoF belongs - // to. toss a coin in case it is on - // an interface - typename DH::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - Assert (cell->is_artificial() == false, - ExcMessage ("You can't call this function for meshes that " - "have artificial cells.")); - - const types::subdomain_id_t subdomain_id = cell->subdomain_id(); - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); - - // set subdomain ids. if dofs - // already have their values - // set then they must be on - // partition interfaces. in - // that case randomly assign - // them to either the previous - // association or the current - // one, where we take "random" - // to be "once this way once - // that way" - for (unsigned int i=0; i + void + extract_locally_relevant_dofs (const DH & dof_handler, + IndexSet & dof_set) + { + // collect all the locally owned dofs + dof_set = dof_handler.locally_owned_dofs(); + + // add the DoF on the adjacent ghost cells + // to the IndexSet, cache them in a + // set. need to check each dof manually + // because we can't be sure that the dof + // range of locally_owned_dofs is really + // contiguous. + std::vector dof_indices; + std::set global_dof_indices; + + typename DH::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_ghost()) + { + dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(dof_indices); + + for (std::vector::iterator it=dof_indices.begin(); + it!=dof_indices.end(); + ++it) + if (!dof_set.is_element(*it)) + global_dof_indices.insert(*it); + } - Assert (std::find (subdomain_association.begin(), - subdomain_association.end(), - types::invalid_subdomain_id) - == subdomain_association.end(), - ExcInternalError()); -} + dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); + dof_set.compress(); + } -template -unsigned int -DoFTools::count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain) -{ - std::vector subdomain_association (dof_handler.n_dofs()); - get_subdomain_association (dof_handler, subdomain_association); - return std::count (subdomain_association.begin(), - subdomain_association.end(), - subdomain); -} + template + void + extract_constant_modes (const DH &dof_handler, + const std::vector &component_select, + std::vector > &constant_modes) + { + const unsigned int n_components = dof_handler.get_fe().n_components(); + Assert (n_components == component_select.size(), + ExcDimensionMismatch(n_components, + component_select.size())); + std::vector localized_component (n_components, + numbers::invalid_unsigned_int); + unsigned int n_components_selected = 0; + for (unsigned int i=0; i dofs_by_component (dof_handler.n_locally_owned_dofs()); + internal::extract_dofs_by_component (dof_handler, component_select, false, + dofs_by_component); + unsigned int n_selected_dofs = 0; + for (unsigned int i=0; i(n_selected_dofs, + false)); + std::vector component_list (n_components, 0); + for (unsigned int d=0; d -IndexSet -DoFTools::dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain) -{ -#ifdef DEAL_II_USE_P4EST - // if the DoFHandler is distributed, the - // only thing we can usefully ask is for - // its locally owned subdomain - Assert ((dynamic_cast*> - (&dof_handler.get_tria()) == 0) - || - (subdomain == - dynamic_cast*> - (&dof_handler.get_tria())->locally_owned_subdomain()), - ExcMessage ("For parallel::distributed::Triangulation objects and " - "associated DoF handler objects, asking for any subdomain other " - "than the locally owned one does not make sense.")); -#endif + template + void + get_active_fe_indices (const DH &dof_handler, + std::vector &active_fe_indices) + { + Assert (active_fe_indices.size() == dof_handler.get_tria().n_active_cells(), + ExcWrongSize (active_fe_indices.size(), + dof_handler.get_tria().n_active_cells())); + + typename DH::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + active_fe_indices[index] = cell->active_fe_index(); + } - IndexSet index_set (dof_handler.n_dofs()); - std::vector local_dof_indices; - local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - typename DH::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - if ((cell->is_artificial() == false) - && - (cell->subdomain_id() == subdomain)) + template + void + get_subdomain_association (const DH &dof_handler, + std::vector &subdomain_association) + { + // if the Triangulation is distributed, the + // only thing we can usefully ask is for + // its locally owned subdomain + Assert ((dynamic_cast*> + (&dof_handler.get_tria()) == 0), + ExcMessage ("For parallel::distributed::Triangulation objects and " + "associated DoF handler objects, asking for any subdomain other " + "than the locally owned one does not make sense.")); + + Assert(subdomain_association.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(subdomain_association.size(), + dof_handler.n_dofs())); + + // preset all values by an invalid value + std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(), + types::invalid_subdomain_id); + + std::vector local_dof_indices; + local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); + + // pseudo-randomly assign variables + // which lie on the interface + // between subdomains to each of + // the two or more + bool coin_flip = true; + + // loop over all cells and record + // which subdomain a DoF belongs + // to. toss a coin in case it is on + // an interface + typename DH::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) { + Assert (cell->is_artificial() == false, + ExcMessage ("You can't call this function for meshes that " + "have artificial cells.")); + + const types::subdomain_id_t subdomain_id = cell->subdomain_id(); const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); - index_set.add_indices (local_dof_indices.begin(), - local_dof_indices.end()); + + // set subdomain ids. if dofs + // already have their values + // set then they must be on + // partition interfaces. in + // that case randomly assign + // them to either the previous + // association or the current + // one, where we take "random" + // to be "once this way once + // that way" + for (unsigned int i=0; i -void -DoFTools::count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain, - std::vector &n_dofs_on_subdomain) -{ - Assert (n_dofs_on_subdomain.size() == dof_handler.get_fe().n_components(), - ExcDimensionMismatch (n_dofs_on_subdomain.size(), - dof_handler.get_fe().n_components())); - std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0); - - // in debug mode, make sure that there are - // some cells at least with this subdomain - // id -#ifdef DEBUG + template + unsigned int + count_dofs_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain) { - bool found = false; - for (typename Triangulation::active_cell_iterator - cell=dof_handler.get_tria().begin_active(); - cell!=dof_handler.get_tria().end(); ++cell) - if (cell->subdomain_id() == subdomain) - { - found = true; - break; - } - Assert (found == true, - ExcMessage ("There are no cells for the given subdomain!")); + std::vector subdomain_association (dof_handler.n_dofs()); + get_subdomain_association (dof_handler, subdomain_association); + + return std::count (subdomain_association.begin(), + subdomain_association.end(), + subdomain); } + + + + template + IndexSet + dof_indices_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain) + { +#ifdef DEAL_II_USE_P4EST + // if the DoFHandler is distributed, the + // only thing we can usefully ask is for + // its locally owned subdomain + Assert ((dynamic_cast*> + (&dof_handler.get_tria()) == 0) + || + (subdomain == + dynamic_cast*> + (&dof_handler.get_tria())->locally_owned_subdomain()), + ExcMessage ("For parallel::distributed::Triangulation objects and " + "associated DoF handler objects, asking for any subdomain other " + "than the locally owned one does not make sense.")); #endif - std::vector subdomain_association (dof_handler.n_dofs()); - get_subdomain_association (dof_handler, subdomain_association); + IndexSet index_set (dof_handler.n_dofs()); - std::vector component_association (dof_handler.n_dofs()); - internal::extract_dofs_by_component (dof_handler, std::vector(), false, - component_association); + std::vector local_dof_indices; + local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - for (unsigned int c=0; c(c))) - ++n_dofs_on_subdomain[c]; - } -} + typename DH::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if ((cell->is_artificial() == false) + && + (cell->subdomain_id() == subdomain)) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + index_set.add_indices (local_dof_indices.begin(), + local_dof_indices.end()); + } + index_set.compress (); + return index_set; + } -namespace internal -{ - template + + template void - resolve_components (const FiniteElement&fe, - const std::vector &dofs_by_component, - const std::vector &target_component, - const bool only_once, - std::vector &dofs_per_component, - unsigned int &component) + count_dofs_with_subdomain_association (const DH &dof_handler, + const types::subdomain_id_t subdomain, + std::vector &n_dofs_on_subdomain) { - for (unsigned int b=0;b& base = fe.base_element(b); - // Dimension of base element - unsigned int d = base.n_components(); - - for (unsigned int m=0;m::active_cell_iterator + cell=dof_handler.get_tria().begin_active(); + cell!=dof_handler.get_tria().end(); ++cell) + if (cell->subdomain_id() == subdomain) { - if (base.n_base_elements() > 1) - resolve_components(base, dofs_by_component, target_component, - only_once, dofs_per_component, component); - else - { - for (unsigned int dd=0;dd subdomain_association (dof_handler.n_dofs()); + get_subdomain_association (dof_handler, subdomain_association); + + std::vector component_association (dof_handler.n_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), false, + component_association); + + for (unsigned int c=0; c(c))) + ++n_dofs_on_subdomain[c]; } } -} - -template -void -DoFTools::count_dofs_per_component ( - const DoFHandler& dof_handler, - std::vector& dofs_per_component, - bool only_once, - std::vector target_component) -{ - const FiniteElement& fe = dof_handler.get_fe(); - std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U); - // If the empty vector was given as - // default argument, set up this - // vector as identity. - if (target_component.size()==0) + namespace internal + { + template + void + resolve_components (const FiniteElement&fe, + const std::vector &dofs_by_component, + const std::vector &target_component, + const bool only_once, + std::vector &dofs_per_component, + unsigned int &component) { - target_component.resize(fe.n_components()); - for (unsigned int i=0; i& base = fe.base_element(b); + // Dimension of base element + unsigned int d = base.n_components(); + + for (unsigned int m=0;m 1) + resolve_components(base, dofs_by_component, target_component, + only_once, dofs_per_component, component); + else + { + for (unsigned int dd=0;dd + void + count_dofs_per_component ( + const DoFHandler& dof_handler, + std::vector& dofs_per_component, + bool only_once, + std::vector target_component) + { + const FiniteElement& fe = dof_handler.get_fe(); - AssertDimension (dofs_per_component.size(), n_target_components); + std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U); - // special case for only one - // component. treat this first - // since it does not require any - // computations - if (n_components == 1) - { - dofs_per_component[0] = dof_handler.n_locally_owned_dofs(); - return; - } + // If the empty vector was given as + // default argument, set up this + // vector as identity. + if (target_component.size()==0) + { + target_component.resize(fe.n_components()); + for (unsigned int i=0; i dofs_by_component (dof_handler.n_locally_owned_dofs()); - internal::extract_dofs_by_component (dof_handler, std::vector(), false, - dofs_by_component); - - // next count what we got - unsigned int component = 0; - internal::resolve_components(fe, dofs_by_component, target_component, - only_once, dofs_per_component, component); - Assert (n_components == component, ExcInternalError()); - - // finally sanity check. this is - // only valid if the finite element - // is actually primitive, so - // exclude other elements from this - Assert (!dof_handler.get_fe().is_primitive() - || - (std::accumulate (dofs_per_component.begin(), - dofs_per_component.end(), 0U) - == dof_handler.n_locally_owned_dofs()), - ExcInternalError()); + const unsigned int max_component + = *std::max_element (target_component.begin(), + target_component.end()); + const unsigned int n_target_components = max_component + 1; + const unsigned int n_components = fe.n_components(); + + AssertDimension (dofs_per_component.size(), n_target_components); + + // special case for only one + // component. treat this first + // since it does not require any + // computations + if (n_components == 1) + { + dofs_per_component[0] = dof_handler.n_locally_owned_dofs(); + return; + } + + + // otherwise determine the number + // of dofs in each component + // separately. do so in parallel + std::vector dofs_by_component (dof_handler.n_locally_owned_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), false, + dofs_by_component); + + // next count what we got + unsigned int component = 0; + internal::resolve_components(fe, dofs_by_component, target_component, + only_once, dofs_per_component, component); + Assert (n_components == component, ExcInternalError()); + + // finally sanity check. this is + // only valid if the finite element + // is actually primitive, so + // exclude other elements from this + Assert (!dof_handler.get_fe().is_primitive() + || + (std::accumulate (dofs_per_component.begin(), + dofs_per_component.end(), 0U) + == dof_handler.n_locally_owned_dofs()), + ExcInternalError()); // reduce information from all CPUs #ifdef DEAL_II_USE_P4EST #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - if (const parallel::distributed::Triangulation * tria - = (dynamic_cast*> - (&dof_handler.get_tria()))) - { - std::vector local_dof_count = dofs_per_component; + if (const parallel::distributed::Triangulation * tria + = (dynamic_cast*> + (&dof_handler.get_tria()))) + { + std::vector local_dof_count = dofs_per_component; - MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, - MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); - } + MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, + MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); + } #endif #endif -} + } -template -void -DoFTools:: -count_dofs_per_block (const DoFHandler& dof_handler, - std::vector &dofs_per_block, - std::vector target_block) -{ - const FiniteElement& fe = dof_handler.get_fe(); + template + void + count_dofs_per_block (const DoFHandler& dof_handler, + std::vector &dofs_per_block, + std::vector target_block) + { + const FiniteElement& fe = dof_handler.get_fe(); - std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U); + std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U); - // If the empty vector was given as - // default argument, set up this - // vector as identity. - if (target_block.size()==0) - { - target_block.resize(fe.n_blocks()); - for (unsigned int i=0; i dofs_by_block (dof_handler.n_locally_owned_dofs()); - internal::extract_dofs_by_component (dof_handler, std::vector(), - true, dofs_by_block); - - // next count what we got - for (unsigned int block=0; block dofs_by_block (dof_handler.n_locally_owned_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), + true, dofs_by_block); + + // next count what we got + for (unsigned int block=0; block * tria - = (dynamic_cast*> - (&dof_handler.get_tria()))) - { - std::vector local_dof_count = dofs_per_block; - MPI_Allreduce ( &local_dof_count[0], &dofs_per_block[0], n_target_blocks, - MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); - } + // if we are working on a parallel + // mesh, we now need to collect + // this information from all + // processors + if (const parallel::distributed::Triangulation * tria + = (dynamic_cast*> + (&dof_handler.get_tria()))) + { + std::vector local_dof_count = dofs_per_block; + MPI_Allreduce ( &local_dof_count[0], &dofs_per_block[0], n_target_blocks, + MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); + } #endif #endif -} + } -template -void -DoFTools:: -count_dofs_per_component (const DoFHandler &dof_handler, - std::vector &dofs_per_component, - std::vector target_component) -{ - count_dofs_per_component (dof_handler, dofs_per_component, - false, target_component); -} + template + void + count_dofs_per_component (const DoFHandler &dof_handler, + std::vector &dofs_per_component, + std::vector target_component) + { + count_dofs_per_component (dof_handler, dofs_per_component, + false, target_component); + } -namespace internal -{ - namespace + namespace internal { - /** - * This is a helper function that - * is used in the computation of - * integrid constraints. See the - * function for a thorough - * description of how it works. - */ - template - unsigned int - compute_intergrid_weights_1 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const dealii::DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - std::vector > &weights, - std::vector &weight_mapping) + namespace { - // aliases to the finite elements - // used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(), - &fine_fe = fine_grid.get_fe(); - - // global numbers of dofs - const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), - n_fine_dofs = fine_grid.n_dofs(); - - // local numbers of dofs - const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell; - - // alias the number of dofs per - // cell belonging to the - // coarse_component which is to be - // the restriction of the fine - // grid: - const unsigned int coarse_dofs_per_cell_component - = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell; - - - // Try to find out whether the - // grids stem from the same coarse - // grid. This is a rather crude - // test, but better than nothing - Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0), - dealii::DoFTools::ExcGridsDontMatch()); - - // check whether the map correlates - // the right objects - Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, - dealii::DoFTools::ExcGridsDontMatch ()); - Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, - dealii::DoFTools::ExcGridsDontMatch ()); - - - // check whether component numbers - // are valid - Assert (coarse_component < coarse_fe.n_components(), - dealii::DoFTools::ExcInvalidComponent (coarse_component, coarse_fe.n_components())); - Assert (fine_component < fine_fe.n_components(), - dealii::DoFTools::ExcInvalidComponent (fine_component, fine_fe.n_components())); - // check whether respective finite - // elements are equal - Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first) - == - fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first), - dealii::DoFTools::ExcFiniteElementsDontMatch()); + /** + * This is a helper function that + * is used in the computation of + * integrid constraints. See the + * function for a thorough + * description of how it works. + */ + template + unsigned int + compute_intergrid_weights_1 ( + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const dealii::DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + std::vector > &weights, + std::vector &weight_mapping) + { + // aliases to the finite elements + // used by the dof handlers: + const FiniteElement &coarse_fe = coarse_grid.get_fe(), + &fine_fe = fine_grid.get_fe(); + + // global numbers of dofs + const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), + n_fine_dofs = fine_grid.n_dofs(); + + // local numbers of dofs + const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell; + + // alias the number of dofs per + // cell belonging to the + // coarse_component which is to be + // the restriction of the fine + // grid: + const unsigned int coarse_dofs_per_cell_component + = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell; + + + // Try to find out whether the + // grids stem from the same coarse + // grid. This is a rather crude + // test, but better than nothing + Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0), + ExcGridsDontMatch()); + + // check whether the map correlates + // the right objects + Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, + ExcGridsDontMatch ()); + Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, + ExcGridsDontMatch ()); + + + // check whether component numbers + // are valid + Assert (coarse_component < coarse_fe.n_components(), + ExcInvalidComponent (coarse_component, coarse_fe.n_components())); + Assert (fine_component < fine_fe.n_components(), + ExcInvalidComponent (fine_component, fine_fe.n_components())); + // check whether respective finite + // elements are equal + Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first) + == + fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first), + ExcFiniteElementsDontMatch()); #ifdef DEBUG - // if in debug mode, check whether - // the coarse grid is indeed - // coarser everywhere than the fine - // grid - for (typename dealii::DoFHandler::active_cell_iterator - cell=coarse_grid.begin_active(); - cell != coarse_grid.end(); ++cell) - Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(), - dealii::DoFTools::ExcGridNotCoarser()); + // if in debug mode, check whether + // the coarse grid is indeed + // coarser everywhere than the fine + // grid + for (typename dealii::DoFHandler::active_cell_iterator + cell=coarse_grid.begin_active(); + cell != coarse_grid.end(); ++cell) + Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(), + ExcGridNotCoarser()); #endif @@ -4773,999 +4752,1000 @@ namespace internal - // set up vectors of cell-local - // data; each vector represents one - // degree of freedom of the - // coarse-grid variable in the - // fine-grid element - std::vector > - parameter_dofs (coarse_dofs_per_cell_component, - dealii::Vector(fine_dofs_per_cell)); - // for each coarse dof: find its - // position within the fine element - // and set this value to one in the - // respective vector (all other values - // are zero by construction) - for (unsigned int local_coarse_dof=0; - local_coarse_dof > + parameter_dofs (coarse_dofs_per_cell_component, + dealii::Vector(fine_dofs_per_cell)); + // for each coarse dof: find its + // position within the fine element + // and set this value to one in the + // respective vector (all other values + // are zero by construction) + for (unsigned int local_coarse_dof=0; + local_coarse_dof dof_is_interesting (fine_grid.n_dofs(), false); - std::vector local_dof_indices (fine_fe.dofs_per_cell); - - for (typename dealii::DoFHandler::active_cell_iterator - cell=fine_grid.begin_active(); - cell!=fine_grid.end(); ++cell) - { - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i dof_is_interesting (fine_grid.n_dofs(), false); + std::vector local_dof_indices (fine_fe.dofs_per_cell); + + for (typename dealii::DoFHandler::active_cell_iterator + cell=fine_grid.begin_active(); + cell!=fine_grid.end(); ++cell) + { + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i local_dof_indices(fine_fe.dofs_per_cell); - unsigned int next_free_index=0; - for (typename dealii::DoFHandler::active_cell_iterator - cell=fine_grid.begin_active(); - cell != fine_grid.end(); ++cell) - { - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i local_dof_indices(fine_fe.dofs_per_cell); + unsigned int next_free_index=0; + for (typename dealii::DoFHandler::active_cell_iterator + cell=fine_grid.begin_active(); + cell != fine_grid.end(); ++cell) + { + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i1) && (sum==0)), ExcInternalError()); - }; + for (unsigned int col=0; col1) && (sum==0)), ExcInternalError()); + }; #endif - return n_parameters_on_fine_grid; - } + return n_parameters_on_fine_grid; + } - /** - * This is a function that is - * called by the _2 function and - * that operates on a range of - * cells only. It is used to - * split up the whole range of - * cells into chunks which are - * then worked on in parallel, if - * multithreading is available. - */ - template - void - compute_intergrid_weights_3 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap > &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights, - const typename dealii::DoFHandler::active_cell_iterator &begin, - const typename dealii::DoFHandler::active_cell_iterator &end) - { - // aliases to the finite elements - // used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(); - - // for each cell on the parameter grid: - // find out which degrees of freedom on the - // fine grid correspond in which way to - // the degrees of freedom on the parameter - // grid - // - // since for continuous FEs some - // dofs exist on more than one - // cell, we have to track which - // ones were already visited. the - // problem is that if we visit a - // dof first on one cell and - // compute its weight with respect - // to some global dofs to be - // non-zero, and later visit the - // dof again on another cell and - // (since we are on another cell) - // recompute the weights with - // respect to the same dofs as - // above to be zero now, we have to - // preserve them. we therefore - // overwrite all weights if they - // are nonzero and do not enforce - // zero weights since that might be - // only due to the fact that we are - // on another cell. - // - // example: - // coarse grid - // | | | - // *-----*-----* - // | cell|cell | - // | 1 | 2 | - // | | | - // 0-----1-----* - // - // fine grid - // | | | | | - // *--*--*--*--* - // | | | | | - // *--*--*--*--* - // | | | | | - // *--x--y--*--* - // - // when on cell 1, we compute the - // weights of dof 'x' to be 1/2 - // from parameter dofs 0 and 1, - // respectively. however, when - // later we are on cell 2, we again - // compute the prolongation of - // shape function 1 restricted to - // cell 2 to the globla grid and - // find that the weight of global - // dof 'x' now is zero. however, we - // should not overwrite the old - // value. - // - // we therefore always only set - // nonzero values. why adding up is - // not useful: dof 'y' would get - // weight 1 from parameter dof 1 on - // both cells 1 and 2, but the - // correct weight is nevertheless - // only 1. - - // vector to hold the representation of - // a single degree of freedom on the - // coarse grid (for the selected fe) - // on the fine grid - const unsigned int n_fine_dofs = weight_mapping.size(); - dealii::Vector global_parameter_representation (n_fine_dofs); - - typename dealii::DoFHandler::active_cell_iterator cell; - std::vector parameter_dof_indices (coarse_fe.dofs_per_cell); - - for (cell=begin; cell!=end; ++cell) - { - // get the global indices of the - // parameter dofs on this parameter - // grid cell - cell->get_dof_indices (parameter_dof_indices); - - // loop over all dofs on this - // cell and check whether they - // are interesting for us - for (unsigned int local_dof=0; - local_dof - set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof], - global_parameter_representation); - // now that we've got the global - // representation of each parameter - // dof, we've only got to clobber the - // non-zero entries in that vector and - // store the result - // - // what we have learned: if entry @p{i} - // of the global vector holds the value - // @p{v[i]}, then this is the weight with - // which the present dof contributes - // to @p{i}. there may be several such - // @p{i}s and their weights' sum should - // be one. Then, @p{v[i]} should be - // equal to @p{\sum_j w_{ij} p[j]} with - // @p{p[j]} be the values of the degrees - // of freedom on the coarse grid. we - // can thus compute constraints which - // link the degrees of freedom @p{v[i]} - // on the fine grid to those on the - // coarse grid, @p{p[j]}. Now to use - // these as real constraints, rather - // than as additional equations, we - // have to identify representants - // among the @p{i} for each @p{j}. this will - // be done by simply taking the first - // @p{i} for which @p{w_{ij}==1}. - // - // guard modification of - // the weights array by a - // Mutex. since it should - // happen rather rarely - // that there are several - // threads operating on - // different intergrid - // weights, have only one - // mutex for all of them - static Threads::ThreadMutex mutex; - Threads::ThreadMutex::ScopedLock lock (mutex); - for (unsigned int i=0; i + void + compute_intergrid_weights_3 ( + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap > &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights, + const typename dealii::DoFHandler::active_cell_iterator &begin, + const typename dealii::DoFHandler::active_cell_iterator &end) + { + // aliases to the finite elements + // used by the dof handlers: + const FiniteElement &coarse_fe = coarse_grid.get_fe(); + + // for each cell on the parameter grid: + // find out which degrees of freedom on the + // fine grid correspond in which way to + // the degrees of freedom on the parameter + // grid + // + // since for continuous FEs some + // dofs exist on more than one + // cell, we have to track which + // ones were already visited. the + // problem is that if we visit a + // dof first on one cell and + // compute its weight with respect + // to some global dofs to be + // non-zero, and later visit the + // dof again on another cell and + // (since we are on another cell) + // recompute the weights with + // respect to the same dofs as + // above to be zero now, we have to + // preserve them. we therefore + // overwrite all weights if they + // are nonzero and do not enforce + // zero weights since that might be + // only due to the fact that we are + // on another cell. + // + // example: + // coarse grid + // | | | + // *-----*-----* + // | cell|cell | + // | 1 | 2 | + // | | | + // 0-----1-----* + // + // fine grid + // | | | | | + // *--*--*--*--* + // | | | | | + // *--*--*--*--* + // | | | | | + // *--x--y--*--* + // + // when on cell 1, we compute the + // weights of dof 'x' to be 1/2 + // from parameter dofs 0 and 1, + // respectively. however, when + // later we are on cell 2, we again + // compute the prolongation of + // shape function 1 restricted to + // cell 2 to the globla grid and + // find that the weight of global + // dof 'x' now is zero. however, we + // should not overwrite the old + // value. + // + // we therefore always only set + // nonzero values. why adding up is + // not useful: dof 'y' would get + // weight 1 from parameter dof 1 on + // both cells 1 and 2, but the + // correct weight is nevertheless + // only 1. + + // vector to hold the representation of + // a single degree of freedom on the + // coarse grid (for the selected fe) + // on the fine grid + const unsigned int n_fine_dofs = weight_mapping.size(); + dealii::Vector global_parameter_representation (n_fine_dofs); + + typename dealii::DoFHandler::active_cell_iterator cell; + std::vector parameter_dof_indices (coarse_fe.dofs_per_cell); + + for (cell=begin; cell!=end; ++cell) + { + // get the global indices of the + // parameter dofs on this parameter + // grid cell + cell->get_dof_indices (parameter_dof_indices); + + // loop over all dofs on this + // cell and check whether they + // are interesting for us + for (unsigned int local_dof=0; + local_dof + set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof], + global_parameter_representation); + // now that we've got the global + // representation of each parameter + // dof, we've only got to clobber the + // non-zero entries in that vector and + // store the result + // + // what we have learned: if entry @p{i} + // of the global vector holds the value + // @p{v[i]}, then this is the weight with + // which the present dof contributes + // to @p{i}. there may be several such + // @p{i}s and their weights' sum should + // be one. Then, @p{v[i]} should be + // equal to @p{\sum_j w_{ij} p[j]} with + // @p{p[j]} be the values of the degrees + // of freedom on the coarse grid. we + // can thus compute constraints which + // link the degrees of freedom @p{v[i]} + // on the fine grid to those on the + // coarse grid, @p{p[j]}. Now to use + // these as real constraints, rather + // than as additional equations, we + // have to identify representants + // among the @p{i} for each @p{j}. this will + // be done by simply taking the first + // @p{i} for which @p{w_{ij}==1}. + // + // guard modification of + // the weights array by a + // Mutex. since it should + // happen rather rarely + // that there are several + // threads operating on + // different intergrid + // weights, have only one + // mutex for all of them + static Threads::ThreadMutex mutex; + Threads::ThreadMutex::ScopedLock lock (mutex); + for (unsigned int i=0; i - void - compute_intergrid_weights_2 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap > &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights) - { - // simply distribute the range of - // cells to different threads - typedef typename dealii::DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > - cell_intervals = Threads::split_range (coarse_grid.begin_active(), - coarse_grid.end(), - multithread_info.n_default_threads); + /** + * This is a helper function that + * is used in the computation of + * integrid constraints. See the + * function for a thorough + * description of how it works. + */ + template + void + compute_intergrid_weights_2 ( + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap > &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights) + { + // simply distribute the range of + // cells to different threads + typedef typename dealii::DoFHandler::active_cell_iterator active_cell_iterator; + std::vector > + cell_intervals = Threads::split_range (coarse_grid.begin_active(), + coarse_grid.end(), + multithread_info.n_default_threads); //TODO: use WorkStream here - Threads::TaskGroup<> tasks; - void (*fun_ptr) (const dealii::DoFHandler &, - const unsigned int , - const InterGridMap > &, - const std::vector > &, - const std::vector &, - std::vector > &, - const typename dealii::DoFHandler::active_cell_iterator &, - const typename dealii::DoFHandler::active_cell_iterator &) - = &compute_intergrid_weights_3; - for (unsigned int i=0; i tasks; + void (*fun_ptr) (const dealii::DoFHandler &, + const unsigned int , + const InterGridMap > &, + const std::vector > &, + const std::vector &, + std::vector > &, + const typename dealii::DoFHandler::active_cell_iterator &, + const typename dealii::DoFHandler::active_cell_iterator &) + = &compute_intergrid_weights_3; + for (unsigned int i=0; i -void -DoFTools::compute_intergrid_constraints ( - const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - ConstraintMatrix &constraints) -{ - // store the weights with which a dof - // on the parameter grid contributes - // to a dof on the fine grid. see the - // long doc below for more info - // - // allocate as many rows as there are - // parameter dofs on the coarse grid - // and as many columns as there are - // parameter dofs on the fine grid. - // - // weight_mapping is used to map the - // global (fine grid) parameter dof - // indices to the columns - // - // in the original implementation, - // the weights array was actually - // of FullMatrix type. this - // wasted huge amounts of memory, - // but was fast. nonetheless, since - // the memory consumption was - // quadratic in the number of - // degrees of freedom, this was not - // very practical, so we now use a - // vector of rows of the matrix, - // and in each row a vector of - // pairs (colnum,value). this seems - // like the best tradeoff between - // memory and speed, as it is now - // linear in memory and still fast - // enough. - // - // to save some memory and since - // the weights are usually - // (negative) powers of 2, we - // choose the value type of the - // matrix to be @p{float} rather - // than @p{double}. - std::vector > weights; - - // this is this mapping. there is one - // entry for each dof on the fine grid; - // if it is a parameter dof, then its - // value is the column in weights for - // that parameter dof, if it is any - // other dof, then its value is -1, - // indicating an error - std::vector weight_mapping; - - const unsigned int n_parameters_on_fine_grid - = internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, - fine_grid, fine_component, - coarse_to_fine_grid_map, - weights, weight_mapping); - - // global numbers of dofs - const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), - n_fine_dofs = fine_grid.n_dofs(); - - - // get an array in which we store - // which dof on the coarse grid is - // a parameter and which is not - std::vector coarse_dof_is_parameter (coarse_grid.n_dofs()); - if (true) - { - std::vector mask (coarse_grid.get_fe().n_components(), - false); - mask[coarse_component] = true; - extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); - }; - - // now we know that the weights in - // each row constitute a - // constraint. enter this into the - // constraints object - // - // first task: for each parameter - // dof on the parameter grid, find - // a representant on the fine, - // global grid. this is possible - // since we use conforming finite - // element. we take this - // representant to be the first - // element in this row with weight - // identical to one. the - // representant will become an - // unconstrained degree of freedom, - // while all others will be - // constrained to this dof (and - // possibly others) - std::vector representants(n_coarse_dofs, -1); - for (unsigned int parameter_dof=0; parameter_dof 0, ExcInternalError()); - - // find the column where the - // representant is mentioned - std::map::const_iterator i = weights[parameter_dof].begin(); - for (; i!=weights[parameter_dof].end(); ++i) - if (i->second == 1) - break; - Assert (i!=weights[parameter_dof].end(), ExcInternalError()); - const unsigned int column = i->first; - - // now we know in which column of - // weights the representant is, but - // we don't know its global index. get - // it using the inverse operation of - // the weight_mapping - unsigned int global_dof=0; - for (; global_dof(column)) - break; - Assert (global_dof < weight_mapping.size(), ExcInternalError()); - - // now enter the representants global - // index into our list - representants[parameter_dof] = global_dof; - } - else + template + void + compute_intergrid_constraints ( + const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + ConstraintMatrix &constraints) + { + // store the weights with which a dof + // on the parameter grid contributes + // to a dof on the fine grid. see the + // long doc below for more info + // + // allocate as many rows as there are + // parameter dofs on the coarse grid + // and as many columns as there are + // parameter dofs on the fine grid. + // + // weight_mapping is used to map the + // global (fine grid) parameter dof + // indices to the columns + // + // in the original implementation, + // the weights array was actually + // of FullMatrix type. this + // wasted huge amounts of memory, + // but was fast. nonetheless, since + // the memory consumption was + // quadratic in the number of + // degrees of freedom, this was not + // very practical, so we now use a + // vector of rows of the matrix, + // and in each row a vector of + // pairs (colnum,value). this seems + // like the best tradeoff between + // memory and speed, as it is now + // linear in memory and still fast + // enough. + // + // to save some memory and since + // the weights are usually + // (negative) powers of 2, we + // choose the value type of the + // matrix to be @p{float} rather + // than @p{double}. + std::vector > weights; + + // this is this mapping. there is one + // entry for each dof on the fine grid; + // if it is a parameter dof, then its + // value is the column in weights for + // that parameter dof, if it is any + // other dof, then its value is -1, + // indicating an error + std::vector weight_mapping; + + const unsigned int n_parameters_on_fine_grid + = internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, + fine_grid, fine_component, + coarse_to_fine_grid_map, + weights, weight_mapping); + + // global numbers of dofs + const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), + n_fine_dofs = fine_grid.n_dofs(); + + + // get an array in which we store + // which dof on the coarse grid is + // a parameter and which is not + std::vector coarse_dof_is_parameter (coarse_grid.n_dofs()); + if (true) { - // consistency check: if this - // is no parameter dof on the - // coarse grid, then the - // respective row must be - // empty! - Assert (weights[parameter_dof].size() == 0, ExcInternalError()); + std::vector mask (coarse_grid.get_fe().n_components(), + false); + mask[coarse_component] = true; + extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); }; + // now we know that the weights in + // each row constitute a + // constraint. enter this into the + // constraints object + // + // first task: for each parameter + // dof on the parameter grid, find + // a representant on the fine, + // global grid. this is possible + // since we use conforming finite + // element. we take this + // representant to be the first + // element in this row with weight + // identical to one. the + // representant will become an + // unconstrained degree of freedom, + // while all others will be + // constrained to this dof (and + // possibly others) + std::vector representants(n_coarse_dofs, -1); + for (unsigned int parameter_dof=0; parameter_dof 0, ExcInternalError()); + + // find the column where the + // representant is mentioned + std::map::const_iterator i = weights[parameter_dof].begin(); + for (; i!=weights[parameter_dof].end(); ++i) + if (i->second == 1) + break; + Assert (i!=weights[parameter_dof].end(), ExcInternalError()); + const unsigned int column = i->first; + + // now we know in which column of + // weights the representant is, but + // we don't know its global index. get + // it using the inverse operation of + // the weight_mapping + unsigned int global_dof=0; + for (; global_dof(column)) + break; + Assert (global_dof < weight_mapping.size(), ExcInternalError()); + // now enter the representants global + // index into our list + representants[parameter_dof] = global_dof; + } + else + { + // consistency check: if this + // is no parameter dof on the + // coarse grid, then the + // respective row must be + // empty! + Assert (weights[parameter_dof].size() == 0, ExcInternalError()); + }; - // note for people that want to - // optimize this function: the - // largest part of the computing - // time is spent in the following, - // rather innocent block of - // code. basically, it must be the - // ConstraintMatrix::add_entry call - // which takes the bulk of the - // time, but it is not known to the - // author how to make it faster... - std::vector > constraint_line; - for (unsigned int global_dof=0; global_dof > constraint_line; + for (unsigned int global_dof=0; global_dof 0, ExcInternalError()); - std::map::const_iterator - col_entry = weights[0].end(); - for (; first_used_rowsecond == 1) && - (representants[first_used_row] == static_cast(global_dof))) - // dof unconstrained or - // constrained to itself - // (in case this cell is - // mapped to itself, rather - // than to children of - // itself) - continue; - } + const unsigned int col = weight_mapping[global_dof]; + Assert (col < n_parameters_on_fine_grid, ExcInternalError()); + unsigned int first_used_row=0; - // otherwise enter all constraints - constraints.add_line (global_dof); - - constraint_line.clear (); - for (unsigned int row=first_used_row; row::const_iterator - j = weights[row].find(col); - if ((j != weights[row].end()) && (j->second != 0)) - constraint_line.push_back (std::make_pair(representants[row], - j->second)); - }; + Assert (weights.size() > 0, ExcInternalError()); + std::map::const_iterator + col_entry = weights[0].end(); + for (; first_used_rowsecond == 1) && + (representants[first_used_row] == static_cast(global_dof))) + // dof unconstrained or + // constrained to itself + // (in case this cell is + // mapped to itself, rather + // than to children of + // itself) + continue; + } + // otherwise enter all constraints + constraints.add_line (global_dof); -template -void -DoFTools:: -compute_intergrid_transfer_representation ( - const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - std::vector > &transfer_representation) -{ - // store the weights with which a dof - // on the parameter grid contributes - // to a dof on the fine grid. see the - // long doc below for more info - // - // allocate as many rows as there are - // parameter dofs on the coarse grid - // and as many columns as there are - // parameter dofs on the fine grid. - // - // weight_mapping is used to map the - // global (fine grid) parameter dof - // indices to the columns - // - // in the original implementation, - // the weights array was actually - // of FullMatrix type. this - // wasted huge amounts of memory, - // but was fast. nonetheless, since - // the memory consumption was - // quadratic in the number of - // degrees of freedom, this was not - // very practical, so we now use a - // vector of rows of the matrix, - // and in each row a vector of - // pairs (colnum,value). this seems - // like the best tradeoff between - // memory and speed, as it is now - // linear in memory and still fast - // enough. - // - // to save some memory and since - // the weights are usually - // (negative) powers of 2, we - // choose the value type of the - // matrix to be @p{float} rather - // than @p{double}. - std::vector > weights; - - // this is this mapping. there is one - // entry for each dof on the fine grid; - // if it is a parameter dof, then its - // value is the column in weights for - // that parameter dof, if it is any - // other dof, then its value is -1, - // indicating an error - std::vector weight_mapping; - - internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, - fine_grid, fine_component, - coarse_to_fine_grid_map, - weights, weight_mapping); - - // now compute the requested - // representation - const unsigned int n_global_parm_dofs - = std::count_if (weight_mapping.begin(), weight_mapping.end(), - std::bind2nd (std::not_equal_to (), -1)); - - // first construct the inverse - // mapping of weight_mapping - std::vector inverse_weight_mapping (n_global_parm_dofs, - DoFHandler::invalid_dof_index); - for (unsigned int i=0; i::invalid_dof_index), - ExcInternalError()); + constraint_line.clear (); + for (unsigned int row=first_used_row; row::const_iterator + j = weights[row].find(col); + if ((j != weights[row].end()) && (j->second != 0)) + constraint_line.push_back (std::make_pair(representants[row], + j->second)); + }; - inverse_weight_mapping[parameter_dof] = i; + constraints.add_entries (global_dof, constraint_line); }; - }; + } - // next copy over weights array - // and replace respective - // numbers - const unsigned int n_rows = weight_mapping.size(); - transfer_representation.clear (); - transfer_representation.resize (n_rows); - const unsigned int n_coarse_dofs = coarse_grid.n_dofs(); - for (unsigned int i=0; i::const_iterator j = weights[i].begin(); - for (; j!=weights[i].end(); ++j) - { - const unsigned int p = inverse_weight_mapping[j->first]; - Assert (p + void + compute_intergrid_transfer_representation ( + const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + std::vector > &transfer_representation) + { + // store the weights with which a dof + // on the parameter grid contributes + // to a dof on the fine grid. see the + // long doc below for more info + // + // allocate as many rows as there are + // parameter dofs on the coarse grid + // and as many columns as there are + // parameter dofs on the fine grid. + // + // weight_mapping is used to map the + // global (fine grid) parameter dof + // indices to the columns + // + // in the original implementation, + // the weights array was actually + // of FullMatrix type. this + // wasted huge amounts of memory, + // but was fast. nonetheless, since + // the memory consumption was + // quadratic in the number of + // degrees of freedom, this was not + // very practical, so we now use a + // vector of rows of the matrix, + // and in each row a vector of + // pairs (colnum,value). this seems + // like the best tradeoff between + // memory and speed, as it is now + // linear in memory and still fast + // enough. + // + // to save some memory and since + // the weights are usually + // (negative) powers of 2, we + // choose the value type of the + // matrix to be @p{float} rather + // than @p{double}. + std::vector > weights; + + // this is this mapping. there is one + // entry for each dof on the fine grid; + // if it is a parameter dof, then its + // value is the column in weights for + // that parameter dof, if it is any + // other dof, then its value is -1, + // indicating an error + std::vector weight_mapping; + + internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, + fine_grid, fine_component, + coarse_to_fine_grid_map, + weights, weight_mapping); + + // now compute the requested + // representation + const unsigned int n_global_parm_dofs + = std::count_if (weight_mapping.begin(), weight_mapping.end(), + std::bind2nd (std::not_equal_to (), -1)); + + // first construct the inverse + // mapping of weight_mapping + std::vector inverse_weight_mapping (n_global_parm_dofs, + DoFHandler::invalid_dof_index); + for (unsigned int i=0; i::invalid_dof_index), + ExcInternalError()); - transfer_representation[p][i] = j->second; - }; - }; -} + inverse_weight_mapping[parameter_dof] = i; + }; + }; + // next copy over weights array + // and replace respective + // numbers + const unsigned int n_rows = weight_mapping.size(); + transfer_representation.clear (); + transfer_representation.resize (n_rows); -template -void -DoFTools::map_dof_to_boundary_indices (const DH &dof_handler, - std::vector &mapping) -{ - Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); - - mapping.clear (); - mapping.insert (mapping.end(), dof_handler.n_dofs(), - DH::invalid_dof_index); - - std::vector dofs_on_face; - dofs_on_face.reserve (max_dofs_per_face(dof_handler)); - unsigned int next_boundary_index = 0; - - // now loop over all cells and - // check whether their faces are at - // the boundary. note that we need - // not take special care of single - // lines being at the boundary - // (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every isolated boundary - // line is also part of a boundary - // face which we will be visiting - // sooner or later - typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->at_boundary(f)) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_face, - cell->active_fe_index()); - for (unsigned int i=0; i::const_iterator j = weights[i].begin(); + for (; j!=weights[i].end(); ++j) + { + const unsigned int p = inverse_weight_mapping[j->first]; + Assert (psecond; + }; + }; + } -template -void DoFTools::map_dof_to_boundary_indices ( - const DH &dof_handler, - const std::set &boundary_indicators, - std::vector &mapping) -{ - Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); - Assert (boundary_indicators.find (255) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); - - mapping.clear (); - mapping.insert (mapping.end(), dof_handler.n_dofs(), - DH::invalid_dof_index); - - // return if there is nothing to do - if (boundary_indicators.size() == 0) - return; - - std::vector dofs_on_face; - dofs_on_face.reserve (max_dofs_per_face(dof_handler)); - unsigned int next_boundary_index = 0; - - typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (boundary_indicators.find (cell->face(f)->boundary_indicator()) != - boundary_indicators.end()) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index()); - for (unsigned int i=0; i + void + map_dof_to_boundary_indices (const DH &dof_handler, + std::vector &mapping) + { + Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); + + mapping.clear (); + mapping.insert (mapping.end(), dof_handler.n_dofs(), + DH::invalid_dof_index); + + std::vector dofs_on_face; + dofs_on_face.reserve (max_dofs_per_face(dof_handler)); + unsigned int next_boundary_index = 0; + + // now loop over all cells and + // check whether their faces are at + // the boundary. note that we need + // not take special care of single + // lines being at the boundary + // (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every isolated boundary + // line is also part of a boundary + // face which we will be visiting + // sooner or later + typename DH::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_face, + cell->active_fe_index()); + for (unsigned int i=0; i -void -DoFTools::map_dofs_to_support_points (const Mapping &mapping, - const DoFHandler &dof_handler, - std::vector > &support_points) -{ - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - - // check whether fe has support - // points - Assert (dof_handler.get_fe().has_support_points(), - ExcFEHasNoSupportPoints()); - Assert (support_points.size() == dof_handler.n_dofs(), - ExcWrongSize (support_points.size(), dof_handler.n_dofs())); - - // now loop over all cells and - // enquire the support points on - // each of these. use a dummy - // quadrature formula where the - // quadrature points are located at - // the unit support points to - // enquire the location of the - // support points in real space - // - // the weights of the quadrature - // rule are set to invalid values - // by the used constructor. - Quadrature q_dummy(dof_handler.get_fe().get_unit_support_points()); - FEValues fe_values (mapping, dof_handler.get_fe(), - q_dummy, update_quadrature_points); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - std::vector local_dof_indices (dofs_per_cell); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell->get_dof_indices (local_dof_indices); - const std::vector > & points - = fe_values.get_quadrature_points (); - for (unsigned int i=0; i + void map_dof_to_boundary_indices ( + const DH &dof_handler, + const std::set &boundary_indicators, + std::vector &mapping) + { + Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); + Assert (boundary_indicators.find (255) == boundary_indicators.end(), + ExcInvalidBoundaryIndicator()); + mapping.clear (); + mapping.insert (mapping.end(), dof_handler.n_dofs(), + DH::invalid_dof_index); -template -void -DoFTools::convert_couplings_to_blocks ( - const DoFHandler& dof_handler, - const Table<2, Coupling>& table, - std::vector >& tables_by_block) -{ - const FiniteElement& fe = dof_handler.get_fe(); - const unsigned int nb = fe.n_blocks(); + // return if there is nothing to do + if (boundary_indicators.size() == 0) + return; - tables_by_block.resize(1); - tables_by_block[0].reinit(nb, nb); - tables_by_block[0].fill(none); + std::vector dofs_on_face; + dofs_on_face.reserve (max_dofs_per_face(dof_handler)); + unsigned int next_boundary_index = 0; - for (unsigned int i=0;i::faces_per_cell; ++f) + if (boundary_indicators.find (cell->face(f)->boundary_indicator()) != + boundary_indicators.end()) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index()); + for (unsigned int i=0; i -void -DoFTools::convert_couplings_to_blocks ( - const hp::DoFHandler& dof_handler, - const Table<2, Coupling>& table, - std::vector >& tables_by_block) -{ - const hp::FECollection& fe_collection = dof_handler.get_fe(); - tables_by_block.resize(fe_collection.size()); - for (unsigned int f=0;f& fe = fe_collection[f]; - const unsigned int nb = fe.n_blocks(); - tables_by_block[f].reinit(nb, nb); - tables_by_block[f].fill(none); - for (unsigned int i=0;i + void + map_dofs_to_support_points (const Mapping &mapping, + const DoFHandler &dof_handler, + std::vector > &support_points) + { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + + // check whether fe has support + // points + Assert (dof_handler.get_fe().has_support_points(), + ExcFEHasNoSupportPoints()); + Assert (support_points.size() == dof_handler.n_dofs(), + ExcWrongSize (support_points.size(), dof_handler.n_dofs())); + + // now loop over all cells and + // enquire the support points on + // each of these. use a dummy + // quadrature formula where the + // quadrature points are located at + // the unit support points to + // enquire the location of the + // support points in real space + // + // the weights of the quadrature + // rule are set to invalid values + // by the used constructor. + Quadrature q_dummy(dof_handler.get_fe().get_unit_support_points()); + FEValues fe_values (mapping, dof_handler.get_fe(), + q_dummy, update_quadrature_points); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + std::vector local_dof_indices (dofs_per_cell); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell->get_dof_indices (local_dof_indices); + const std::vector > & points + = fe_values.get_quadrature_points (); + for (unsigned int i=0; i + void + convert_couplings_to_blocks ( + const DoFHandler& dof_handler, + const Table<2, Coupling>& table, + std::vector >& tables_by_block) + { + const FiniteElement& fe = dof_handler.get_fe(); + const unsigned int nb = fe.n_blocks(); -template class DH> -void -DoFTools::make_zero_boundary_constraints (const DH &dof, - ConstraintMatrix &zero_boundary_constraints, - const std::vector &component_mask_) -{ - Assert ((component_mask_.size() == 0) || - (component_mask_.size() == dof.get_fe().n_components()), - ExcMessage ("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.")); - - const unsigned int n_components = DoFTools::n_components(dof); - - // set the component mask to either - // the original value or a vector - // of trues - const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (n_components, true) : - component_mask_); - Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, - VectorTools::ExcNoComponentSelected()); - - // a field to store the indices - std::vector face_dofs; - face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); - - typename DH::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) + tables_by_block.resize(1); + tables_by_block[0].reinit(nb, nb); + tables_by_block[0].fill(none); + + for (unsigned int i=0;i &fe = cell->get_fe(); + const unsigned int ib = fe.component_to_block_index(i); + for (unsigned int j=0;j + void + convert_couplings_to_blocks ( + const hp::DoFHandler& dof_handler, + const Table<2, Coupling>& table, + std::vector >& tables_by_block) + { + const hp::FECollection& fe_collection = dof_handler.get_fe(); + tables_by_block.resize(fe_collection.size()); - typename DH::face_iterator face = cell->face(face_no); - if (face->boundary_indicator () == 0) - // face is of the right component + for (unsigned int f=0;f& fe = fe_collection[f]; + + const unsigned int nb = fe.n_blocks(); + tables_by_block[f].reinit(nb, nb); + tables_by_block[f].fill(none); + for (unsigned int i=0;iget_dof_indices (face_dofs, cell->active_fe_index()); - - // enter those dofs into the list - // that match the component - // signature. - for (unsigned int i=0; i &nonzero_component_array - = cell->get_fe().get_nonzero_components (i); - bool nonzero = false; - for (unsigned int c=0; c class DH> + void + make_zero_boundary_constraints (const DH &dof, + ConstraintMatrix &zero_boundary_constraints, + const std::vector &component_mask_) + { + Assert ((component_mask_.size() == 0) || + (component_mask_.size() == dof.get_fe().n_components()), + ExcMessage ("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); + + const unsigned int n_components = DoFTools::n_components(dof); + + // set the component mask to either + // the original value or a vector + // of trues + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (n_components, true) : + component_mask_); + Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, + VectorTools::ExcNoComponentSelected()); + + // a field to store the indices + std::vector face_dofs; + face_dofs.reserve (max_dofs_per_face(dof)); + + typename DH::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const FiniteElement &fe = cell->get_fe(); + + typename DH::face_iterator face = cell->face(face_no); + if (face->boundary_indicator () == 0) + // face is of the right component + { + // get indices and physical + // location on this face + face_dofs.resize (fe.dofs_per_face); + face->get_dof_indices (face_dofs, cell->active_fe_index()); + + // enter those dofs into the list + // that match the component + // signature. + for (unsigned int i=0; i &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + bool nonzero = false; + for (unsigned int c=0; c