RCM: Update lac, part 2: constrained_linear_operator

diff --git a/include/deal.II/lac/constrained_linear_operator.h b/include/deal.II/lac/constrained_linear_operator.h
index 0560f20a3ad33b41a2c75a5a59d316692e67bfb4..96eccd4fc845213e3781907a18412ebd2feee29c 100644 (file)
--- a/include/deal.II/lac/constrained_linear_operator.h
+++ b/include/deal.II/lac/constrained_linear_operator.h
@@ -16,7 +16,7 @@
 #ifndef dealii_constrained_linear_operator_h
 #define dealii_constrained_linear_operator_h
 
-#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/linear_operator.h>
 #include <deal.II/lac/packaged_operation.h>
 
@@ -31,16 +31,16 @@ DEAL_II_NAMESPACE_OPEN
 
 
 /**
- * This function takes a ConstraintMatrix @p constraint_matrix and an operator
- * exemplar @p exemplar (this exemplar is usually a linear operator that
- * describes the system matrix - it is only used to create domain and range
- * vectors of appropriate sizes, its action <tt>vmult</tt> is never used). A
- * LinearOperator object associated with the "homogeneous action" of the
- * underlying ConstraintMatrix object is returned:
+ * This function takes an AffineConstraints object @p constraints and
+ * an operator exemplar @p exemplar (this exemplar is usually a linear
+ * operator that describes the system matrix - it is only used to create
+ * domain and range vectors of appropriate sizes, its action <tt>vmult</tt>
+ * is never used). A LinearOperator object associated with the "homogeneous
+ * action" of the underlying AffineConstraints object is returned:
  *
  * Applying the LinearOperator object on a vector <code>u</code> results in a
  * vector <code>v</code> that stores the result of calling
- * ConstraintMatrix::distribute() on <code>u</code> - with one important
+ * AffineConstraints::distribute() on <code>u</code> - with one important
  * difference: inhomogeneities are not applied, but always treated as 0
  * instead.
  *
@@ -62,12 +62,12 @@ DEAL_II_NAMESPACE_OPEN
 template <typename Range, typename Domain>
 LinearOperator<Range, Domain>
 distribute_constraints_linear_operator(
-  const ConstraintMatrix &             constraint_matrix,
-  const LinearOperator<Range, Domain> &exemplar)
+  const AffineConstraints<typename Range::value_type> &constraints,
+  const LinearOperator<Range, Domain> &                exemplar)
 {
   LinearOperator<Range, Domain> return_op = exemplar;
 
-  return_op.vmult_add = [&constraint_matrix](Range &v, const Domain &u) {
+  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
     Assert(!dealii::PointerComparison::equal(&v, &u),
            dealii::ExcMessage("The domain and range vectors must be different "
                               "storage locations"));
@@ -77,7 +77,7 @@ distribute_constraints_linear_operator(
     v += u;
 
     const auto &locally_owned_elements = v.locally_owned_elements();
-    for (const auto &line : constraint_matrix.get_lines())
+    for (const auto &line : constraints.get_lines())
       {
         const auto i = line.index;
         if (locally_owned_elements.is_element(i))
@@ -95,7 +95,7 @@ distribute_constraints_linear_operator(
     v.compress(VectorOperation::add);
   };
 
-  return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u) {
+  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
     Assert(!dealii::PointerComparison::equal(&v, &u),
            dealii::ExcMessage("The domain and range vectors must be different "
                               "storage locations"));
@@ -105,7 +105,7 @@ distribute_constraints_linear_operator(
     v += u;
 
     const auto &locally_owned_elements = v.locally_owned_elements();
-    for (const auto &line : constraint_matrix.get_lines())
+    for (const auto &line : constraints.get_lines())
       {
         const auto i = line.index;
 
@@ -145,7 +145,7 @@ distribute_constraints_linear_operator(
 
 
 /**
- * Given a ConstraintMatrix @p constraint_matrix and an operator exemplar @p
+ * Given a AffineConstraints @p constraints and an operator exemplar @p
  * exemplar, return a LinearOperator that is the projection to the subspace of
  * constrained degrees of freedom, i.e. all entries of the result vector that
  * correspond to unconstrained degrees of freedom are set to zero.
@@ -158,14 +158,14 @@ distribute_constraints_linear_operator(
 template <typename Range, typename Domain>
 LinearOperator<Range, Domain>
 project_to_constrained_linear_operator(
-  const ConstraintMatrix &             constraint_matrix,
-  const LinearOperator<Range, Domain> &exemplar)
+  const AffineConstraints<typename Range::value_type> &constraints,
+  const LinearOperator<Range, Domain> &                exemplar)
 {
   LinearOperator<Range, Domain> return_op = exemplar;
 
-  return_op.vmult_add = [&constraint_matrix](Range &v, const Domain &u) {
+  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
     const auto &locally_owned_elements = v.locally_owned_elements();
-    for (const auto &line : constraint_matrix.get_lines())
+    for (const auto &line : constraints.get_lines())
       {
         const auto i = line.index;
         if (locally_owned_elements.is_element(i))
@@ -177,9 +177,9 @@ project_to_constrained_linear_operator(
     v.compress(VectorOperation::add);
   };
 
-  return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u) {
+  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
     const auto &locally_owned_elements = v.locally_owned_elements();
-    for (const auto &line : constraint_matrix.get_lines())
+    for (const auto &line : constraints.get_lines())
       {
         const auto i = line.index;
         if (locally_owned_elements.is_element(i))
@@ -210,7 +210,7 @@ project_to_constrained_linear_operator(
 
 
 /**
- * Given a ConstraintMatrix object @p constraint_matrix and a LinearOperator
+ * Given a AffineConstraints object @p constraints and a LinearOperator
  * @p linop, this function creates a LinearOperator object consisting of the
  * composition of three operations and a regularization:
  * @code
@@ -218,9 +218,9 @@ project_to_constrained_linear_operator(
  * @endcode
  * with
  * @code
- *   C = distribute_constraints_linear_operator(constraint_matrix, linop);
+ *   C = distribute_constraints_linear_operator(constraints, linop);
  *   Ct = transpose_operator(C);
- *   Id_c = project_to_constrained_linear_operator(constraint_matrix, linop);
+ *   Id_c = project_to_constrained_linear_operator(constraints, linop);
  * @endcode
  * and <code>Id_c</code> is the projection to the subspace consisting of all
  * vector entries associated with constrained degrees of freedom.
@@ -248,20 +248,19 @@ project_to_constrained_linear_operator(
  */
 template <typename Range, typename Domain>
 LinearOperator<Range, Domain>
-constrained_linear_operator(const ConstraintMatrix &constraint_matrix,
-                            const LinearOperator<Range, Domain> &linop)
+constrained_linear_operator(
+  const AffineConstraints<typename Range::value_type> &constraints,
+  const LinearOperator<Range, Domain> &                linop)
 {
-  const auto C =
-    distribute_constraints_linear_operator(constraint_matrix, linop);
-  const auto Ct = transpose_operator(C);
-  const auto Id_c =
-    project_to_constrained_linear_operator(constraint_matrix, linop);
+  const auto C    = distribute_constraints_linear_operator(constraints, linop);
+  const auto Ct   = transpose_operator(C);
+  const auto Id_c = project_to_constrained_linear_operator(constraints, linop);
   return Ct * linop * C + Id_c;
 }
 
 
 /**
- * Given a ConstraintMatrix object @p constraint_matrix, a LinearOperator @p
+ * Given a AffineConstraints object @p constraints, a LinearOperator @p
  * linop and a right-hand side @p right_hand_side, this function creates a
  * PackagedOperation that stores the following computation:
  * @code
@@ -269,7 +268,7 @@ constrained_linear_operator(const ConstraintMatrix &constraint_matrix,
  * @endcode
  * with
  * @code
- *   C = distribute_constraints_linear_operator(constraint_matrix, linop);
+ *   C = distribute_constraints_linear_operator(constraints, linop);
  *   Ct = transpose_operator(C);
  * @endcode
  *
@@ -296,27 +295,26 @@ constrained_linear_operator(const ConstraintMatrix &constraint_matrix,
  */
 template <typename Range, typename Domain>
 PackagedOperation<Range>
-constrained_right_hand_side(const ConstraintMatrix &constraint_matrix,
-                            const LinearOperator<Range, Domain> &linop,
-                            const Range &right_hand_side)
+constrained_right_hand_side(
+  const AffineConstraints<typename Range::value_type> &constraints,
+  const LinearOperator<Range, Domain> &                linop,
+  const Range &                                        right_hand_side)
 {
   PackagedOperation<Range> return_comp;
 
   return_comp.reinit_vector = linop.reinit_range_vector;
 
-  return_comp.apply_add =
-    [&constraint_matrix, &linop, &right_hand_side](Range &v) {
-      const auto C =
-        distribute_constraints_linear_operator(constraint_matrix, linop);
-      const auto Ct = transpose_operator(C);
+  return_comp.apply_add = [&constraints, &linop, &right_hand_side](Range &v) {
+    const auto C  = distribute_constraints_linear_operator(constraints, linop);
+    const auto Ct = transpose_operator(C);
 
-      GrowingVectorMemory<Domain>            vector_memory;
-      typename VectorMemory<Domain>::Pointer k(vector_memory);
-      linop.reinit_domain_vector(*k, /*bool fast=*/false);
-      constraint_matrix.distribute(*k);
+    GrowingVectorMemory<Domain>            vector_memory;
+    typename VectorMemory<Domain>::Pointer k(vector_memory);
+    linop.reinit_domain_vector(*k, /*bool fast=*/false);
+    constraints.distribute(*k);
 
-      v += Ct * (right_hand_side - linop * *k);
-    };
+    v += Ct * (right_hand_side - linop * *k);
+  };
 
   // lambda capture expressions are a C++14 feature...
   const auto apply_add = return_comp.apply_add;