]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add consistent strategy to merge constraints in parallel 14905/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Thu, 16 Mar 2023 16:56:39 +0000 (17:56 +0100)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 9 Jul 2024 15:03:17 +0000 (17:03 +0200)
include/deal.II/lac/affine_constraints.templates.h
tests/lac/constraints_make_consistent_in_parallel_03.cc [new file with mode: 0644]
tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output [new file with mode: 0644]

index 98c1d654413b3a62f4b8ae5e6b393c241266a84a..5bc85a5ab74d7dad4201b614a71a76fe9f840269 100644 (file)
@@ -321,23 +321,125 @@ namespace internal
       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());
 
diff --git a/tests/lac/constraints_make_consistent_in_parallel_03.cc b/tests/lac/constraints_make_consistent_in_parallel_03.cc
new file mode 100644 (file)
index 0000000..37fa665
--- /dev/null
@@ -0,0 +1,75 @@
+// ------------------------------------------------------------------------
+//
+// 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;
+}
diff --git a/tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output b/tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output
new file mode 100644 (file)
index 0000000..452ceeb
--- /dev/null
@@ -0,0 +1,14 @@
+
+    0 1:  0.500000
+    0 3:  1.25000
+    2 3:  1.50000
+
+    0 1:  0.500000
+    0 3:  1.25000
+    2 3:  1.50000
+
+
+    0 1:  0.500000
+    0 3:  1.25000
+    2 3:  1.50000
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.