Fix incorrect size in overlapping EPetra_Map to IndexSet

diff --git a/doc/news/changes/minor/20170815Heister b/doc/news/changes/minor/20170815Heister
new file mode 100644 (file)
index 0000000..0d08221
--- /dev/null
+++ b/doc/news/changes/minor/20170815Heister
@@ -0,0 +1,3 @@
+Fixed: Constructing an IndexSet from an Epetra_Map generated the wrong size.
+<br>
+(Timo Heister, 2017/08/15)

diff --git a/source/base/index_set.cc b/source/base/index_set.cc
index 4ed7b59b3a420e62a19dc8a1f5471df906610669..9543531e65946fce6691fb95a5a52dab193a85dc 100644 (file)
--- a/source/base/index_set.cc
+++ b/source/base/index_set.cc
@@ -41,16 +41,21 @@ DEAL_II_NAMESPACE_OPEN
 IndexSet::IndexSet (const Epetra_Map &map)
   :
   is_compressed (true),
-  index_space_size (map.NumGlobalElements64()),
+  index_space_size (1+map.MaxAllGID64()),
   largest_range (numbers::invalid_unsigned_int)
 {
+  Assert(
+    map.MinAllGID64() == 0,
+    ExcMessage("The Epetra_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.LinearMap())
     add_range(size_type(map.MinMyGID64()), size_type(map.MaxMyGID64()+1));
   else
     {
       const size_type n_indices = map.NumMyElements();
-      size_type *indices = (size_type *)map.MyGlobalElements64();
+      size_type *indices = reinterpret_cast<size_type *>(map.MyGlobalElements64());
       add_indices(indices, indices+n_indices);
     }
   compress();
@@ -63,16 +68,21 @@ IndexSet::IndexSet (const Epetra_Map &map)
 IndexSet::IndexSet (const Epetra_Map &map)
   :
   is_compressed (true),
-  index_space_size (map.NumGlobalElements()),
+  index_space_size (1+map.MaxAllGID()),
   largest_range (numbers::invalid_unsigned_int)
 {
+  Assert(
+    map.MinAllGID() == 0,
+    ExcMessage("The Epetra_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.LinearMap())
     add_range(size_type(map.MinMyGID()), size_type(map.MaxMyGID()+1));
   else
     {
       const size_type n_indices = map.NumMyElements();
-      unsigned int *indices = (unsigned int *)map.MyGlobalElements();
+      unsigned int *indices = reinterpret_cast<unsigned int *>(map.MyGlobalElements());
       add_indices(indices, indices+n_indices);
     }
   compress();
@@ -566,7 +576,7 @@ IndexSet::make_trilinos_map (const MPI_Comm &communicator,
     }
 #endif
 
-  // Find out if the IndexSet is ascending and 1:1. This correspinds to a
+  // Find out if the IndexSet is ascending and 1:1. This corresponds to a
   // linear EpetraMap. Overlapping IndexSets are never 1:1.
   const bool linear = overlapping ? false : is_ascending_and_one_to_one(communicator);
 

diff --git a/tests/trilinos/index_set_02.cc b/tests/trilinos/index_set_02.cc
new file mode 100644 (file)
index 0000000..49fd75f
--- /dev/null
+++ b/tests/trilinos/index_set_02.cc
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test conversion from epetra map back to indexset, this was broken for
+// overlapping IndexSets
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+
+void test ()
+{
+  IndexSet set_my(100);
+  IndexSet set_ghost(100);
+
+  unsigned int n_proc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  if (myid == 0)
+    {
+      set_my.add_range(0,50);
+      set_ghost.add_range(0,50);
+      set_ghost.add_range(55,60);
+    }
+  else if (myid == 1)
+    {
+      set_my.add_range(50,100);
+      set_ghost.add_range(45,100);
+    }
+  else
+    {
+    }
+
+  auto check = [&] (IndexSet & idxset)
+  {
+    deallog << "IndexSet before size=" << idxset.size() << " values: ";
+    idxset.print(deallog);
+    IndexSet back(idxset.make_trilinos_map(MPI_COMM_WORLD,true));
+    deallog << "IndexSet after size=" << back.size() << " values: ";
+    back.print(deallog);
+  };
+
+  deallog << "without overlap:" << std::endl;
+  check(set_my);
+
+  deallog << "with overlap:" << std::endl;
+  check(set_ghost);
+
+  TrilinosWrappers::MPI::Vector v;
+
+  v.reinit(set_my, set_ghost, MPI_COMM_WORLD);
+  IndexSet from_partitioner(v.vector_partitioner());
+  deallog << "vec size: " << v.size()
+          << " from_partitioner: " << from_partitioner.size()
+          << std::endl;
+
+  from_partitioner.print(deallog);
+}
+
+
+
+int main (int argc,char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+  MPILogInitAll log;
+  test();
+}

diff --git a/tests/trilinos/index_set_02.mpirun=2.output b/tests/trilinos/index_set_02.mpirun=2.output
new file mode 100644 (file)
index 0000000..22a1b65
--- /dev/null
+++ b/tests/trilinos/index_set_02.mpirun=2.output
@@ -0,0 +1,19 @@
+
+DEAL:0::without overlap:
+DEAL:0::IndexSet before size=100 values: {[0,49]}
+DEAL:0::IndexSet after size=100 values: {[0,49]}
+DEAL:0::with overlap:
+DEAL:0::IndexSet before size=100 values: {[0,49], [55,59]}
+DEAL:0::IndexSet after size=100 values: {[0,49], [55,59]}
+DEAL:0::vec size: 100 from_partitioner: 100
+DEAL:0::{[0,49], [55,59]}
+
+DEAL:1::without overlap:
+DEAL:1::IndexSet before size=100 values: {[50,99]}
+DEAL:1::IndexSet after size=100 values: {[50,99]}
+DEAL:1::with overlap:
+DEAL:1::IndexSet before size=100 values: {[45,99]}
+DEAL:1::IndexSet after size=100 values: {[45,99]}
+DEAL:1::vec size: 100 from_partitioner: 100
+DEAL:1::{[45,99]}
+