Utilities::MPI::this_mpi_process(mpi_communicator);
// helper function
- const auto sort_and_make_unique =
- [](std::vector<ConstraintType> &constraints) {
- std::sort(constraints.begin(),
- constraints.end(),
- [](const ConstraintType &l1, const ConstraintType &l2) {
- return l1.index < l2.index;
+ const auto sort_and_make_unique = [&constraints_in, &locally_owned_dofs](
+ std::vector<ConstraintType>
+ &locally_relevant_constraints) {
+ if (locally_relevant_constraints.empty())
+ return;
+
+ for (auto &entry : locally_relevant_constraints)
+ std::sort(entry.entries.begin(),
+ entry.entries.end(),
+ [](const auto &l1, const auto &l2) {
+ return l1.first < l2.first;
});
- constraints.erase(std::unique(constraints.begin(),
- constraints.end(),
- [](const ConstraintType &l1,
- const ConstraintType &l2) {
- return l1.index == l2.index;
- }),
- constraints.end());
+ std::sort(locally_relevant_constraints.begin(),
+ locally_relevant_constraints.end(),
+ [](const ConstraintType &l1, const ConstraintType &l2) {
+ return l1.index < l2.index;
+ });
+
+ auto equal_with_tol = [](const number d, const number e) {
+ if (std::abs(std::real(d - e)) <
+ std::sqrt(
+ std::numeric_limits<
+ typename numbers::NumberTraits<number>::real_type>::epsilon()))
+ return true;
+ else
+ return false;
};
+ auto read_ptr = locally_relevant_constraints.begin();
+ auto write_ptr = locally_relevant_constraints.begin();
+ // go through sorted locally relevant constraints and look out for
+ // duplicates (same constrained index, same entries) or cases that need
+ // to be augmented (same index, different entries)
+ while (++read_ptr != locally_relevant_constraints.end())
+ {
+ if (read_ptr->index != write_ptr->index)
+ {
+ // if we have a different index, use it here
+ if (++write_ptr != read_ptr)
+ *write_ptr = std::move(*read_ptr);
+ }
+ else // equal global dof index
+ {
+ auto &a = *write_ptr;
+ const auto &b = *read_ptr;
+ Assert(a.index == b.index, ExcInternalError());
+ if (!equal_with_tol(a.inhomogeneity, b.inhomogeneity) &&
+ locally_owned_dofs.is_element(b.index))
+ {
+ Assert(equal_with_tol(b.inhomogeneity,
+ constraints_in.get_inhomogeneity(
+ b.index)),
+ ExcInternalError());
+ a.inhomogeneity = b.inhomogeneity;
+ }
+
+ auto &av = a.entries;
+ const auto &bv = b.entries;
+ // check if entries vectors are equal
+ bool vectors_are_equal = (av.size() == bv.size());
+ for (unsigned int i = 0; vectors_are_equal && i < av.size(); ++i)
+ {
+ if (av[i].first != bv[i].first)
+ vectors_are_equal = false;
+ else
+ Assert(equal_with_tol(av[i].second, bv[i].second),
+ ExcInternalError());
+ }
+
+ // merge entries vectors if different, otherwise ignore the
+ // second entry
+ if (!vectors_are_equal)
+ {
+ typename dealii::AffineConstraints<
+ number>::ConstraintLine::Entries cv;
+ cv.reserve(av.size() + bv.size());
+
+ unsigned int i = 0;
+ unsigned int j = 0;
+
+ for (; (i < av.size()) && (j < bv.size());)
+ {
+ if (av[i].first == bv[j].first)
+ {
+ Assert(equal_with_tol(av[i].second, bv[j].second),
+ ExcInternalError());
+
+ cv.push_back(av[i]);
+ ++i;
+ ++j;
+ }
+ else if (av[i].first < bv[j].first)
+ {
+ cv.push_back(av[i]);
+ ++i;
+ }
+ else
+ {
+ cv.push_back(bv[j]);
+ ++j;
+ }
+ }
+
+ for (; i < av.size(); ++i)
+ cv.push_back(av[i]);
+
+ for (; j < bv.size(); ++j)
+ cv.push_back(bv[j]);
+
+ std::swap(av, cv);
+ }
+ }
+ }
+ ++write_ptr;
+ locally_relevant_constraints.erase(write_ptr,
+ locally_relevant_constraints.end());
+ };
+
// 0) collect constrained indices of the current object
IndexSet constrained_indices(locally_owned_dofs.size());
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2018 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// test AffineConstraints<double>::make_consistent_in_parallel for case
+// where constraints need to be combined between different ranks
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv);
+ MPILogInitAll all;
+
+ const unsigned int my_proc = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+ AssertThrow(n_procs == 3, ExcMessage("Only working on 3 ranks"));
+
+ const unsigned int n_indices = 5;
+ IndexSet owned_elements(n_indices);
+ IndexSet relevant_elements(n_indices);
+ AffineConstraints<double> constraints;
+ if (my_proc == 0)
+ {
+ owned_elements.add_index(0);
+ relevant_elements.add_range(0, 3);
+ constraints.reinit(owned_elements, relevant_elements);
+ }
+ else if (my_proc == 1)
+ {
+ owned_elements.add_range(1, 3);
+ relevant_elements.add_range(0, n_indices);
+ constraints.reinit(owned_elements, relevant_elements);
+ constraints.add_line(0);
+ constraints.add_entry(0, 1, 0.5);
+ constraints.add_entry(0, 2, 0.5);
+ constraints.add_line(2);
+ constraints.add_entry(2, 3, 1.5);
+ }
+ else if (my_proc == 2)
+ {
+ owned_elements.add_range(3, n_indices);
+ relevant_elements.add_range(0, n_indices);
+ constraints.reinit(owned_elements, relevant_elements);
+ constraints.add_line(0);
+ constraints.add_entry(0, 1, 0.5);
+ constraints.add_entry(0, 3, 0.5);
+ }
+
+ constraints.make_consistent_in_parallel(owned_elements,
+ relevant_elements,
+ MPI_COMM_WORLD);
+ constraints.close();
+ constraints.print(deallog.get_file_stream());
+
+ return 0;
+}
projects / dealii.git / commitdiff