Reference documentation for deal.II version 8.1.0
Classes | Enumerations | Functions
Constraints on degrees of freedom
Collaboration diagram for Constraints on degrees of freedom:

Classes

class  ConstraintMatrix
 
struct  ConstraintMatrix::ConstraintLine
 

Enumerations

enum  ConstraintMatrix::MergeConflictBehavior { ConstraintMatrix::no_conflicts_allowed, ConstraintMatrix::left_object_wins, ConstraintMatrix::right_object_wins }
 

Functions

template<int dim, int spacedim, template< int, int > class DH>
void DoFTools::make_zero_boundary_constraints (const DH< dim, spacedim > &dof, const types::boundary_id boundary_indicator, ConstraintMatrix &zero_boundary_constraints, const ComponentMask &component_mask=ComponentMask())
 
template<int dim, int spacedim, template< int, int > class DH>
void DoFTools::make_zero_boundary_constraints (const DH< dim, spacedim > &dof, ConstraintMatrix &zero_boundary_constraints, const ComponentMask &component_mask=ComponentMask())
 

Sparsity Pattern Generation

template<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern (const DH &dof, SparsityPattern &sparsity_pattern, const ConstraintMatrix &constraints=ConstraintMatrix(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
 
template<class DH , class SparsityPattern >
void DoFTools::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 subdomain_id=numbers::invalid_subdomain_id)
 
template<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern (const DH &dof_handler, SparsityPattern &sparsity_pattern, const ConstraintMatrix &constraints, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_unsigned_int)
 
template<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern (const DH &dof, const std::vector< std::vector< bool > > &mask, SparsityPattern &sparsity_pattern) DEAL_II_DEPRECATED
 
template<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern (const DH &dof_row, const DH &dof_col, SparsityPattern &sparsity)
 
template<class DH , class SparsityPattern >
void DoFTools::make_boundary_sparsity_pattern (const DH &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPattern &sparsity_pattern)
 
template<class DH , class SparsityPattern >
void DoFTools::make_boundary_sparsity_pattern (const DH &dof, const typename FunctionMap< DH::space_dimension >::type &boundary_indicators, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPattern &sparsity)
 
template<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern (const DH &dof_handler, SparsityPattern &sparsity_pattern)
 
template<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern (const DH &dof, SparsityPattern &sparsity, const Table< 2, Coupling > &int_mask, const Table< 2, Coupling > &flux_mask)
 

Hanging Nodes

template<class DH >
void DoFTools::make_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints)
 
template<int dim, int spacedim>
void DoFTools::extract_hanging_node_dofs (const DoFHandler< dim, spacedim > &dof_handler, std::vector< bool > &selected_dofs)
 

Interpolation and projection

template<class DH >
void VectorTools::interpolate_boundary_values (const Mapping< DH::dimension, DH::space_dimension > &mapping, const DH &dof, const typename FunctionMap< DH::space_dimension >::type &function_map, ConstraintMatrix &constraints, const ComponentMask &component_mask=ComponentMask())
 
template<class DH >
void VectorTools::interpolate_boundary_values (const Mapping< DH::dimension, DH::space_dimension > &mapping, const DH &dof, const types::boundary_id boundary_component, const Function< DH::space_dimension > &boundary_function, ConstraintMatrix &constraints, const ComponentMask &component_mask=ComponentMask())
 
template<class DH >
void VectorTools::interpolate_boundary_values (const DH &dof, const types::boundary_id boundary_component, const Function< DH::space_dimension > &boundary_function, ConstraintMatrix &constraints, const ComponentMask &component_mask=ComponentMask())
 
template<class DH >
void VectorTools::interpolate_boundary_values (const DH &dof, const typename FunctionMap< DH::space_dimension >::type &function_map, ConstraintMatrix &constraints, const ComponentMask &component_mask=ComponentMask())
 
template<int dim, int spacedim>
void VectorTools::project_boundary_values (const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const typename FunctionMap< spacedim >::type &boundary_functions, const Quadrature< dim-1 > &q, ConstraintMatrix &constraints, std::vector< unsigned int > component_mapping=std::vector< unsigned int >())
 
template<int dim, int spacedim>
void VectorTools::project_boundary_values (const DoFHandler< dim, spacedim > &dof, const typename FunctionMap< spacedim >::type &boundary_function, const Quadrature< dim-1 > &q, ConstraintMatrix &constraints, std::vector< unsigned int > component_mapping=std::vector< unsigned int >())
 
template<int dim>
void VectorTools::project_boundary_values_curl_conforming (const DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const Mapping< dim > &mapping=StaticMappingQ1< dim >::mapping)
 
template<int dim>
void VectorTools::project_boundary_values_curl_conforming (const hp::DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const hp::MappingCollection< dim, dim > &mapping_collection=hp::StaticMappingQ1< dim >::mapping_collection)
 
template<int dim>
void VectorTools::project_boundary_values_div_conforming (const DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const Mapping< dim > &mapping=StaticMappingQ1< dim >::mapping)
 
template<int dim>
void VectorTools::project_boundary_values_div_conforming (const hp::DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const hp::MappingCollection< dim, dim > &mapping_collection=hp::StaticMappingQ1< dim >::mapping_collection)
 
template<int dim, template< int, int > class DH, int spacedim>
void VectorTools::compute_no_normal_flux_constraints (const DH< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, ConstraintMatrix &constraints, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim >::mapping)
 

Detailed Description

This module deals with constraints on degrees of freedom. The central class to deal with constraints is the ConstraintMatrix class.

Constraints typically come from several sources, for example:

In all of these examples, constraints on degrees of freedom are linear and possibly inhomogeneous. In other words, the always have the form $x_{i_1} = \sum_{j=2}^M a_{i_j} x_{i_j} + b_i$. The deal.II class that deals with storing and using these constraints is ConstraintMatrix. The naming stems from the fact that the class originally only stored the (sparse) matrix $a_{i_j}$. The class name component "matrix" no longer makes much sense today since the class has learned to also deal with inhomogeneities $b_i$.

Eliminating constraints

When building the global system matrix and the right hand sides, one can build them without taking care of the constraints, i.e. by simply looping over cells and adding the local contributions to the global matrix and right hand side objects. In order to do actual calculations, you have to 'condense' the linear system: eliminate constrained degrees of freedom and distribute the appropriate values to the unconstrained dofs. This changes the sparsity pattern of the sparse matrices used in finite element calculations and is thus a quite expensive operation. The general scheme of things is then that you build your system, you eliminate (condense) away constrained nodes using the ConstraintMatrix::condense() functions, then you solve the remaining system, and finally you compute the values of constrained nodes from the values of the unconstrained ones using the ConstraintMatrix::distribute() function. Note that the ConstraintMatrix::condense() function is applied to matrix and right hand side of the linear system, while the ConstraintMatrix::distribute() function is applied to the solution vector. This is the method used in the first few tutorial programs, see for example step-6.

This scheme of first building a linear system and then eliminating constrained degrees of freedom is inefficient, and a bottleneck if there are many constraints and matrices are full, i.e. especially for 3d and/or higher order or hp finite elements. Furthermore, it is impossible to implement for parallel computations where a process may not have access to elements of the matrix. We therefore offer a second way of building linear systems, using the ConstraintMatrix::add_entries_local_to_global() and ConstraintMatrix::distribute_local_to_global() functions discussed below. The resulting linear systems are equivalent to those one gets after calling the ConstraintMatrix::condense() functions.

Note
Both ways of applying constraints set the value of the matrix diagonals to constrained entries to a positive entry of the same magnitude as the other entries in the matrix. As a consequence, you need to set up your problem such that the weak form describing the main matrix contribution is not negative definite. Otherwise, iterative solvers such as CG will break down or be considerably slower as GMRES.
While these two ways are equivalent, i.e., the solution of linear systems computed via either approach is the same, the linear systems themselves do not necessarily have the same matrix and right hand side vector entries. Specifically, the matrix diagonal and right hand side entries corresponding to constrained degrees of freedom may be different as a result of the way in which we compute them; they are, however, always chosen in such a way that the solution to the linear system is the same.

Condensing matrices and sparsity patterns

As mentioned above, the first way of using constraints is to build linear systems without regards to constraints and then "condensing" them away. Condensation of a matrix is done in four steps:

To do these steps, you have (at least) two possibilities:

The ConstraintMatrix class provides two sets of condense functions: those taking two arguments refer to the first possibility above, those taking only one do their job in-place and refer to the second possibility.

The condensation functions exist for different argument types. The in-place functions (i.e. those following the second way) exist for arguments of type SparsityPattern, SparseMatrix and BlockSparseMatrix. Note that there are no versions for arguments of type PETScWrappers::SparseMatrix() or any of the other PETSc or Trilinos matrix wrapper classes. This is due to the fact that it is relatively hard to get a representation of the sparsity structure of PETSc matrices, and to modify them; this holds in particular, if the matrix is actually distributed across a cluster of computers. If you want to use PETSc/Trilinos matrices, you can either copy an already condensed deal.II matrix, or assemble the PETSc/Trilinos matrix in the already condensed form, see the discussion below.

Condensing vectors

Condensing vectors works exactly as described above for matrices. Note that condensation is an idempotent operation, i.e. doing it more than once on a vector or matrix yields the same result as doing it only once: once an object has been condensed, further condensation operations don't change it any more.

In contrast to the matrix condensation functions, the vector condensation functions exist in variants for PETSc and Trilinos vectors. However, using them is typically expensive, and should be avoided. You should use the same techniques as mentioned above to avoid their use.

Avoiding explicit condensation

Sometimes, one wants to avoid explicit condensation of a linear system after it has been built at all. There are two main reasons for wanting to do so:

In this case, one possibility is to distribute local entries to the final destinations right at the moment of transferring them into the global matrices and vectors, and similarly build a sparsity pattern in the condensed form at the time it is set up originally.

The ConstraintMatrix class offers support for these operations as well. For example, the ConstraintMatrix::add_entries_local_to_global() function adds nonzero entries to a sparsity pattern object. It not only adds a given entry, but also all entries that we will have to write to if the current entry corresponds to a constrained degree of freedom that will later be eliminated. Similarly, one can use the ConstraintMatrix::distribute_local_to_global() functions to directly distribute entries in vectors and matrices when copying local contributions into a global matrix or vector. These calls make a subsequent call to ConstraintMatrix::condense() unnecessary. For examples of their use see the tutorial programs referenced above.

Note that, despite their name which describes what the function really does, the ConstraintMatrix::distribute_local_to_global() functions has to be applied to matrices and right hand side vectors, whereas the ConstraintMatrix::distribute() function discussed below is applied to the solution vector after solving the linear system.

Distributing constraints

After solving the condensed system of equations, the solution vector has to be "distributed": the modification to the original linear system that results from calling ConstraintMatrix::condense leads to a linear system that solves correctly for all degrees of freedom that are unconstrained but leaves the values of constrained degrees of freedom undefined. To get the correct values also for these degrees of freedom, you need to "distribute" the unconstrained values also to their constrained colleagues. This is done by the two ConstraintMatrix::distribute() functions, one working with two vectors, one working in-place. The operation of distribution undoes the condensation process in some sense, but it should be noted that it is not the inverse operation. Basically, distribution sets the values of the constrained nodes to the value that is computed from the constraint given the values of the unconstrained nodes plus possible inhomogeneities.

Treatment of inhomogeneous constraints

In case some constraint lines have inhomogeneities (which is typically the case if the constraint comes from implementation of inhomogeneous boundary conditions), the situation is a bit more complicated than if the only constraints were due to hanging nodes alone. This is because the elimination of the non-diagonal values in the matrix generate contributions in the eliminated rows in the vector. This means that inhomogeneities can only be handled with functions that act simultaneously on a matrix and a vector. This means that all inhomogeneities are ignored in case the respective condense function is called without any matrix (or if the matrix has already been condensed before).

The use of ConstraintMatrix for implementing Dirichlet boundary conditions is discussed in the step-22 tutorial program. A further example that applies the ConstraintMatrix is step-41. The situation here is little more complicated, because there we have some constraints which are not at the boundary. There are two ways to apply inhomogeneous constraints after creating the ConstraintMatrix:

First approach:

Second approach:

Of course, both approaches lead to the same final answer but in different ways. Using approach (i.e., when using use_inhomogeneities_for_rhs = false in ConstraintMatrix::distribute_local_to_global()), the linear system we build has zero entries in the right hand side in all those places where a degree of freedom is constrained, and some positive value on the matrix diagonal of these lines. Consequently, the solution vector of the linear system will have a zero value for inhomogeneously constrained degrees of freedom and we need to call ConstraintMatrix::distribute() to give these degrees of freedom their correct nonzero values.

On the other hand, in the second approach, the matrix diagonal element and corresponding right hand side entry for inhomogeneously constrained degrees of freedom are so that the solution of the linear system already has the correct value (e.g., if the constraint is that $x_{13}=42$ then row $13$ if the matrix is empty with the exception of the diagonal entry, and $b_{13}/A_{13,13}=42$ so that the solution of $Ax=b$ must satisfy $x_{13}=42$ as desired). As a consequence, we do not need to call ConstraintMatrix::distribute() after solving to fix up inhomogeneously constrained components of the solution, though there is also no harm in doing so.

There remains the question of which of the approaches to take and why we need to set to zero the values of the solution vector in the first approach. The answer to both questions has to do with how iterative solvers solve the linear system. To this end, consider that we typically stop iterations when the residual has dropped below a certain fraction of the norm of the right hand side, or, alternatively, a certain fraction of the norm of the initial residual. Now consider this:

In addition to these considerations, consider the case where we have inhomogeneous constraints of the kind $x_{3}=\tfrac 12 x_1 + \tfrac 12$, e.g., from a hanging node constraint of the form $x_{3}=\tfrac 12 (x_1 + x_2)$ where $x_2$ is itself constrained by boundary values to $x_2=1$. In this case, the ConstraintMatrix can of course not figure out what the final value of $x_3$ should be and, consequently, can not set the solution vector's third component correctly. Thus, the second approach will not work and you should take the first.

Dealing with conflicting constraints

There are situations where degrees of freedom are constrained in more than one way, and sometimes in conflicting ways. Consider, for example the following situation:

conflicting_constraints.png

Here, degree of freedom $x_0$ marked in blue is a hanging node. If we used trilinear finite elements, i.e. FE_Q(1), then it would carry the constraint $x_0=\frac 12 (x_{1}+x_{2})$. On the other hand, it is at the boundary, and if we have imposed boundary conditions $u|_{\partial\Omega}=g$ then we will have the constraint $x_0=g_0$ where $g_0$ is the value of the boundary function $g(\mathbf x)$ at the location of this degree of freedom.

So, which one will win? Or maybe: which one should win? There is no good answer to this question:

That said, what should you do if you know what you want is this:

Either behavior can also be achieved by building two separate ConstraintMatrix objects and calling ConstraintMatrix::merge function with a particular second argument.

Enumeration Type Documentation

An enum that describes what should happen if the two ConstraintMatrix objects involved in a call to the merge() function happen to have constraints on the same degrees of freedom.

Enumerator
no_conflicts_allowed 

Throw an exception if the two objects concerned have conflicting constraints on the same degree of freedom.

left_object_wins 

In an operation cm1.merge(cm2), if cm1 and cm2 have constraints on the same degree of freedom, take the one from cm1.

right_object_wins 

In an operation cm1.merge(cm2), if cm1 and cm2 have constraints on the same degree of freedom, take the one from cm2.

Definition at line 158 of file constraint_matrix.h.

Function Documentation

template<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern ( const DH &  dof,
SparsityPattern sparsity_pattern,
const ConstraintMatrix constraints = ConstraintMatrix(),
const bool  keep_constrained_dofs = true,
const types::subdomain_id  subdomain_id = numbers::invalid_subdomain_id 
)

Locate non-zero entries 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, 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 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.

template<class DH , class SparsityPattern >
void DoFTools::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  subdomain_id = numbers::invalid_subdomain_id 
)

Locate non-zero entries for vector valued finite elements. This function does mostly the same as the previous make_sparsity_pattern(), but it is specialized for vector finite elements and allows to specify which variables couple in which equation. For example, if wanted to solve the Stokes equations,

\begin{align*} -\Delta \mathbf u + \nabla p &= 0,\\ \text{div}\ u &= 0 \end{align*}

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:

\[ \left[ \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 0 \end{array} \right] \]

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

Todo:
Not implemented for hp::DoFHandler.

As mentioned before, the creation of the sparsity pattern is a purely local process and the sparsity pattern does not provide for entries introduced by the elimination of hanging nodes. They have to be taken care of by a call to ConstraintMatrix::condense() afterwards.

Alternatively, the constraints on degrees of freedom can already be taken into account at the time of creating the sparsity pattern. For this, pass the ConstraintMatrix object as the third argument to the current function. No call to ConstraintMatrix::condense() is then necessary. This process is explained in 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 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.

template<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern ( const DH &  dof_handler,
SparsityPattern sparsity_pattern,
const ConstraintMatrix constraints,
const bool  keep_constrained_dofs = true,
const types::subdomain_id  subdomain_id = numbers::invalid_unsigned_int 
)

This function does the same as the other with the same name, but it gets a ConstraintMatrix additionally. This is for the case where you have fluxes but constraints as well.

template<class DH >
void DoFTools::make_hanging_node_constraints ( const DH &  dof_handler,
ConstraintMatrix constraints 
)

Compute the constraints resulting from the presence of hanging nodes. Hanging nodes are best explained using a small picture:

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 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 DoFHandler or hp::DoFHandler.

template<int dim, int spacedim, template< int, int > class DH>
void DoFTools::make_zero_boundary_constraints ( const DH< dim, spacedim > &  dof,
const types::boundary_id  boundary_indicator,
ConstraintMatrix zero_boundary_constraints,
const ComponentMask component_mask = ComponentMask() 
)

Make a constraint matrix for the constraints that result from zero boundary values on the given boundary indicator.

This function constrains all degrees of freedom on the given part of the boundary.

A variant of this function with different arguments is used in step-36.

Parameters
dofThe DoFHandler to work on.
boundary_indicatorThe indicator of that part of the boundary for which constraints should be computed. If this number equals numbers::invalid_boundary_id then all boundaries of the domain will be treated.
zero_boundary_constraintsThe constraint object into which the constraints will be written. The new constraints due to zero boundary values will simply be added, preserving any other constraints previously present. However, this will only work if the previous content of that object consists of constraints on degrees of freedom that are not located on the boundary treated here. If there are previously existing constraints for degrees of freedom located on the boundary, then this would constitute a conflict. See the Constraints on degrees of freedom module for handling the case where there are conflicting constraints on individual degrees of freedom.
component_maskAn optional component mask that restricts the functionality of this function to a subset of an FESystem. For non-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 by the component mask (see GlossComponentMask). If this argument is omitted, all components of the finite element with degrees of freedom at the boundary will be considered.
See also
Glossary entry on boundary indicators
template<int dim, int spacedim, template< int, int > class DH>
void DoFTools::make_zero_boundary_constraints ( const DH< dim, spacedim > &  dof,
ConstraintMatrix zero_boundary_constraints,
const ComponentMask component_mask = ComponentMask() 
)

Do the same as the previous function, except do it for all parts of the boundary, not just those with a particular boundary indicator. This function is then equivalent to calling the previous one with numbers::invalid_boundary_id as second argument.

This function is used in step-36, for example.

template<class DH >
void VectorTools::interpolate_boundary_values ( const Mapping< DH::dimension, DH::space_dimension > &  mapping,
const DH &  dof,
const typename FunctionMap< DH::space_dimension >::type &  function_map,
ConstraintMatrix constraints,
const ComponentMask component_mask = ComponentMask() 
)

Insert the (algebraic) constraints due to Dirichlet boundary conditions into a ConstraintMatrix constraints. This function identifies the degrees of freedom subject to Dirichlet boundary conditions, adds them to the list of constrained DoFs in constraints and sets the respective inhomogeneity to the value interpolated around the boundary. If this routine encounters a DoF that already is constrained (for instance by a hanging node constraint, see below, or any other type of constraint, e.g. from periodic boundary conditions), the old setting of the constraint (dofs the entry is constrained to, inhomogeneities) is kept and nothing happens.

Note
When combining adaptively refined meshes with hanging node constraints and boundary conditions like from the current function within one ConstraintMatrix object, the hanging node constraints should always be set first, and then the boundary conditions since boundary conditions are not set in the second operation on degrees of freedom that are already constrained. This makes sure that the discretization remains conforming as is needed. See the discussion on conflicting constraints in the module on Constraints on degrees of freedom .

The parameter boundary_component corresponds to the number boundary_indicator of the face.

The flags in the last parameter, component_mask denote which components of the finite element space shall be interpolated. If it is left as specified by the default value (i.e. an empty array), all components are interpolated. If it is different from the default value, it is assumed that the number of entries equals the number of components in the boundary functions and the finite element, and those components in the given boundary function will be used for which the respective flag was set in the component mask. See also GlossComponentMask. As an example, assume that you are solving the Stokes equations in 2d, with variables $(u,v,p)$ and that you only want to interpolate boundary values for the pressure, then the component mask should correspond to (true,true,false).

Note
Whether a component mask has been specified or not, the number of components of the functions in function_map must match that of the finite element used by dof. In other words, for the example above, you need to provide a Function object that has 3 components (the two velocities and the pressure), even though you are only interested in the first two of them. interpolate_boundary_values() will then call this function to obtain a vector of 3 values at each interpolation point but only take the first two and discard the third. In other words, you are free to return whatever you like in the third component of the vector returned by Function::vector_value, but the Function object must state that it has 3 components.

If the finite element used has shape functions that are non-zero in more than one component (in deal.II speak: they are non-primitive), then these components can presently not be used for interpolating boundary values. Thus, the elements in the component mask corresponding to the components of these non-primitive shape functions must be false.

See the general documentation of this class for more information.

Definition at line 1888 of file vector_tools.templates.h.

template<class DH >
void VectorTools::interpolate_boundary_values ( const Mapping< DH::dimension, DH::space_dimension > &  mapping,
const DH &  dof,
const types::boundary_id  boundary_component,
const Function< DH::space_dimension > &  boundary_function,
ConstraintMatrix constraints,
const ComponentMask component_mask = ComponentMask() 
)

Same function as above, but taking only one pair of boundary indicator and corresponding boundary function. The same comments apply as for the previous function, in particular about the use of the component mask and the requires size of the function object.

See also
Glossary entry on boundary indicators

Definition at line 1917 of file vector_tools.templates.h.

template<class DH >
void VectorTools::interpolate_boundary_values ( const DH &  dof,
const types::boundary_id  boundary_component,
const Function< DH::space_dimension > &  boundary_function,
ConstraintMatrix constraints,
const ComponentMask component_mask = ComponentMask() 
)

Calls the other interpolate_boundary_values() function, see above, with mapping=MappingQ1<dim>(). The same comments apply as for the previous function, in particular about the use of the component mask and the requires size of the function object.

See also
Glossary entry on boundary indicators

Definition at line 1935 of file vector_tools.templates.h.

template<class DH >
void VectorTools::interpolate_boundary_values ( const DH &  dof,
const typename FunctionMap< DH::space_dimension >::type &  function_map,
ConstraintMatrix constraints,
const ComponentMask component_mask = ComponentMask() 
)

Calls the other interpolate_boundary_values() function, see above, with mapping=MappingQ1<dim>(). The same comments apply as for the previous function, in particular about the use of the component mask and the requires size of the function object.

Definition at line 1951 of file vector_tools.templates.h.

template<int dim, int spacedim>
void VectorTools::project_boundary_values ( const Mapping< dim, spacedim > &  mapping,
const DoFHandler< dim, spacedim > &  dof,
const typename FunctionMap< spacedim >::type &  boundary_functions,
const Quadrature< dim-1 > &  q,
ConstraintMatrix constraints,
std::vector< unsigned int component_mapping = std::vector<unsigned int>() 
)

Project a function to the boundary of the domain, using the given quadrature formula for the faces. This function identifies the degrees of freedom subject to Dirichlet boundary conditions, adds them to the list of constrained DoFs in constraints and sets the respective inhomogeneity to the value resulting from the projection operation. If this routine encounters a DoF that already is constrained (for instance by a hanging node constraint, see below, or any other type of constraint, e.g. from periodic boundary conditions), the old setting of the constraint (dofs the entry is constrained to, inhomogeneities) is kept and nothing happens.

Note
When combining adaptively refined meshes with hanging node constraints and boundary conditions like from the current function within one ConstraintMatrix object, the hanging node constraints should always be set first, and then the boundary conditions since boundary conditions are not set in the second operation on degrees of freedom that are already constrained. This makes sure that the discretization remains conforming as is needed. See the discussion on conflicting constraints in the module on Constraints on degrees of freedom .

If component_mapping is empty, it is assumed that the number of components of boundary_function matches that of the finite element used by dof.

In 1d, projection equals interpolation. Therefore, interpolate_boundary_values is called.

  • component_mapping: if the components in boundary_functions and dof do not coincide, this vector allows them to be remapped. If the vector is not empty, it has to have one entry for each component in dof. This entry is the component number in boundary_functions that should be used for this component in dof. By default, no remapping is applied.

Definition at line 2215 of file vector_tools.templates.h.

template<int dim, int spacedim>
void VectorTools::project_boundary_values ( const DoFHandler< dim, spacedim > &  dof,
const typename FunctionMap< spacedim >::type &  boundary_function,
const Quadrature< dim-1 > &  q,
ConstraintMatrix constraints,
std::vector< unsigned int component_mapping = std::vector<unsigned int>() 
)

Calls the project_boundary_values() function, see above, with mapping=MappingQ1<dim>().

Definition at line 2242 of file vector_tools.templates.h.

template<int dim>
void VectorTools::project_boundary_values_curl_conforming ( const DoFHandler< dim > &  dof_handler,
const unsigned int  first_vector_component,
const Function< dim > &  boundary_function,
const types::boundary_id  boundary_component,
ConstraintMatrix constraints,
const Mapping< dim > &  mapping = StaticMappingQ1<dim>::mapping 
)

Compute constraints that correspond to boundary conditions of the form $\vec{n}\times\vec{u}=\vec{n}\times\vec{f}$, i.e. the tangential components of $u$ and $f$ shall coincide.

If the ConstraintMatrix constraints contained values or other constraints before, the new ones are added or the old ones overwritten, if a node of the boundary part to be used was already in the list of constraints. This is handled by using inhomogeneous constraints. Please note that when combining adaptive meshes and this kind of constraints, the Dirichlet conditions should be set first, and then completed by hanging node constraints, in order to make sure that the discretization remains consistent. See the discussion on conflicting constraints in the module on Constraints on degrees of freedom .

This function is explecitly written to use with the FE_Nedelec elements. Thus it throws an exception, if it is called with other finite elements.

The second argument of this function denotes the first vector component in the finite element that corresponds to the vector function that you want to constrain. For example, if we want to solve Maxwell's equations in 3d and the finite element has components $(E_x,E_y,E_z,B_x,B_y,B_z)$ and we want the boundary conditions $\vec{n}\times\vec{B}=\vec{n}\times\vec{f}$, then first_vector_component would be 3. Vectors are implicitly assumed to have exactly dim components that are ordered in the same way as we usually order the coordinate directions, i.e. $x$-, $y$-, and finally $z$-component.

The parameter boundary_component corresponds to the number boundary_indicator of the face. numbers::internal_face_boundary_id is an illegal value, since it is reserved for interior faces.

The last argument is denoted to compute the normal vector $\vec{n}$ at the boundary points.

Computing constraints

To compute the constraints we use projection-based interpolation as proposed in Solin, Segeth and Dolezel (Higher order finite elements, Chapman&Hall, 2004) on every face located at the boundary.

First one projects $\vec{f}$ on the lowest-order edge shape functions. Then the remaining part $(I-P_0)\vec{f}$ of the function is projected on the remaining higher-order edge shape functions. In the last step we project $(I-P_0-P_e)\vec{f}$ on the bubble shape functions defined on the face.

See also
Glossary entry on boundary indicators

Definition at line 3224 of file vector_tools.templates.h.

template<int dim>
void VectorTools::project_boundary_values_curl_conforming ( const hp::DoFHandler< dim > &  dof_handler,
const unsigned int  first_vector_component,
const Function< dim > &  boundary_function,
const types::boundary_id  boundary_component,
ConstraintMatrix constraints,
const hp::MappingCollection< dim, dim > &  mapping_collection = hp::StaticMappingQ1< dim >::mapping_collection 
)

Same as above for the hp-namespace.

See also
Glossary entry on boundary indicators
template<int dim>
void VectorTools::project_boundary_values_div_conforming ( const DoFHandler< dim > &  dof_handler,
const unsigned int  first_vector_component,
const Function< dim > &  boundary_function,
const types::boundary_id  boundary_component,
ConstraintMatrix constraints,
const Mapping< dim > &  mapping = StaticMappingQ1<dim>::mapping 
)

Compute constraints that correspond to boundary conditions of the form $\vec{n}^T\vec{u}=\vec{n}^T\vec{f}$, i.e. the normal components of $u$ and $f$ shall coincide.

If the ConstraintMatrix constraints contained values or other constraints before, the new ones are added or the old ones overwritten, if a node of the boundary part to be used was already in the list of constraints. This is handled by using inhomogeneous constraints. Please note that when combining adaptive meshes and this kind of constraints, the Dirichlet conditions should be set first, and then completed by hanging node constraints, in order to make sure that the discretization remains consistent. See the discussion on conflicting constraints in the module on Constraints on degrees of freedom .

This function is explecitly written to use with the FE_RaviartThomas elements. Thus it throws an exception, if it is called with other finite elements.

The second argument of this function denotes the first vector component in the finite element that corresponds to the vector function that you want to constrain. Vectors are implicitly assumed to have exactly dim components that are ordered in the same way as we usually order the coordinate directions, i.e. $x$-, $y$-, and finally $z$-component.

The parameter boundary_component corresponds to the number boundary_indicator of the face. numbers::internal_face_boundary_id is an illegal value, since it is reserved for interior faces.

The last argument is denoted to compute the normal vector $\vec{n}$ at the boundary points.

Computing constraints

To compute the constraints we use interpolation operator proposed in Brezzi, Fortin (Mixed and Hybrid (Finite Element Methods, Springer, 1991) on every face located at the boundary.

See also
Glossary entry on boundary indicators

Definition at line 3805 of file vector_tools.templates.h.

template<int dim>
void VectorTools::project_boundary_values_div_conforming ( const hp::DoFHandler< dim > &  dof_handler,
const unsigned int  first_vector_component,
const Function< dim > &  boundary_function,
const types::boundary_id  boundary_component,
ConstraintMatrix constraints,
const hp::MappingCollection< dim, dim > &  mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection 
)

Same as above for the hp-namespace.

See also
Glossary entry on boundary indicators

Definition at line 3971 of file vector_tools.templates.h.

template<int dim, template< int, int > class DH, int spacedim>
void VectorTools::compute_no_normal_flux_constraints ( const DH< dim, spacedim > &  dof_handler,
const unsigned int  first_vector_component,
const std::set< types::boundary_id > &  boundary_ids,
ConstraintMatrix constraints,
const Mapping< dim, spacedim > &  mapping = StaticMappingQ1<dim>::mapping 
)

This function computes the constraints that correspond to boundary conditions of the form $\vec n \cdot \vec u=0$, i.e. no normal flux if $\vec u$ is a vector-valued quantity. These conditions have exactly the form handled by the ConstraintMatrix class, so instead of creating a map between boundary degrees of freedom and corresponding value, we here create a list of constraints that are written into a ConstraintMatrix. This object may already have some content, for example from hanging node constraints, that remains untouched. These constraints have to be applied to the linear system like any other such constraints, i.e. you have to condense the linear system with the constraints before solving, and you have to distribute the solution vector afterwards.

The use of this function is explained in more detail in step-31. It doesn't make much sense in 1d, so the function throws an exception in that case.

The second argument of this function denotes the first vector component in the finite element that corresponds to the vector function that you want to constrain. For example, if we were solving a Stokes equation in 2d and the finite element had components $(u,v,p)$, then first_vector_component would be zero. On the other hand, if we solved the Maxwell equations in 3d and the finite element has components $(E_x,E_y,E_z,B_x,B_y,B_z)$ and we want the boundary condition $\vec n\cdot \vec B=0$, then first_vector_component would be 3. Vectors are implicitly assumed to have exactly dim components that are ordered in the same way as we usually order the coordinate directions, i.e. $x$-, $y$-, and finally $z$-component. The function assumes, but can't check, that the vector components in the range [first_vector_component,first_vector_component+dim) come from the same base finite element. For example, in the Stokes example above, it would not make sense to use a FESystem<dim>(FE_Q<dim>(2), 1, FE_Q<dim>(1), dim) (note that the first velocity vector component is a $Q_2$ element, whereas all the other ones are $Q_1$ elements) as there would be points on the boundary where the $x$-velocity is defined but no corresponding $y$- or $z$-velocities.

The third argument denotes the set of boundary indicators on which the boundary condition is to be enforced. Note that, as explained below, this is one of the few functions where it makes a difference where we call the function multiple times with only one boundary indicator, or whether we call the function onces with the whole set of boundary indicators at once.

The mapping argument is used to compute the boundary points where the function needs to request the normal vector $\vec n$ from the boundary description.

Note
When combining adaptively refined meshes with hanging node constraints and boundary conditions like from the current function within one ConstraintMatrix object, the hanging node constraints should always be set first, and then the boundary conditions since boundary conditions are not set in the second operation on degrees of freedom that are already constrained. This makes sure that the discretization remains conforming as is needed. See the discussion on conflicting constraints in the module on Constraints on degrees of freedom .

Computing constraints in 2d

Computing these constraints requires some smarts. The main question revolves around the question what the normal vector is. Consider the following situation:

no_normal_flux_1.png

Here, we have two cells that use a bilinear mapping (i.e. MappingQ1). Consequently, for each of the cells, the normal vector is perpendicular to the straight edge. If the two edges at the top and right are meant to approximate a curved boundary (as indicated by the dashed line), then neither of the two computed normal vectors are equal to the exact normal vector (though they approximate it as the mesh is refined further). What is worse, if we constrain $\vec n \cdot \vec u=0$ at the common vertex with the normal vector from both cells, then we constrain the vector $\vec u$ with respect to two linearly independent vectors; consequently, the constraint would be $\vec u=0$ at this point (i.e. all components of the vector), which is not what we wanted.

To deal with this situation, the algorithm works in the following way: at each point where we want to constrain $\vec u$, we first collect all normal vectors that adjacent cells might compute at this point. We then do not constrain $\vec n \cdot \vec u=0$ for each of these normal vectors but only for the average of the normal vectors. In the example above, we therefore record only a single constraint $\vec n \cdot \vec {\bar u}=0$, where $\vec {\bar u}$ is the average of the two indicated normal vectors.

Unfortunately, this is not quite enough. Consider the situation here:

no_normal_flux_2.png

If again the top and right edges approximate a curved boundary, and the left boundary a separate boundary (for example straight) so that the exact boundary has indeed a corner at the top left vertex, then the above construction would not work: here, we indeed want the constraint that $\vec u$ at this point (because the normal velocities with respect to both the left normal as well as the top normal vector should be zero), not that the velocity in the direction of the average normal vector is zero.

Consequently, we use the following heuristic to determine whether all normal vectors computed at one point are to be averaged: if two normal vectors for the same point are computed on different cells, then they are to be averaged. This covers the first example above. If they are computed from the same cell, then the fact that they are different is considered indication that they come from different parts of the boundary that might be joined by a real corner, and must not be averaged.

There is one problem with this scheme. If, for example, the same domain we have considered above, is discretized with the following mesh, then we get into trouble:

no_normal_flux_3.png

Here, the algorithm assumes that the boundary does not have a corner at the point where faces $F1$ and $F2$ join because at that point there are two different normal vectors computed from different cells. If you intend for there to be a corner of the exact boundary at this point, the only way to deal with this is to assign the two parts of the boundary different boundary indicators and call this function twice, once for each boundary indicators; doing so will yield only one normal vector at this point per invocation (because we consider only one boundary part at a time), with the result that the normal vectors will not be averaged. This situation also needs to be taken into account when using this function around reentrant corners on Cartesian meshes. If no-normal-flux boundary conditions are to be enforced on non-Cartesian meshes around reentrant corners, one may even get cycles in the constraints as one will in general constrain different components from the two sides. In that case, set a no-slip constraint on the reentrant vertex first.

Computing constraints in 3d

The situation is more complicated in 3d. Consider the following case where we want to compute the constraints at the marked vertex:

no_normal_flux_4.png

Here, we get four different normal vectors, one from each of the four faces that meet at the vertex. Even though they may form a complete set of vectors, it is not our intent to constrain all components of the vector field at this point. Rather, we would like to still allow tangential flow, where the term "tangential" has to be suitably defined.

In a case like this, the algorithm proceeds as follows: for each cell that has computed two tangential vectors at this point, we compute the unconstrained direction as the outer product of the two tangential vectors (if necessary multiplied by minus one). We then average these tangential vectors. Finally, we compute constraints for the two directions perpendicular to this averaged tangential direction.

There are cases where one cell contributes two tangential directions and another one only one; for example, this would happen if both top and front faces of the left cell belong to the boundary selected whereas only the top face of the right cell belongs to it, maybe indicating the the entire front part of the domain is a smooth manifold whereas the top really forms two separate manifolds that meet in a ridge, and that no-flux boundary conditions are only desired on the front manifold and the right one on top. In cases like these, it's difficult to define what should happen. The current implementation simply ignores the one contribution from the cell that only contributes one normal vector. In the example shown, this is acceptable because the normal vector for the front face of the left cell is the same as the normal vector provided by the front face of the right cell (the surface is planar) but it would be a problem if the front manifold would be curved. Regardless, it is unclear how one would proceed in this case and ignoring the single cell is likely the best one can do.

Results

Because it makes for good pictures, here are two images of vector fields on a circle and on a sphere to which the constraints computed by this function have been applied:

no_normal_flux_5.png
no_normal_flux_6.png

The vectors fields are not physically reasonable but the tangentiality constraint is clearly enforced. The fact that the vector fields are zero at some points on the boundary is an artifact of the way it is created, it is not constrained to be zero at these points.

See also
Glossary entry on boundary indicators

Definition at line 4109 of file vector_tools.templates.h.

template<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern ( const DH &  dof,
const std::vector< std::vector< bool > > &  mask,
SparsityPattern sparsity_pattern 
)
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<class DH , class SparsityPattern >
void DoFTools::make_sparsity_pattern ( const DH &  dof_row,
const DH &  dof_col,
SparsityPattern sparsity 
)

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<class DH , class SparsityPattern >
void DoFTools::make_boundary_sparsity_pattern ( const DH &  dof,
const std::vector< types::global_dof_index > &  dof_to_boundary_mapping,
SparsityPattern sparsity_pattern 
)

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<class DH , class SparsityPattern >
void DoFTools::make_boundary_sparsity_pattern ( const DH &  dof,
const typename FunctionMap< DH::space_dimension >::type &  boundary_indicators,
const std::vector< types::global_dof_index > &  dof_to_boundary_mapping,
SparsityPattern sparsity 
)

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 set of boundary indicators, a 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<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern ( const DH &  dof_handler,
SparsityPattern sparsity_pattern 
)

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 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<class DH , class SparsityPattern >
void DoFTools::make_flux_sparsity_pattern ( const DH &  dof,
SparsityPattern sparsity,
const Table< 2, Coupling > &  int_mask,
const Table< 2, Coupling > &  flux_mask 
)

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.

Todo:
Not implemented for hp::DoFHandler.
template<int dim, int spacedim>
void DoFTools::extract_hanging_node_dofs ( const DoFHandler< dim, spacedim > &  dof_handler,
std::vector< bool > &  selected_dofs 
)

Select all dofs that will be constrained by interface constraints, i.e. all hanging nodes.

The size of selected_dofs shall equal dof_handler.n_dofs(). Previous contents of this array or overwritten.