// ---------------------------------------------------------------------
//
-// Copyright (C) 2006 - 2015 by the deal.II authors
+// Copyright (C) 2006 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// example
#include <deal.II/lac/block_sparsity_pattern.h>
-#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
DoFRenumbering::Cuthill_McKee(dof);
DoFRenumbering::component_wise(dof);
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
DoFTools::make_hanging_node_constraints(dof, constraints);
constraints.close();
// Dirichlet conditions are to be set for the velocity only. To this end,
// we use a ComponentMask that only selects the velocity components. The
// component mask is obtained from the finite element by specifying the
- // particular components we want. Since we use adaptively refined grids
- // the constraint matrix needs to be first filled with hanging node
+ // particular components we want. Since we use adaptively refined grids,
+ // the affine constraints object needs to be first filled with hanging node
// constraints generated from the DoF handler. Note the order of the two
// functions — we first compute the hanging node constraints, and
- // then insert the boundary values into the constraint matrix. This makes
+ // then insert the boundary values into the constraints object. This makes
// sure that we respect H<sup>1</sup> conformity on boundaries with
// hanging nodes (in three space dimensions), where the hanging node needs
// to dominate the Dirichlet boundary values.
using namespace DoFTools;
make_periodicity_constraints(face_1,
face_2,
- constraint_matrix,
+ affine_constraints,
component_mask = <default value>;
face_orientation = <default value>,
face_flip = <default value>,
@p face_orientation, @p face_flip and @p face_orientation. For a closer description
have a look at the documentation of DoFTools::make_periodicity_constraints.
The remaining parameters are the same as for the high level interface apart
-from the self-explaining @p component_mask and @p constraint_matrix.
+from the self-explaining @p component_mask and @p affine_constraints.
<a name="problem"></a>
/**
* This function does essentially the same as the previous
- * make_flux_sparsity_pattern() function but allows the application of
- * a constraint matrix. This is useful in the case where some components
- * of a finite element are continuous and some discontinuous, allowing
- * constraints to be imposed on the continuous part while also building
- * the flux terms needed for the discontinuous part.
+ * make_flux_sparsity_pattern() function but allows the application of an
+ * AffineConstraints object. This is useful in the case where some
+ * components of a finite element are continuous and some discontinuous,
+ * allowing constraints to be imposed on the continuous part while also
+ * building the flux terms needed for the discontinuous part.
*/
template <typename DoFHandlerType,
typename SparsityPatternType,
* @ref step_6 "step-6"
* tutorial program and is used in almost all following programs as well.
*
- * This function does not clear the constraint matrix object before use, in
+ * This function does not clear the AffineConstraints 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 AffineConstraints::clear() function.
/**
- * Make a constraint matrix for the constraints that result from zero
- * boundary values on the given boundary indicator.
+ * Add constraints to @p zero_boundary_constraints corresponding to
+ * enforcing a zero boundary condition on the given boundary indicator.
*
* This function constrains all degrees of freedom on the given part of the
* boundary.
*
* <h4>Interpolation matrices in two dimensions</h4>
*
- * In addition to the fields discussed above for 1D, a constraint matrix is
- * needed to describe hanging node constraints if the finite element has
+ * In addition to the fields discussed above for 1D, a constraint matrix
+ * is needed to describe hanging node constraints if the finite element has
* degrees of freedom located on edges or vertices. These constraints are
* represented by an $m\times n$-matrix #interface_constraints, where <i>m</i>
* is the number of degrees of freedom on the refined side without the corner
* weights. This is sometimes tricky since the nodes on the edge may have
* different local numbers.
*
- * For the constraint matrix this means the following: if a degree of freedom
- * on one edge of a face is constrained by some other nodes on the same edge
- * with some weights, then the weights have to be exactly the same as those
+ * For the constraint matrix this means the following: if a degree of
+ * freedom on one edge of a face is constrained by some other nodes on the same
+ * edge with some weights, then the weights have to be exactly the same as those
* for constrained nodes on the three other edges with respect to the
* corresponding nodes on these edges. If this isn't the case, you will get
* into trouble with the AffineConstraints class that is the primary consumer
* numbering of all (global) DoFs that are elements in the filter.
*
* If, for example, the filter represents the range <tt>[10,20)</tt>, and
- * the constraint matrix @p constraints_in includes the global indices
+ * the constraints object @p constraints_in includes the global indices
* <tt>{7,13,14}</tt>, the indices <tt>{3,4}</tt> are added to the calling
- * constraint matrix (since 13 and 14 are elements in the filter and element
+ * constraints object (since 13 and 14 are elements in the filter and element
* 13 is the fourth element in the index, and 14 is the fifth).
*
* This function provides an easy way to create a AffineConstraints for
* system.
*
* Since this function adds new nonzero entries to the sparsity pattern, the
- * given sparsity pattern must not be compressed. The constraint matrix
- * (i.e., the current object) must be closed. The sparsity pattern is
- * compressed at the end of the function.
+ * given sparsity pattern must not be compressed. The current object must be
+ * closed. The sparsity pattern is compressed at the end of the function.
*/
void
condense(SparsityPattern &sparsity) const;
* condensed and compressed. This function is the appropriate choice for
* applying inhomogeneous constraints.
*
- * The constraint matrix object must be closed to call this function.
+ * The current object object must be closed to call this function.
*
* See the general documentation of this class for more detailed
* information.
* of freedom.
*
* This function is not typically called from user code, but is used in the
- * DoFTools::make_sparsity_pattern() function when passed a constraint
- * matrix object.
+ * DoFTools::make_sparsity_pattern() function when passed an
+ * AffineConstraints object.
*
* @note This function in itself is thread-safe, i.e., it works properly
* also when several threads call it simultaneously. However, the function
*/
/**
- * This class represents one line of a constraint matrix.
+ * This class represents one constraint in an AffineConstraints object.
*/
struct ConstraintLine
{
*/
DeclException1(ExcRowNotStoredHere,
size_type,
- << "The index set given to this constraint matrix indicates "
+ << "The index set given to this constraints object indicates "
<< "constraints for degree of freedom " << arg1
<< " should not be stored by this object, but a constraint "
<< "is being added.");
DeclException2(ExcColumnNotStoredHere,
size_type,
size_type,
- << "The index set given to this constraint matrix indicates "
+ << "The index set given to this constraints object indicates "
<< "constraints using degree of freedom " << arg2
<< " should not be stored by this object, but a constraint "
<< "for degree of freedom " << arg1 << " uses it.");
return;
Assert(filter.size() > constraints.lines.back().index,
- ExcMessage("Filter needs to be larger than constraint matrix size."));
+ ExcMessage(
+ "The filter must be larger than the given constraints object."));
for (const ConstraintLine &line : constraints.lines)
if (filter.is_element(line.index))
{
// store for each index whether it must be distributed or not. If entry
// is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
+ // otherwise, the number states which line in the AffineConstraints object
// handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
// store for each index whether it must be distributed or not. If entry
// is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
+ // otherwise, the number states which line in the AffineConstraints object
// handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
// store for each index whether it must be distributed or not. If entry
// is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
+ // otherwise, the number states which line in the constraints object
// handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
// store for each index whether it must be distributed or not. If entry
// is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
+ // otherwise, the number states which line in the constraints object
// handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
// store for each index whether it must be distributed or not. If entry
// is invalid_size_type, no distribution is necessary. otherwise, the
- // number states which line in the constraint matrix handles this index
+ // number states which line in the constraints object handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
// store for each index whether it must be distributed or not. If entry
// is numbers::invalid_size_type, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
+ // otherwise, the number states which line in the constraints object
// handles this index
std::vector<size_type> distribute(sparsity.n_rows(),
numbers::invalid_size_type);
* then added to the initial guess. This setup can be realized by using a
* vector of two pointers pointing to the same DoFHandler object and a
* vector of two pointers to the two AffineConstraints objects. If the first
- * constraint matrix is the one including the zero Dirichlet constraints,
- * one would give a std::vector<unsigned int>(1, 0) to the initialize()
- * function, i.e., a vector of length 1 that selects exactly the first
- * constraint matrix with index 0.
+ * AffineConstraints object is the one including the zero Dirichlet
+ * constraints, one would give a std::vector<unsigned int>(1, 0) to the
+ * initialize() function, i.e., a vector of length 1 that selects exactly the
+ * first AffineConstraints object with index 0.
*
* For systems of PDEs where the different blocks of MatrixFree are
* associated with different physical components of the equations, adding
* another block with a different AffineConstraints argument solely for the
* purpose of boundary conditions might lead to cumbersome index
* handling. Instead, one could set up a second MatrixFree instance with the
- * different constraint matrix but the same interpretation of blocks, and
- * use that for interpolating inhomogeneous boundary conditions (see also
+ * different AffineConstraints object but the same interpretation of blocks,
+ * and use that for interpolating inhomogeneous boundary conditions (see also
* the discussion in the results section of the step-37 tutorial program):
*
* @code
*
* Furthermore, this call sets up an AffineConstraints object on each
* level that contains possible periodicity constraints in case those
- * have been added to the underlying triangulation. The constraint matrix
- * can be queried by get_level_constraint_matrix(level). Note that the
+ * have been added to the underlying triangulation. The AffineConstraints
+ * object can be queried by get_level_constraints(level). Note that the
* current implementation of periodicity constraints in this class does
* not support rotation matrices in the periodicity definition, i.e., the
* respective argument in the GridTools::collect_periodic_faces() may not
have_boundary_indices() const;
/**
- * Return the level constraint matrix for a given level, containing
+ * Return the AffineConstraints object for a given level, containing
* periodicity constraints (if enabled on the triangulation).
*/
const AffineConstraints<double> &
+ get_level_constraints(const unsigned int level) const;
+
+ /**
+ * Return the AffineConstraints object for a given level, containing
+ * periodicity constraints (if enabled on the triangulation).
+ *
+ * @deprecated Use get_level_constraints instead, which has a more descriptive name.
+ */
+ DEAL_II_DEPRECATED
+ const AffineConstraints<double> &
get_level_constraint_matrix(const unsigned int level) const;
private:
->face(cell->periodic_neighbor_face_no(f))
->get_mg_dof_indices(l, dofs_1, 0);
cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
- // Store periodicity information in the level constraint
- // matrix Skip DoFs for which we've previously entered
- // periodicity constraints already; this can happen, for
- // example, for a vertex dof at a periodic boundary that we
- // visit from more than one cell
+ // Store periodicity information in the level
+ // AffineConstraints object. Skip DoFs for which we've
+ // previously entered periodicity constraints already; this
+ // can happen, for example, for a vertex dof at a periodic
+ // boundary that we visit from more than one cell
for (unsigned int i = 0; i < dofs_per_face; ++i)
if (level_constraints[l].can_store_line(dofs_2[i]) &&
level_constraints[l].can_store_line(dofs_1[i]) &&
inline const AffineConstraints<double> &
-MGConstrainedDoFs::get_level_constraint_matrix(const unsigned int level) const
+MGConstrainedDoFs::get_level_constraints(const unsigned int level) const
{
AssertIndexRange(level, level_constraints.size());
return level_constraints[level];
+inline const AffineConstraints<double> &
+MGConstrainedDoFs::get_level_constraint_matrix(const unsigned int level) const
+{
+ return get_level_constraints(level);
+}
+
+
+
DEAL_II_NAMESPACE_CLOSE
#endif
* distribute_local_to_global into the global matrix and vector. This way,
* adding several matrices from different sources is more complicated and
* you should make sure that you do not mix different ways of applying
- * constraints. Particular caution is necessary when the given constraint
- * matrix contains inhomogeneous constraints: In that case, the matrix
- * assembled this way must be the only matrix (or you need to assemble the
- * <b>same</b> right hand side for <b>every</b> matrix you generate and add
- * together).
+ * constraints. Particular caution is necessary when the given
+ * AffineConstraints object contains inhomogeneous constraints: In that case,
+ * the matrix assembled this way must be the only matrix (or you need to
+ * assemble the <b>same</b> right hand side for <b>every</b> matrix you
+ * generate and add together).
*
* If you want to use boundary conditions with the matrices generated by the
* functions of this namespace in addition to the ones in a possible
- * constraint matrix, you have to use a function like
+ * AffineConstraints object, you have to use a function like
* <tt>apply_boundary_values</tt> with the matrix, solution, and right hand
* side.
*
* Whether this is a problem you need to worry about or not depends on your
* application. The situation is easily corrected, of course, by applying
* AffineConstraints::distribute() to your solution vector after transfer,
- * using a constraint matrix object computed on the new DoFHandler object (you
+ * using a constraints object computed on the new DoFHandler object (you
* probably need to create this object anyway if you have hanging nodes). This
* is also what is done, for example, in step-15.
*
SparseMatrix<number> mass_matrix(sparsity);
Vector<number> tmp(mass_matrix.n());
- // If the constraint matrix does not conflict with the given boundary
+ // If the constraints object does not conflict with the given boundary
// values (i.e., it either does not contain boundary values or it contains
// the same as boundary_values), we can let it call
// distribute_local_to_global straight away, otherwise we need to first
/**
* Add the constraint $(\vec u-\vec u_\Gamma) \| \vec t$ to the list of
* constraints. In 2d, this is a single constraint, in 3d these are two
- * constraints
+ * constraints.
*
* Here, $\vec u$ is represented by the set of given DoF indices, and
* $\vec t$ by the vector specified as the second argument.
cell->face(face)->get_dof_indices(
face_dof_indices, cell->active_fe_index());
- // Add the computed
- // constraints to the
- // constraint matrix,
- // if the degree of
- // freedom is not
- // already constrained.
+ // Add the computed constraints to the constraints
+ // object, if the degree of freedom is not already
+ // constrained.
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
if (dofs_processed[dof] &&
constraints.can_store_line(face_dof_indices[dof]) &&
cell->face(face)->get_dof_indices(
face_dof_indices, cell->active_fe_index());
- // Add the computed constraints to the constraint
- // matrix, assuming the degree of freedom is not
- // already constrained.
+ // Add the computed constraints to the
+ // AffineConstraints object, assuming the degree
+ // of freedom is not already constrained.
for (unsigned int dof = 0; dof < dofs_per_face;
++dof)
{
}
// so all components of this vector dof are constrained. enter
- // this into the constraint matrix
+ // this into the AffineConstraints object
//
// ignore dofs already constrained
const internal::VectorDoFTuple<dim> &dof_indices =
/**
- * Copy constraints into a constraint matrix object.
+ * Copy constraints into an AffineConstraints object.
*
* This function removes zero constraints and those, which constrain a DoF
* which was already eliminated in one of the previous steps of the hp
* hanging node procedure.
*
- * It also suppresses very small entries in the constraint matrix to avoid
- * making the sparsity pattern fuller than necessary.
+ * It also suppresses very small entries in the AffineConstraints object
+ * to avoid making the sparsity pattern fuller than necessary.
*/
template <typename number1, typename number2>
void
Assert(dofs_on_children.size() <= n_dofs_on_children,
ExcInternalError());
- // for each row in the constraint matrix for this line:
+ // for each row in the AffineConstraints object for this line:
for (unsigned int row = 0; row != dofs_on_children.size();
++row)
{
Assert(dofs_on_children.size() <= n_dofs_on_children,
ExcInternalError());
- // for each row in the constraint matrix for this line:
+ // for each row in the AffineConstraints object for this line:
for (unsigned int row = 0; row != dofs_on_children.size();
++row)
{
cell->face(face)->get_dof_indices(
master_dofs, cell->active_fe_index());
- // Now create constraint matrix for the subfaces and
+ // Now create constraints for the subfaces and
// assemble it. ignore all interfaces with artificial
// cells because we can only get to such interfaces if
// the current cell is a ghost cell
subface_interpolation_matrices
[cell->active_fe_index()][subface_fe_index][c]);
- // Add constraints to global constraint matrix.
+ // Add constraints to global AffineConstraints
+ // object.
filter_constraints(master_dofs,
slave_dofs,
*(subface_interpolation_matrices
const FaceIterator & face_1,
const typename identity<FaceIterator>::type &face_2,
const FullMatrix<double> & transformation,
- AffineConstraints<double> & constraint_matrix,
+ AffineConstraints<double> & affine_constraints,
const ComponentMask & component_mask,
const bool face_orientation,
const bool face_flip,
set_periodicity_constraints(face_1,
face_2->child(c),
child_transformation,
- constraint_matrix,
+ affine_constraints,
component_mask,
face_orientation,
face_flip,
if (dofs_2[i] == dofs_1[j])
if (!(is_identity_constrained && target == i))
{
- constraint_matrix.add_line(dofs_2[i]);
+ affine_constraints.add_line(dofs_2[i]);
constraint_set = true;
}
}
face_flip,
face_rotation)];
- if (constraint_matrix.is_constrained(dofs_2[i]))
+ if (affine_constraints.is_constrained(dofs_2[i]))
{
// if the two aren't already identity constrained
// (whichever way around) or already identical (in
bool enter_constraint = false;
// see if this would add an identity constraint
// cycle
- if (!constraint_matrix.is_constrained(dofs_1[j]))
+ if (!affine_constraints.is_constrained(dofs_1[j]))
{
types::global_dof_index new_dof = dofs_2[i];
while (new_dof != dofs_1[j])
- if (constraint_matrix.is_constrained(new_dof))
+ if (affine_constraints.is_constrained(
+ new_dof))
{
const std::vector<
std::pair<types::global_dof_index,
double>> *constraint_entries =
- constraint_matrix
+ affine_constraints
.get_constraint_entries(new_dof);
if (constraint_entries->size() == 1)
new_dof =
if (enter_constraint)
{
- constraint_matrix.add_line(dofs_1[j]);
- constraint_matrix.add_entry(
+ affine_constraints.add_line(dofs_1[j]);
+ affine_constraints.add_entry(
dofs_1[j],
dofs_2[i],
is_identity_constrained ? 1.0 : -1.0);
// enter the constraint. otherwise there is nothing
// for us still to do
bool enter_constraint = false;
- if (!constraint_matrix.is_constrained(dofs_1[j]))
+ if (!affine_constraints.is_constrained(dofs_1[j]))
{
if (dofs_2[i] != dofs_1[j])
enter_constraint = true;
const std::vector<
std::pair<types::global_dof_index, double>>
*constraint_entries =
- constraint_matrix.get_constraint_entries(
+ affine_constraints.get_constraint_entries(
dofs_1[j]);
if (constraint_entries->size() == 1 &&
(*constraint_entries)[0].first == dofs_2[i])
// this pair of constraints means that
// both dofs have to be constrained to
// 0.
- constraint_matrix.add_line(dofs_2[i]);
+ affine_constraints.add_line(dofs_2[i]);
}
}
else
// constraint cycle
types::global_dof_index new_dof = dofs_1[j];
while (new_dof != dofs_2[i])
- if (constraint_matrix.is_constrained(
+ if (affine_constraints.is_constrained(
new_dof))
{
const std::vector<std::pair<
types::global_dof_index,
double>> *constraint_entries =
- constraint_matrix
+ affine_constraints
.get_constraint_entries(new_dof);
if (constraint_entries->size() == 1)
new_dof =
if (enter_constraint)
{
- constraint_matrix.add_line(dofs_2[i]);
- constraint_matrix.add_entry(
+ affine_constraints.add_line(dofs_2[i]);
+ affine_constraints.add_entry(
dofs_2[i],
dofs_1[j],
is_identity_constrained ? 1.0 : -1.0);
}
}
}
- else if (!constraint_matrix.is_constrained(dofs_2[i]))
+ else if (!affine_constraints.is_constrained(dofs_2[i]))
{
// this is just a regular constraint. enter it piece by
// piece
- constraint_matrix.add_line(dofs_2[i]);
+ affine_constraints.add_line(dofs_2[i]);
for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
{
// Query the correct face_index on face_1 respecting
// And finally constrain the two DoFs respecting
// component_mask:
if (transformation(i, jj) != 0)
- constraint_matrix.add_entry(
+ affine_constraints.add_entry(
dofs_2[i], dofs_1[j], transformation(i, jj));
}
}
make_periodicity_constraints(
const FaceIterator & face_1,
const typename identity<FaceIterator>::type &face_2,
- AffineConstraints<double> & constraint_matrix,
+ AffineConstraints<double> & affine_constraints,
const ComponentMask & component_mask,
const bool face_orientation,
const bool face_flip,
make_periodicity_constraints(face_1->child(i),
face_2->child(j),
- constraint_matrix,
+ affine_constraints,
component_mask,
face_orientation,
face_flip,
set_periodicity_constraints(face_2,
face_1,
transformation,
- constraint_matrix,
+ affine_constraints,
component_mask,
face_orientation,
face_flip,
set_periodicity_constraints(face_2,
face_1,
inverse,
- constraint_matrix,
+ affine_constraints,
component_mask,
face_orientation,
face_flip,
set_periodicity_constraints(face_1,
face_2,
transformation,
- constraint_matrix,
+ affine_constraints,
component_mask,
face_orientation,
face_orientation ?
constraint_point(k) = 1. * coord_int / interval;
// The following lines of code should eliminate the problems with
- // the Constraint-Matrix, which appeared for P>=4. The
+ // the constraints object which appeared for P>=4. The
// AffineConstraints class complained about different constraints
- // for the same entry of the Constraint-Matrix. Actually this
+ // for the same entry: Actually, this
// difference could be attributed to FP errors, as it was in the
// range of 1.0e-16. These errors originate in the loss of
// symmetry in the FP approximation of the shape-functions.
std::vector<types::global_dof_index> &dof_indices)
{
if (mg_constrained_dofs != nullptr &&
- mg_constrained_dofs->get_level_constraint_matrix(level)
- .n_constraints() > 0)
+ mg_constrained_dofs->get_level_constraints(level).n_constraints() >
+ 0)
for (auto &ind : dof_indices)
- if (mg_constrained_dofs->get_level_constraint_matrix(level)
+ if (mg_constrained_dofs->get_level_constraints(level)
.is_identity_constrained(ind))
{
- Assert(mg_constrained_dofs->get_level_constraint_matrix(level)
+ Assert(mg_constrained_dofs->get_level_constraints(level)
.get_constraint_entries(ind)
->size() == 1,
ExcInternalError());
- ind = mg_constrained_dofs->get_level_constraint_matrix(level)
+ ind = mg_constrained_dofs->get_level_constraints(level)
.get_constraint_entries(ind)
->front()
.first;
std::vector<types::global_dof_index> &dof_indices)
{
if (mg_constrained_dofs != nullptr &&
- mg_constrained_dofs->get_level_constraint_matrix(level)
- .n_constraints() > 0)
+ mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
for (auto &ind : dof_indices)
- if (mg_constrained_dofs->get_level_constraint_matrix(level)
+ if (mg_constrained_dofs->get_level_constraints(level)
.is_identity_constrained(ind))
{
- Assert(mg_constrained_dofs->get_level_constraint_matrix(level)
+ Assert(mg_constrained_dofs->get_level_constraints(level)
.get_constraint_entries(ind)
->size() == 1,
ExcInternalError());
- ind = mg_constrained_dofs->get_level_constraint_matrix(level)
+ ind = mg_constrained_dofs->get_level_constraints(level)
.get_constraint_entries(ind)
->front()
.first;