constrained_linear_operator: use ConstraintMatrix iterators

diff --git a/include/deal.II/lac/constrained_linear_operator.h b/include/deal.II/lac/constrained_linear_operator.h
index a73df3e7508f761ea4fc48b6dd5b1fcf58d182b7..e72cd5b48fa9d6ae88d297fcf9e5b651b6c89912 100644 (file)
--- a/include/deal.II/lac/constrained_linear_operator.h
+++ b/include/deal.II/lac/constrained_linear_operator.h
@@ -72,20 +72,27 @@ LinearOperator<Range, Domain> distribute_constraints_linear_operator(
            dealii::ExcMessage("The domain and range vectors must be different "
                               "storage locations"));
 
-    for (auto i : v.locally_owned_elements())
+    // First, add vector u to v unconditionally and clean up constrained
+    // degrees of freedom later.
+    v += u;
+
+    const auto &locally_owned_elements = v.locally_owned_elements();
+    for (const auto &line : constraint_matrix.get_lines())
       {
-        if (constraint_matrix.is_constrained(i))
+        const auto i = line.index;
+        if (locally_owned_elements.is_element(i))
           {
-            const auto *entries = constraint_matrix.get_constraint_entries(i);
-            for (types::global_dof_index j = 0; j < entries->size(); ++j)
+            v(i) -= u(i);
+            const auto &entries = line.entries;
+            for (types::global_dof_index j = 0; j < entries.size(); ++j)
               {
-                const auto pos = (*entries)[j].first;
-                v(i) += u(pos) * (*entries)[j].second;
+                const auto pos = entries[j].first;
+                v(i) += u(pos) * entries[j].second;
               }
           }
-        else
-          v(i) += u(i);
       }
+
+    v.compress(VectorOperation::add);
   };
 
   return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u)
@@ -94,20 +101,30 @@ LinearOperator<Range, Domain> distribute_constraints_linear_operator(
            dealii::ExcMessage("The domain and range vectors must be different "
                               "storage locations"));
 
-    for (auto i : v.locally_owned_elements())
+    // First, add vector u to v unconditionally and clean up constrained
+    // degrees of freedom later.
+    v += u;
+
+    const auto &locally_owned_elements = v.locally_owned_elements();
+    for (const auto &line : constraint_matrix.get_lines())
       {
-        if (constraint_matrix.is_constrained(i))
+        const auto i = line.index;
+
+        if (locally_owned_elements.is_element(i))
           {
-            const auto *entries = constraint_matrix.get_constraint_entries(i);
-            for (types::global_dof_index j = 0; j < entries->size(); ++j)
-              {
-                const auto pos = (*entries)[j].first;
-                v(pos) += u(i) * (*entries)[j].second;
-              }
+            v(i) -= u(i);
+          }
+
+        const auto &entries = line.entries;
+        for (types::global_dof_index j = 0; j < entries.size(); ++j)
+          {
+            const auto pos = entries[j].first;
+            if (locally_owned_elements.is_element(pos))
+              v(pos) += u(i) * entries[j].second;
           }
-        else
-          v(i)+=u(i);
       }
+
+    v.compress(VectorOperation::add);
   };
 
   // lambda capture expressions are a C++14 feature...
@@ -150,39 +167,48 @@ LinearOperator<Range, Domain> project_to_constrained_linear_operator(
 
   return_op.vmult_add = [&constraint_matrix](Range &v, const Domain &u)
   {
-    for (auto i : v.locally_owned_elements())
-      if (constraint_matrix.is_constrained(i))
-        v(i) += u(i);
+    const auto &locally_owned_elements = v.locally_owned_elements();
+    for (const auto &line : constraint_matrix.get_lines())
+      {
+        const auto i = line.index;
+        if (locally_owned_elements.is_element(i))
+          {
+            v(i) += u(i);
+          }
+      }
+
+    v.compress(VectorOperation::add);
   };
 
   return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u)
   {
-    for (auto i : v.locally_owned_elements())
-      if (constraint_matrix.is_constrained(i))
-        v(i) += u(i);
+    const auto &locally_owned_elements = v.locally_owned_elements();
+    for (const auto &line : constraint_matrix.get_lines())
+      {
+        const auto i = line.index;
+        if (locally_owned_elements.is_element(i))
+          {
+            v(i) += u(i);
+          }
+      }
+
+    v.compress(VectorOperation::add);
   };
 
-  return_op.vmult = [&constraint_matrix](Range &v, const Domain &u)
+  // lambda capture expressions are a C++14 feature...
+  const auto vmult_add = return_op.vmult_add;
+  return_op.vmult = [vmult_add](Range &v, const Domain &u)
   {
-    Assert(!dealii::PointerComparison::equal(&v, &u),
-           dealii::ExcMessage("The domain and range vectors must be different "
-                              "storage locations"));
-
     v = 0.;
-    for (auto i : v.locally_owned_elements())
-      if (constraint_matrix.is_constrained(i))
-        v(i) = u(i);
+    vmult_add(v, u);
   };
 
-  return_op.Tvmult = [&constraint_matrix](Domain &v, const Range &u)
+  // lambda capture expressions are a C++14 feature...
+  const auto Tvmult_add = return_op.Tvmult_add;
+  return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u)
   {
-    Assert(!dealii::PointerComparison::equal(&v, &u),
-           dealii::ExcMessage("The domain and range vectors must be different "
-                              "storage locations"));
     v = 0.;
-    for (auto i : v.locally_owned_elements())
-      if (constraint_matrix.is_constrained(i))
-        v(i) = u(i);
+    Tvmult_add(v, u);
   };
 
   return return_op;