In this example, the mean value along the boundary allows just such a
representation, with $C$ being a matrix with just one row (i.e. there is only
-one constraint). In the implementation, we will create a AffineConstraints
+one constraint). In the implementation, we will create an AffineConstraints
object, add one constraint (i.e. add another row to the matrix) referring to the
first boundary node $i_0$, and insert the weights with which all the other nodes
contribute, which in this example happens to be just $-1$.
// interpolated into the finite element space in which we have solved
// the dual problem: But, again as in the
// <code>WeightedResidual::output_solution</code> function we first need
- // to create a AffineConstraints object including the hanging node
+ // to create an AffineConstraints object including the hanging node
// constraints, but this time of the dual finite element space.
AffineConstraints<double> dual_hanging_node_constraints;
DoFTools::make_hanging_node_constraints(DualSolver<dim>::dof_handler,
FESystem<dim> fe;
- ConstraintMatrix hanging_node_constraints;
+ AffineConstraints<double> hanging_node_constraints;
PETScWrappers::MPI::SparseMatrix system_matrix;
// this is that we should not try to first assemble the matrix and
// right hand side as if there were no hanging node constraints and
// boundary values, and then eliminate these in a second step
- // (using, for example, ConstraintMatrix::condense()). Rather, we
+ // (using, for example, AffineConstraints::condense()). Rather, we
// should try to eliminate hanging node constraints before handing
// these entries over to PETSc. This is easy: instead of copying
// elements by hand into the global matrix (as we do in step-4), we
- // use the ConstraintMatrix::distribute_local_to_global() functions
+ // use the AffineConstraints::distribute_local_to_global() functions
// to take care of hanging nodes at the same time. We also already
// did this in step-6. The second step, elimination of boundary
// nodes, could also be done this way by putting the boundary values
- // into the same ConstraintMatrix object as hanging nodes (see the
+ // into the same AffineConstraints object as hanging nodes (see the
// way it is done in step-6, for example); however, it is not
// strictly necessary to do this here because eliminating boundary
// values can be done with only the data stored on each process
// show how to do this better in step-40.) On the other hand,
// distributing hanging node constraints is simple on this local
// copy, using the usual function
- // ConstraintMatrix::distributed(). In particular, we can compute
+ // AffineConstraints::distributed(). In particular, we can compute
// the values of <i>all</i> constrained degrees of freedom,
// whether the current process owns them or not:
hanging_node_constraints.distribute(localized_solution);
// program (we compute on a uniformly refined mesh), but the function
// requires the argument anyway, of course. So we have to create a
// constraint object. In its original state, constraint objects are
- // unsorted, and have to be sorted (using the ConstraintMatrix::close
+ // unsorted, and have to be sorted (using the AffineConstraints::close
// function) before they can be used. This is what we do here, and which is
// why we can't simply call the VectorTools::project function with an
- // anonymous temporary object <code>ConstraintMatrix()</code> as the second
- // argument.
+ // anonymous temporary object <code>AffineConstraints<double>()</code> as the
+ // second argument.
//
// The second point worth mentioning is that we only compute the length of
// the present time step in the middle of solving the linear system
make_grid_and_dofs();
{
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
constraints.close();
VectorTools::project(dof_handler,
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
// function requires a hanging node constraints object, but to be used we
// first need to close it:
{
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
constraints.close();
VectorTools::project(dof_handler,
constraints,
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
// matrix here by simply calling two functions in the library.
//
// Note that we do not take the hanging node constraints into account when
- // assembling the matrices (both functions have a ConstraintMatrix argument
+ // assembling the matrices (both functions have an AffineConstraints argument
// that defaults to an empty object). This is because we are going to
// condense the constraints in run() after combining the matrices for the
// current time-step.
Vector<double> system_rhs;
std::map<types::global_dof_index, double> boundary_values;
- ConstraintMatrix hanging_node_constraints;
+ AffineConstraints<double> hanging_node_constraints;
// @sect5{Private member functions}
Triangulation<dim> triangulation;
double global_Omega_diameter;
- const unsigned int stokes_degree;
- FESystem<dim> stokes_fe;
- DoFHandler<dim> stokes_dof_handler;
- ConstraintMatrix stokes_constraints;
+ const unsigned int stokes_degree;
+ FESystem<dim> stokes_fe;
+ DoFHandler<dim> stokes_dof_handler;
+ AffineConstraints<double> stokes_constraints;
std::vector<IndexSet> stokes_partitioning;
TrilinosWrappers::BlockSparseMatrix stokes_matrix;
TrilinosWrappers::MPI::BlockVector stokes_rhs;
- const unsigned int temperature_degree;
- FE_Q<dim> temperature_fe;
- DoFHandler<dim> temperature_dof_handler;
- ConstraintMatrix temperature_constraints;
+ const unsigned int temperature_degree;
+ FE_Q<dim> temperature_fe;
+ DoFHandler<dim> temperature_dof_handler;
+ AffineConstraints<double> temperature_constraints;
TrilinosWrappers::SparseMatrix temperature_mass_matrix;
TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
// The last step in the loop over all cells is to enter the local
// contributions into the global matrix and vector structures to the
// positions specified in <code>local_dof_indices</code>. Again, we
- // let the ConstraintMatrix class do the insertion of the cell matrix
- // elements to the global matrix, which already condenses the hanging
- // node constraints.
+ // let the AffineConstraints class do the insertion of the cell
+ // matrix elements to the global matrix, which already condenses the
+ // hanging node constraints.
cell->get_dof_indices(local_dof_indices);
if (rebuild_stokes_matrix == true)
const MappingQ<dim> mapping;
- const FESystem<dim> stokes_fe;
- DoFHandler<dim> stokes_dof_handler;
- ConstraintMatrix stokes_constraints;
+ const FESystem<dim> stokes_fe;
+ DoFHandler<dim> stokes_dof_handler;
+ AffineConstraints<double> stokes_constraints;
TrilinosWrappers::BlockSparseMatrix stokes_matrix;
TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
TrilinosWrappers::MPI::BlockVector stokes_rhs;
- FE_Q<dim> temperature_fe;
- DoFHandler<dim> temperature_dof_handler;
- ConstraintMatrix temperature_constraints;
+ FE_Q<dim> temperature_fe;
+ DoFHandler<dim> temperature_dof_handler;
+ AffineConstraints<double> temperature_constraints;
TrilinosWrappers::SparseMatrix temperature_mass_matrix;
TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
// we won't even notice that this part is not parallelized by threads.
//
// Regarding the implementation of inhomogeneous Dirichlet boundary
- // conditions: Since we use the temperature ConstraintMatrix, we could apply
- // the boundary conditions directly when building the respective matrix and
- // right hand side. In this case, the boundary conditions are inhomogeneous,
- // which makes this procedure somewhat tricky since we get the matrix from
- // some other function that uses its own integration and assembly
+ // conditions: Since we use the temperature AffineConstraints object, we
+ // could apply the boundary conditions directly when building the respective
+ // matrix and right hand side. In this case, the boundary conditions are
+ // inhomogeneous, which makes this procedure somewhat tricky since we get the
+ // matrix from some other function that uses its own integration and assembly
// loop. However, the correct imposition of boundary conditions needs the
// matrix data we work on plus the right hand side simultaneously, since the
// right hand side is created by Gaussian elimination on the matrix rows. In
// having the matrix data available, we choose to create a dummy matrix
// <code>matrix_for_bc</code> that we only fill with data when we need it
// for imposing boundary conditions. These positions are exactly those where
- // we have an inhomogeneous entry in the ConstraintMatrix. There are only a
- // few such positions (on the boundary DoFs), so it is still much cheaper to
- // use this function than to create the full matrix here. To implement this,
- // we ask the constraint matrix whether the DoF under consideration is
- // inhomogeneously constrained. In that case, we generate the respective
- // matrix column that we need for creating the correct right hand side. Note
- // that this (manually generated) matrix entry needs to be exactly the entry
- // that we would fill the matrix with — otherwise, this will not work.
+ // we have an inhomogeneous entry in the AffineConstraints<double>. There are
+ // only a few such positions (on the boundary DoFs), so it is still much
+ // cheaper to use this function than to create the full matrix here. To
+ // implement this, we ask the constraint matrix whether the DoF under
+ // consideration is inhomogeneously constrained. In that case, we generate the
+ // respective matrix column that we need for creating the correct right hand
+ // side. Note that this (manually generated) matrix entry needs to be exactly
+ // the entry that we would fill the matrix with — otherwise, this will
+ // not work.
template <int dim>
void BoussinesqFlowProblem<dim>::project_temperature_field()
{
The example uses an ad hoc refinement indicator that shows some usefulness in shock-type problems, and
in the downhill flow example included. We refine according to the squared gradient of the density.
Hanging nodes are handled by computing the numerical flux across cells that are of differing
-refinement levels, rather than using the ConstraintMatrix class as in
+refinement levels, rather than using the AffineConstraints class as in
all other tutorial programs so far. In this way, the example combines
the continuous and DG methodologies. It also simplifies the generation
of the Jacobian because we do not have to track constrained degrees of
// to face terms either if we enforce hanging node constraints
// strongly (as in all previous tutorial programs so far whenever we
// used continuous finite elements -- this enforcement is done by the
- // ConstraintMatrix class together with
+ // AffineConstraints class together with
// DoFTools::make_hanging_node_constraints). In the current program,
// however, we opt to enforce continuity weakly at faces between cells
// of different refinement level, for two reasons: (i) because we can,
// and more importantly (ii) because we would have to thread the
// automatic differentiation we use to compute the elements of the
// Newton matrix from the residual through the operations of the
- // ConstraintMatrix class. This would be possible, but is not trivial,
- // and so we choose this alternative approach.
+ // AffineConstraints class. This would be possible, but is not
+ // trivial, and so we choose this alternative approach.
//
// What needs to be decided is which side of an interface between two
// cells of different refinement level we are sitting on.
later time when we apply boundary values. There are two basic ways
of doing this (either using MatrixTools::apply_boundary_values()
<i>after</i> assembling the linear system, or using
-ConstraintMatrix::distribute_local_to_global() <i>during</i> assembly;
+AffineConstraints::distribute_local_to_global() <i>during</i> assembly;
see the @ref constraints "constraints module" for more information),
but both result in the same: a linear system that has a total
number of rows equal to the number of <i>all</i> degrees of freedom,
// Since this assembler allows us to fill several vectors, the interface is
// a little more complicated as above. The pointers to the vectors have to
- // be stored in a AnyData object. While this seems to cause two extra
+ // be stored in an AnyData object. While this seems to cause two extra
// lines of code here, it actually comes handy in more complex
// applications.
MeshWorker::Assembler::ResidualSimple<Vector<double>> assembler;
F_{\mathcal{F}_k}\\ F_{\mathcal{A}_k}
\end{pmatrix},
@f}
-and then let the ConstraintMatrix class eliminate all constrained degrees of
+and then let the AffineConstraints class eliminate all constrained degrees of
freedom, namely $U^k_{\mathcal{A}_k}=B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k}$,
in the same way as if the dofs in $\mathcal{A}_k$ were Dirichlet data. The
result linear system (the second to last one above) is symmetric and positive
void solve();
void output_results(const unsigned int iteration) const;
- Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
- IndexSet active_set;
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+ AffineConstraints<double> constraints;
+ IndexSet active_set;
TrilinosWrappers::SparseMatrix system_matrix;
TrilinosWrappers::SparseMatrix complete_system_matrix;
// In a sense, this is the central function of this program. It updates the
// active set of constrained degrees of freedom as discussed in the
- // introduction and computes a ConstraintMatrix object from it that can then
+ // introduction and computes an AffineConstraints object from it that can then
// be used to eliminate constrained degrees of freedom from the solution of
// the next iteration. At the same time we set the constrained degrees of
// freedom of the solution to the correct value, namely the height of the
//
// If we decide that the DoF should be part of the active set, we
// add its index to the active set, introduce an inhomogeneous
- // equality constraint in the ConstraintMatrix object, and reset the
+ // equality constraint in the AffineConstraints object, and reset the
// solution value to the height of the obstacle. Finally, the
// residual of the non-contact part of the system serves as an
// additional control (the residual equals the remaining,
A difficulty is handling of the constraints from
the Dirichlet conditions, hanging nodes and the inequality condition that
arises from the contact. For this purpose we create three objects of type
-ConstraintMatrix that describe the various constraints and that we will
+AffineConstraints that describe the various constraints and that we will
combine as appropriate in each iteration.
Compared to step-41, the programs has a few new classes:
// @sect3{The <code>PlasticityContactProblem</code> class template}
- // This is the main class of this program and supplies all functions
- // and variables needed to describe
- // the nonlinear contact problem. It is
- // close to step-41 but with some additional
- // features like handling hanging nodes,
- // a Newton method, using Trilinos and p4est
- // for parallel distributed computing.
- // To deal with hanging nodes makes
- // life a bit more complicated since
- // we need another ConstraintMatrix now.
- // We create a Newton method for the
- // active set method for the contact
- // situation and to handle the nonlinear
+ // This is the main class of this program and supplies all functions and
+ // variables needed to describe the nonlinear contact problem. It is close to
+ // step-41 but with some additional features like handling hanging nodes, a
+ // Newton method, using Trilinos and p4est for parallel distributed computing.
+ // To deal with hanging nodes makes life a bit more complicated since we need
+ // another AffineConstraints object now. We create a Newton method for the
+ // active set method for the contact situation and to handle the nonlinear
// operator for the constitutive law.
//
// The general layout of this class is very much like for most other tutorial
// also step-40 and the @ref distributed documentation module) as
// well as a variety of constraints: those imposed by hanging nodes,
// by Dirichlet boundary conditions, and by the active set of
- // contact nodes. Of the three ConstraintMatrix variables defined
+ // contact nodes. Of the three AffineConstraints variables defined
// here, the first only contains hanging node constraints, the
// second also those associated with Dirichlet boundary conditions,
// and the third these plus the contact constraints.
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
- ConstraintMatrix constraints_hanging_nodes;
- ConstraintMatrix constraints_dirichlet_and_hanging_nodes;
- ConstraintMatrix all_constraints;
+ AffineConstraints<double> constraints_hanging_nodes;
+ AffineConstraints<double> constraints_dirichlet_and_hanging_nodes;
+ AffineConstraints<double> all_constraints;
IndexSet active_set;
Vector<float> fraction_of_plastic_q_points_per_cell;
// between faces), we need to evaluate the gap between the
// deformed object and the obstacle. If the active set
// condition is true, then we add a constraint to the
- // ConstraintMatrix object that the next Newton update needs
+ // AffineConstraints object that the next Newton update needs
// to satisfy, set the solution vector's corresponding element
// to the correct value, and add the index to the IndexSet
// object that stores which degree of freedom is part of the
// At the end of this function, we exchange data between processors updating
// those ghost elements in the <code>solution</code> variable that have been
// written by other processors. We then merge the Dirichlet constraints and
- // those from hanging nodes into the ConstraintMatrix object that already
+ // those from hanging nodes into the AffineConstraints object that already
// contains the active set. We finish the function by outputting the total
// number of actively constrained degrees of freedom for which we sum over
// the number of actively constrained degrees of freedom owned by each
// Newton right hand side and Newton matrix. It looks fairly innocent because
// the heavy lifting happens in the call to
// <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
- // particular in ConstraintMatrix::distribute_local_to_global(), using the
+ // particular in AffineConstraints::distribute_local_to_global(), using the
// constraints we have previously computed.
template <int dim>
void PlasticityContactProblem<dim>::assemble_newton_system(
// mostly it is just setup then solve. Among the complications are:
//
// - For the hanging nodes we have to apply
- // the ConstraintMatrix::set_zero function to newton_rhs.
+ // the AffineConstraints::set_zero function to newton_rhs.
// This is necessary if a hanging node with solution value $x_0$
// has one neighbor with value $x_1$ which is in contact with the
// obstacle and one neighbor $x_2$ which is not in contact. Because
// Darcy solution. We also need a helper function that figures out whether
// we do indeed need to recompute the Darcy solution.
//
- // Unlike step-31, this step uses one more ConstraintMatrix object called
+ // Unlike step-31, this step uses one more AffineConstraints object called
// darcy_preconditioner_constraints. This constraint object is used only for
// assembling the matrix for the Darcy preconditioner and includes hanging
// node constraints as well as Dirichlet boundary value constraints for the
const unsigned int degree;
- const unsigned int darcy_degree;
- FESystem<dim> darcy_fe;
- DoFHandler<dim> darcy_dof_handler;
- ConstraintMatrix darcy_constraints;
+ const unsigned int darcy_degree;
+ FESystem<dim> darcy_fe;
+ DoFHandler<dim> darcy_dof_handler;
+ AffineConstraints<double> darcy_constraints;
- ConstraintMatrix darcy_preconditioner_constraints;
+ AffineConstraints<double> darcy_preconditioner_constraints;
TrilinosWrappers::BlockSparseMatrix darcy_matrix;
TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix;
TrilinosWrappers::MPI::BlockVector second_last_computed_darcy_solution;
- const unsigned int saturation_degree;
- FE_Q<dim> saturation_fe;
- DoFHandler<dim> saturation_dof_handler;
- ConstraintMatrix saturation_constraints;
+ const unsigned int saturation_degree;
+ FE_Q<dim> saturation_fe;
+ DoFHandler<dim> saturation_dof_handler;
+ AffineConstraints<double> saturation_constraints;
TrilinosWrappers::SparseMatrix saturation_matrix;
// applying the constraints (i.e. darcy_preconditioner_constraints) that
// takes care of hanging node and zero Dirichlet boundary condition
// constraints. By doing so, we don't have to do that afterwards, and we
- // later don't have to use ConstraintMatrix::condense and
+ // later don't have to use AffineConstraints::condense and
// MatrixTools::apply_boundary_values, both functions that would need to
// modify matrix and vector entries and so are difficult to write for the
// Trilinos classes where we don't immediately have access to individual
// The last step in the loop over all cells is to enter the local
// contributions into the global matrix and vector structures to the
// positions specified in local_dof_indices. Again, we let the
- // ConstraintMatrix class do the insertion of the cell matrix elements to
- // the global matrix, which already condenses the hanging node
+ // AffineConstraints class do the insertion of the cell matrix
+ // elements to the global matrix, which already condenses the hanging node
// constraints.
typename DoFHandler<dim>::active_cell_iterator cell = darcy_dof_handler
.begin_active(),
// matrix for the left hand side of the saturation linear system by basis
// functions phi_i_s and phi_j_s only. Finally, as usual, we enter the local
// contribution into the global matrix by specifying the position in
- // local_dof_indices. This is done by letting the ConstraintMatrix class do
+ // local_dof_indices. This is done by letting the AffineConstraints class do
// the insertion of the cell matrix elements to the global matrix, which
// already condenses the hanging node constraints.
template <int dim>
const unsigned int n_q_points_f;
// Objects that store the converged solution and right-hand side vectors,
- // as well as the tangent matrix. There is a ConstraintMatrix object used
+ // as well as the tangent matrix. There is an AffineConstraints object used
// to keep track of constraints. We make use of a sparsity pattern
// designed for a block system.
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> tangent_matrix;
BlockVector<double> system_rhs;
// indicating the hanging node constraints. We have none in this program
// So we have to create a constraint object. In its original state, constraint
// objects are unsorted, and have to be sorted (using the
- // ConstraintMatrix::close function) before they can be used. Have a look at
+ // AffineConstraints::close function) before they can be used. Have a look at
// step-21 for more information. We only need to enforce the initial condition
// on the dilatation. In order to do this, we make use of a
// ComponentSelectFunction which acts as a mask and sets the J_component of
make_grid();
system_setup();
{
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
constraints.close();
const ComponentSelectFunction<dim> J_mask(J_component, n_components);
what GridGenerator::parallelogram() documents.
The resulting @p matched_pairs can be used in
-DoFTools::make_periodicity_constraints for populating a AffineConstraints
-with periodicity constraints:
+DoFTools::make_periodicity_constraints for populating an AffineConstraints
+object with periodicity constraints:
@code
DoFTools::make_periodicity_constraints(matched_pairs, constraints);
@endcode
// In order to implement periodic boundary conditions only two functions
// have to be modified:
// - <code>StokesProblem<dim>::setup_dofs()</code>:
-// To populate a AffineConstraints object with periodicity constraints
+// To populate an AffineConstraints object with periodicity constraints
// - <code>StokesProblem<dim>::run()</code>:
// To supply a distributed triangulation with periodicity information.
//
@endcode
The call <code>constraints.add_line(t)</code> tells the
-ConstraintMatrix to start a new constraint for degree of freedom
+AffineConstraints to start a new constraint for degree of freedom
<code>t</code> of the form $x_t=\sum_{l=0}^{N-1} c_{tl} x_l +
b_t$. Typically, one would then proceed to set individual coefficients
-$c_{tl}$ to nonzero values (using ConstraintMatrix::add_entry) or set
+$c_{tl}$ to nonzero values (using AffineConstraints::add_entry) or set
$b_t$ to something nonzero (using
-ConstraintMatrix::set_inhomogeneity); doing nothing as above, funny as
+AffineConstraints::set_inhomogeneity); doing nothing as above, funny as
it looks, simply leaves the constraint to be $x_t=0$, which is exactly
what we need in the current context. The call to
FiniteElement::face_system_to_component_index makes sure that we only set
hp::FECollection<dim> fe_collection;
hp::DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
// Once we have the contributions from cell integrals, we copy them
// into the global matrix (taking care of constraints right away,
- // through the ConstraintMatrix::distribute_local_to_global
+ // through the AffineConstraints::distribute_local_to_global
// function). Note that we have not written anything into the
// <code>local_rhs</code> variable, though we still need to pass it
// along since the elimination of nonzero boundary values requires the
<h3>Implementation of constraints</h3>
The usual way to handle constraints in <code>deal.II</code> is to use
-the ConstraintMatrix class that builds a sparse matrix storing
+the AffineConstraints class that builds a sparse matrix storing
information about which degrees of freedom (DoF) are constrained and
how they are constrained. This format uses an unnecessarily large
amount of memory since there are not so many different types of
needed once during matrix and right hand side assembly, it becomes a
problem when we want to use the matrix-free approach since there this
information has to be accessed every time we apply the operator. Thus,
-instead of a ConstraintMatrix, we use a variable that we call
+instead of an AffineConstraints object, we use a variable that we call
<code>constraint_pool</code> that collects the weights of the
different constraints. Then, we only have to store an identifier of
each constraint in the mesh instead of all the weights. Moreover, we
#else
Triangulation<dim> triangulation;
#endif
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
- IndexSet locally_relevant_dofs;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+ AffineConstraints<double> constraints;
+ IndexSet locally_relevant_dofs;
MatrixFree<dim, double> matrix_free_data;
IndexSet locally_relevant_set;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
vector_t solution;
vector_t system_rhs;
// by calling <code>get_boundary_indices ()</code>. The third
// step is to construct constraints on all those degrees of
// freedom: their value should be zero after each application of
- // the level operators. To this end, we construct ConstraintMatrix
+ // the level operators. To this end, we construct AffineConstraints
// objects for each level, and add to each of these constraints
- // for each degree of freedom. Due to the way the ConstraintMatrix
+ // for each degree of freedom. Due to the way the AffineConstraints class
// stores its data, the function to add a constraint on a single
// degree of freedom and force it to be zero is called
- // ConstraintMatrix::add_line(); doing so for several degrees of
+ // AffineConstraints::add_line(); doing so for several degrees of
// freedom at once can be done using
- // ConstraintMatrix::add_lines():
- std::vector<ConstraintMatrix> boundary_constraints(
+ // AffineConstraints::add_lines():
+ std::vector<AffineConstraints<double>> boundary_constraints(
triangulation.n_global_levels());
- ConstraintMatrix empty_constraints;
+ AffineConstraints<double> empty_constraints;
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
{
// The degrees of freedom corresponding to the skeleton strongly enforce
// Dirichlet boundary conditions, just as in a continuous Galerkin finite
// element method. We can enforce the boundary conditions in an analogous
- // manner through the use of ConstraintMatrix constructs. In
+ // manner through the use of AffineConstraints constructs. In
// addition, hanging nodes are handled in the same way as for
// continuous finite elements: For the face elements which
// only define degrees of freedom on the face, this process sets the
// solution on the refined to be the one from the coarse side.
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
// The usage of the ChunkSparseMatrix class is similar to the usual sparse
// matrices: You need a sparsity pattern of type ChunkSparsityPattern and
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
LA::MPI::BlockSparseMatrix system_matrix;
LA::MPI::BlockSparseMatrix preconditioner_matrix;
DoFHandler<dim> dof_handler;
DoFHandler<dim> velocity_dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
- std::vector<ConstraintMatrix> boundary_constraints(
+ std::vector<AffineConstraints<double>> boundary_constraints(
triangulation.n_levels());
- std::vector<ConstraintMatrix> boundary_interface_constraints(
+ std::vector<AffineConstraints<double>> boundary_interface_constraints(
triangulation.n_levels());
for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
{
// and the update. Additionally, the evaluation point is
// for temporarily holding Newton update in line search. A sparse matrix
// for the pressure mass matrix is created for the operator of a block Schur
- // complement preconditioner. We use one ConstraintMatrix for Dirichlet
- // boundary conditions at the initial step and a zero ConstraintMatrix for the
- // Newton is defined by 1/Re which has been discussed in the introduction.
+ // complement preconditioner. We use one AffineConstraints object for
+ // Dirichlet boundary conditions at the initial step and a zero
+ // AffineConstraints object for the Newton is defined by 1/Re which has been
+ // discussed in the introduction.
template <int dim>
class StationaryNavierStokes
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
- ConstraintMatrix zero_constraints;
- ConstraintMatrix nonzero_constraints;
+ AffineConstraints<double> zero_constraints;
+ AffineConstraints<double> nonzero_constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
// is zero so we do not need to set the right hand side function in this
// tutorial. The number of components of the boundary function is dim+1.
// In practice, the boundary values are
- // applied to our solution through ConstraintMatrix which is obtained by using
- // VectorTools::interpolate_boundary_values. The components of boundary value
- // functions are required to be chosen according to the finite element space.
- // Therefore we have to define the boundary value of pressure even though we
- // actually do not need it.
+ // applied to our solution through an AffineConstraints object which is
+ // obtained by using VectorTools::interpolate_boundary_values. The components
+ // of boundary value functions are required to be chosen according to the
+ // finite element space. Therefore we have to define the boundary value of
+ // pressure even though we actually do not need it.
// The following function represents the boundary values:
template <int dim>
cell->get_dof_indices(local_dof_indices);
- const ConstraintMatrix &constraints_used =
+ const AffineConstraints<double> &constraints_used =
initial_step ? nonzero_constraints : zero_constraints;
if (assemble_matrix)
template <int dim>
void StationaryNavierStokes<dim>::solve(const bool initial_step)
{
- const ConstraintMatrix &constraints_used =
+ const AffineConstraints<double> &constraints_used =
initial_step ? nonzero_constraints : zero_constraints;
SolverControl solver_control(system_matrix.m(),
// The setup function differs in two aspects from step-37. The first is that
// we do not need to interpolate any constraints for the discontinuous
- // ansatz space, and simply pass a dummy ConstraintMatrix object into
+ // ansatz space, and simply pass a dummy AffineConstraints object into
// Matrixfree::reinit(). The second change arises because we need to tell
// MatrixFree to also initialize the data structures for faces. We do this
// by setting update flags for the inner and boundary faces,
<< std::endl;
time.restart();
- ConstraintMatrix dummy;
+ AffineConstraints<double> dummy;
dummy.close();
{
about the details.
The steps you need to make this work are essentially like this:
-- You have to create a ConstraintMatrix object, which (as the name
+- You have to create an AffineConstraints object, which (as the name
suggests) will store all constraints on the finite element space. In
the current context, these are the constraints due to our desire to
keep the solution space continuous even in the presence of hanging
the elements of the finite element space.
- You have to use this object when you copy the local contributions to
the matrix and right hand side into the global objects, by using
- ConstraintMatrix::distribute_local_to_global(). Up until
+ AffineConstraints::distribute_local_to_global(). Up until
now, we have done this ourselves, but now with constraints, this
is where the magic happens and we apply the constraints to the
linear system. What this function does is make sure that the
on hanging nodes get their correct (constrained) value so that the
solution you then visualize or evaluate in other ways in in
fact continuous. This is done by calling
- ConstraintMatrix::distribute() immediately after solving.
+ AffineConstraints::distribute() immediately after solving.
These four steps are really all that is necessary -- it's that simple
from a user perspective. The fact that, in the function calls mentioned
$\mathbf x_j$, then imposing the boundary condition $u=g$ on $\partial\Omega$
simply yields the constraint $U_j=g({\mathbf x}_j)$.
-The ConstraintMatrix can handle such constraints as well, which makes it
+The AffineConstraints class can handle such constraints as well, which makes it
convenient to let the same object we use for hanging node constraints
also deal with these Dirichlet boundary conditions.
This way, we don't need to apply the boundary conditions after assembly
(like we did in the earlier steps).
All that is necessary is that we call the variant of
VectorTools::interpolate_boundary_values() that returns its information
-in a ConstraintMatrix object, rather than the `std::map` we have used
+in an AffineConstraints object, rather than the `std::map` we have used
in previous tutorial programs.
// This is the new variable in the main class. We need an object which holds
// a list of constraints to hold the hanging nodes and the boundary
// conditions.
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
// The sparsity pattern and sparse matrix are deliberately declared in the
// opposite of the order used in step-2 through step-5 to demonstrate the
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
- // We may now populate the ConstraintMatrix with the hanging node
+ // We may now populate the AffineConstraints object with the hanging node
// constraints. Since we will call this function in a loop we first clear
// the current set of constraints from the last system and then compute new
// ones:
// whole boundary) and store the resulting constraints in our
// <code>constraints</code> object. Note that we do not to apply the
// boundary conditions after assembly, like we did in earlier steps: instead
- // we put all constraints on our function space in the ConstraintMatrix. We
- // can add constraints to the ConstraintMatrix in either order: if two
- // constraints conflict then the constraint matrix either abort or throw an
- // exception via the Assert macro.
+ // we put all constraints on our function space in the AffineConstraints
+ // object. We can add constraints to the AffineConstraints object in either
+ // order: if two constraints conflict then the constraint matrix either abort
+ // or throw an exception via the Assert macro.
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(),
// Now we first build our compressed sparsity pattern like we did in the
// previous examples. Nevertheless, we do not copy it to the final sparsity
// pattern immediately. Note that we call a variant of
- // make_sparsity_pattern that takes the ConstraintMatrix as the third
+ // make_sparsity_pattern that takes the AffineConstraints object as the third
// argument. We are letting the routine know that we will never write into
// the locations given by <code>constraints</code> by setting the argument
// <code>keep_constrained_dofs</code> to false (in other words, that we will
//
// Second, to copy the local matrix and vector on each cell into the global
// system, we are no longer using a hand-written loop. Instead, we use
-// ConstraintMatrix::distribute_local_to_global() that internally executes
+// AffineConstraints::distribute_local_to_global() that internally executes
// this loop while performing Gaussian elimination on rows and columns
// corresponding to constrained degrees on freedom.
//
// We continue with gradual improvements. The function that solves the linear
// system again uses the SSOR preconditioner, and is again unchanged except
// that we have to incorporate hanging node constraints. As mentioned above,
-// the degrees of freedom from the ConstraintMatrix corresponding to hanging
-// node constraints and boundary values have been removed from the linear
-// system by giving the rows and columns of the matrix a special
+// the degrees of freedom from the AffineConstraints object corresponding to
+// hanging node constraints and boundary values have been removed from the
+// linear system by giving the rows and columns of the matrix a special
// treatment. This way, the values for these degrees of freedom have wrong,
// but well-defined values after solving the linear system. What we then have
// to do is to use the constraints to assign to them the values that they
SparseMatrix<double> stiffness_matrix;
SparseMatrix<double> coupling_matrix;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
Vector<double> solution;
Vector<double> rhs;
*embedded_dh,
quad,
dsp,
- ConstraintMatrix(),
+ AffineConstraints<double>(),
ComponentMask(),
ComponentMask(),
*embedded_mapping);
*embedded_dh,
quad,
coupling_matrix,
- ConstraintMatrix(),
+ AffineConstraints<double>(),
ComponentMask(),
ComponentMask(),
*embedded_mapping);
projects / dealii.git / commitdiff