Make IndexSet compatible with Teuchos::RCP<Tpetra::Map>>

diff --git a/doc/news/changes/minor/20231017SebastianKinnewig b/doc/news/changes/minor/20231017SebastianKinnewig
new file mode 100644 (file)
index 0000000..e09a098
--- /dev/null
+++ b/doc/news/changes/minor/20231017SebastianKinnewig
@@ -0,0 +1,6 @@
+New: New constructor for IndexSet that takes a 
+Teuchos::RCP<Tpetra::Map> as input. And the 
+function IndexSet::make_tpetra_map_rcp() 
+returning a Teuchos::RCP<Tpetra::Map>.
+<br>
+(Sebastian Kinnewig, 2023/10/17)

diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h
index 3da4d349e42495f0337c0487ac92f1827606a864..28868eac534d505578aa143b38dbf5cb84fe0496 100644 (file)
--- a/include/deal.II/base/index_set.h
+++ b/include/deal.II/base/index_set.h
@@ -131,11 +131,20 @@ public:
   operator=(IndexSet &&is) noexcept;
 
 #ifdef DEAL_II_WITH_TRILINOS
+
+#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
+  /**
+   * Constructor from a Trilinos Teuchos::RCP<Tpetra::Map>.
+   */
+  explicit IndexSet(
+    Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>> map);
+#  endif // DEAL_II_TRILINOS_WITH_TPETRA
+
   /**
    * Constructor from a Trilinos Epetra_BlockMap.
    */
   explicit IndexSet(const Epetra_BlockMap &map);
-#endif
+#endif // DEAL_II_WITH_TRILINOS
 
   /**
    * Remove all indices from this index set. The index set retains its size,
@@ -596,6 +605,10 @@ public:
   Tpetra::Map<int, types::signed_global_dof_index>
   make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
                   const bool     overlapping  = false) const;
+
+  Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
+  make_tpetra_map_rcp(const MPI_Comm communicator = MPI_COMM_WORLD,
+                      const bool     overlapping  = false) const;
 #  endif
 #endif
 

diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
index ee5d051143e4535f6ea56a5f2b3208395b6ae4bf..1d5c78161bea1271def927e649132feb0598118f 100644 (file)
--- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
@@ -72,8 +72,7 @@ namespace LinearAlgebra
                            const MPI_Comm  communicator)
       : Subscriptor()
       , vector(new Tpetra::Vector<Number, int, types::signed_global_dof_index>(
-          Teuchos::rcp(new Tpetra::Map<int, types::signed_global_dof_index>(
-            parallel_partitioner.make_tpetra_map(communicator, false)))))
+          parallel_partitioner.make_tpetra_map_rcp(communicator, false)))
     {}
 
 
@@ -84,13 +83,12 @@ namespace LinearAlgebra
                            const MPI_Comm  communicator,
                            const bool      omit_zeroing_entries)
     {
-      Tpetra::Map<int, types::signed_global_dof_index> input_map =
-        parallel_partitioner.make_tpetra_map(communicator, false);
-      if (vector->getMap()->isSameAs(input_map) == false)
+      Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>> input_map =
+        parallel_partitioner.make_tpetra_map_rcp(communicator, false);
+      if (vector->getMap()->isSameAs(*input_map) == false)
         vector = std::make_unique<
           Tpetra::Vector<Number, int, types::signed_global_dof_index>>(
-          Teuchos::rcp(
-            new Tpetra::Map<int, types::signed_global_dof_index>(input_map)));
+          input_map);
       else if (omit_zeroing_entries == false)
         {
           vector->putScalar(0.);

diff --git a/source/base/index_set.cc b/source/base/index_set.cc
index 64c8b524d8868c93096acdff538c5a1f07013ff0..17108ea82a2f55cefa13884e36d9021d745090b5 100644 (file)
--- a/source/base/index_set.cc
+++ b/source/base/index_set.cc
@@ -27,6 +27,9 @@
 #  endif
 #  include <Epetra_Map.h>
 #  include <Epetra_SerialComm.h>
+#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#    include <Tpetra_Map.hpp>
+#  endif
 #endif
 
 DEAL_II_NAMESPACE_OPEN
@@ -35,6 +38,36 @@ DEAL_II_NAMESPACE_OPEN
 
 #ifdef DEAL_II_WITH_TRILINOS
 
+#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+IndexSet::IndexSet(
+  Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>> map)
+  : is_compressed(true)
+  , index_space_size(1 + map->getMaxAllGlobalIndex())
+  , largest_range(numbers::invalid_unsigned_int)
+{
+  Assert(map->getMinAllGlobalIndex() == 0,
+         ExcMessage(
+           "The Tpetra::Map does not contain the global index 0, "
+           "which means some entries are not present on any processor."));
+
+  // For a contiguous map, we do not need to go through the whole data...
+  if (map->isContiguous())
+    add_range(size_type(map->getMinGlobalIndex()),
+              size_type(map->getMaxGlobalIndex() + 1));
+  else
+    {
+      const size_type n_indices = map->getLocalNumElements();
+      const types::signed_global_dof_index *indices =
+        map->getMyGlobalIndices().data();
+      add_indices(indices, indices + n_indices);
+    }
+  compress();
+}
+
+#  endif // DEAL_II_TRILINOS_WITH_TPETRA
+
+
 // the 64-bit path uses a few different names, so put that into a separate
 // implementation
 
@@ -916,6 +949,15 @@ IndexSet::fill_index_vector(std::vector<size_type> &indices) const
 Tpetra::Map<int, types::signed_global_dof_index>
 IndexSet::make_tpetra_map(const MPI_Comm communicator,
                           const bool     overlapping) const
+{
+  return *make_tpetra_map_rcp(communicator, overlapping);
+}
+
+
+
+Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
+IndexSet::make_tpetra_map_rcp(const MPI_Comm communicator,
+                              const bool     overlapping) const
 {
   compress();
   (void)communicator;
@@ -949,16 +991,17 @@ IndexSet::make_tpetra_map(const MPI_Comm communicator,
   const bool linear =
     overlapping ? false : is_ascending_and_one_to_one(communicator);
   if (linear)
-    return Tpetra::Map<int, types::signed_global_dof_index>(
-      size(),
-      n_elements(),
-      0,
+    return Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>(
+      new Tpetra::Map<int, types::signed_global_dof_index>(
+        size(),
+        n_elements(),
+        0,
 #    ifdef DEAL_II_WITH_MPI
-      Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
+        Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
 #    else
-      Teuchos::rcp(new Teuchos::Comm<int>())
+        Teuchos::rcp(new Teuchos::Comm<int>())
 #    endif
-    );
+          ));
   else
     {
       const std::vector<size_type>                indices = get_index_vector();
@@ -966,16 +1009,17 @@ IndexSet::make_tpetra_map(const MPI_Comm communicator,
       std::copy(indices.begin(), indices.end(), int_indices.begin());
       const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
         int_indices);
-      return Tpetra::Map<int, types::signed_global_dof_index>(
-        size(),
-        arr_view,
-        0,
+      return Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>(
+        new Tpetra::Map<int, types::signed_global_dof_index>(
+          size(),
+          arr_view,
+          0,
 #    ifdef DEAL_II_WITH_MPI
-        Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
+          Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
 #    else
-        Teuchos::rcp(new Teuchos::Comm<int>())
+          Teuchos::rcp(new Teuchos::Comm<int>())
 #    endif
-      );
+            ));
     }
 }
 #  endif

diff --git a/source/lac/trilinos_tpetra_communication_pattern.cc b/source/lac/trilinos_tpetra_communication_pattern.cc
index c1ad9b978abe002073ec7f5c453fdd3c833fe96d..b52585fe46def39476a4998529c75778fc81ef8b 100644 (file)
--- a/source/lac/trilinos_tpetra_communication_pattern.cc
+++ b/source/lac/trilinos_tpetra_communication_pattern.cc
@@ -52,11 +52,9 @@ namespace LinearAlgebra
       comm = std::make_shared<const MPI_Comm>(communicator);
 
       auto vector_space_vector_map =
-        Teuchos::rcp(new Tpetra::Map<int, types::signed_global_dof_index>(
-          locally_owned_indices.make_tpetra_map(*comm, false)));
+        locally_owned_indices.make_tpetra_map_rcp(*comm, false);
       auto read_write_vector_map =
-        Teuchos::rcp(new Tpetra::Map<int, types::signed_global_dof_index>(
-          ghost_indices.make_tpetra_map(*comm, true)));
+        ghost_indices.make_tpetra_map_rcp(*comm, true);
 
       // Target map is read_write_vector_map
       // Source map is vector_space_vector_map. This map must have uniquely