![]() |
Reference documentation for deal.II version 8.1.0
|
Classes | |
class | ConstraintMatrix |
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()) |
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_) |
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) |
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) |
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) |
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) |
template<int dim, int spacedim> | |
void | VectorTools::project_boundary_values (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) |
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) |
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) |
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) |
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) |
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) |
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 . 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
. The class name component "matrix" no longer makes much sense today since the class has learned to also deal with inhomogeneities
.
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.
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:
Use two different sparsity patterns and two different matrices: you may eliminate the rows and columns associated with a constrained degree of freedom, and create a totally new sparsity pattern and a new system matrix. This has the advantage that the resulting system of equations is smaller and free from artifacts of the condensation process and is therefore faster in the solution process since no unnecessary multiplications occur (see below). However, there are two major drawbacks: keeping two matrices at the same time can be quite unacceptable if you're short of memory. Secondly, the condensation process is expensive, since all entries of the matrix have to be copied, not only those which are subject to constraints.
This procedure is therefore not advocated and not discussed in the Tutorial programs. deal.II used to have functions that could perform these shrinking operations, but for the reasons outlined above they were inefficient, rarely used and consequently removed in version 8.0.
Use only one sparsity pattern and one matrix: doing it this way, the condense functions add nonzero entries to the sparsity pattern of the large matrix (with constrained nodes in it) where the condensation process of the matrix will create additional nonzero elements. In the condensation process itself, rows and columns subject to constraints are distributed to the rows and columns of unconstrained nodes. The constrained degrees of freedom remain in place, however, unlike in the first possibility described above. In order not to disturb the solution process, these rows and columns are filled with zeros and an appropriate positive value on the main diagonal (we choose an average of the magnitudes of the other diagonal elements, so as to make sure that the new diagonal entry has the same order of magnitude as the other entries; this preserves the scaling properties of the matrix). The corresponding value in the right hand sides is set to zero. This way, the constrained node will always get the value zero upon solution of the equation system and will not couple to other nodes any more.
This method has the advantage that only one matrix and sparsity pattern is needed thus using less memory. Additionally, the condensation process is less expensive, since not all but only constrained values in the matrix have to be copied. On the other hand, the solution process will take a bit longer, since matrix vector multiplications will incur multiplications with zeroes in the lines subject to constraints. Additionally, the vector size is larger than in the first possibility, resulting in more memory consumption for those iterative solution methods using a larger number of auxiliary vectors (e.g. methods using explicit orthogonalization procedures).
Nevertheless, this process is overall more efficient due to its lower memory consumption and is the one discussed in the first few programs of the Tutorial programs , for example in step-6.
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 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.
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:
Condensation is an expensive operation, in particular if there are many constraints and/or if the matrix has many nonzero entries. Both is typically the case for 3d, or high polynomial degree computations, as well as for hp finite element methods, see for example the hp paper. This is the case discussed in the hp tutorial program, step-27, as well as in step-22 and step-31.
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.
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.
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 then row
if the matrix is empty with the exception of the diagonal entry, and
so that the solution of
must satisfy
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 , e.g., from a hanging node constraint of the form
where
is itself constrained by boundary values to
. In this case, the ConstraintMatrix can of course not figure out what the final value of
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.
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:
Here, degree of freedom marked in blue is a hanging node. If we used trilinear finite elements, i.e. FE_Q(1), then it would carry the constraint
. On the other hand, it is at the boundary, and if we have imposed boundary conditions
then we will have the constraint
where
is the value of the boundary function
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.
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.
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 and that you only want to interpolate boundary values for the pressure, then the component mask should correspond to
(true,true,false)
.
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.
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.
Definition at line 1917 of file vector_tools.templates.h.
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.
Definition at line 1935 of file vector_tools.templates.h.
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.
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.
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.
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.
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 , i.e. the tangential components of
and
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 and we want the boundary conditions
, 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. -,
-, and finally
-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 at the boundary points.
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 on the lowest-order edge shape functions. Then the remaining part
of the function is projected on the remaining higher-order edge shape functions. In the last step we project
on the bubble shape functions defined on the face.
Definition at line 3224 of file vector_tools.templates.h.
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.
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 , i.e. the normal components of
and
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. -,
-, and finally
-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 at the boundary points.
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.
Definition at line 3805 of file vector_tools.templates.h.
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.
Definition at line 3971 of file vector_tools.templates.h.
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 , i.e. no normal flux if
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 , 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 and we want the boundary condition
, 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. -,
-, and finally
-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 element, whereas all the other ones are
elements) as there would be points on the boundary where the
-velocity is defined but no corresponding
- or
-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 from the boundary description.
Computing these constraints requires some smarts. The main question revolves around the question what the normal vector is. Consider the following situation:
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 at the common vertex with the normal vector from both cells, then we constrain the vector
with respect to two linearly independent vectors; consequently, the constraint would be
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 , we first collect all normal vectors that adjacent cells might compute at this point. We then do not constrain
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
, where
is the average of the two indicated normal vectors.
Unfortunately, this is not quite enough. Consider the situation here:
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 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:
Here, the algorithm assumes that the boundary does not have a corner at the point where faces and
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.
The situation is more complicated in 3d. Consider the following case where we want to compute the constraints at the marked vertex:
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.
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:
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.
Definition at line 4109 of file vector_tools.templates.h.
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.
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,
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:
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.
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.
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.
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:
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.
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.
dof | The DoFHandler to work on. |
boundary_indicator | The 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_constraints | The 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_mask | An 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. |
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.
void DoFTools::make_sparsity_pattern | ( | const DH & | dof, |
const std::vector< std::vector< bool > > & | mask, | ||
SparsityPattern & | sparsity_pattern | ||
) |
true
value in the mask is translated into a Coupling::always value in the table) and calls the function above. 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.
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.
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.
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.
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.
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.