From 89071088c086880162bd370613ed0c604b518ff8 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 15 Aug 2017 14:23:15 -0600 Subject: [PATCH] Fix incorrect size in overlapping EPetra_Map to IndexSet --- doc/news/changes/minor/20170815Heister | 3 + source/base/index_set.cc | 20 +++-- tests/trilinos/index_set_02.cc | 84 +++++++++++++++++++++ tests/trilinos/index_set_02.mpirun=2.output | 19 +++++ 4 files changed, 121 insertions(+), 5 deletions(-) create mode 100644 doc/news/changes/minor/20170815Heister create mode 100644 tests/trilinos/index_set_02.cc create mode 100644 tests/trilinos/index_set_02.mpirun=2.output diff --git a/doc/news/changes/minor/20170815Heister b/doc/news/changes/minor/20170815Heister new file mode 100644 index 0000000000..0d0822138a --- /dev/null +++ b/doc/news/changes/minor/20170815Heister @@ -0,0 +1,3 @@ +Fixed: Constructing an IndexSet from an Epetra_Map generated the wrong size. +
+(Timo Heister, 2017/08/15) diff --git a/source/base/index_set.cc b/source/base/index_set.cc index 4ed7b59b3a..9543531e65 100644 --- 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(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(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 index 0000000000..49fd75fe48 --- /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 +#include +#include + + +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 index 0000000000..22a1b652ad --- /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]} + -- 2.39.5