From: Matthias Maier Date: Tue, 8 Aug 2017 21:51:54 +0000 (-0500) Subject: constrained_linear_operator: use ConstraintMatrix iterators X-Git-Tag: v9.0.0-rc1~1309^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4740%2Fhead;p=dealii.git 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 a73df3e750..e72cd5b48f 100644 --- a/include/deal.II/lac/constrained_linear_operator.h +++ b/include/deal.II/lac/constrained_linear_operator.h @@ -72,20 +72,27 @@ LinearOperator 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 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 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;