23#ifdef DEAL_II_WITH_TRILINOS
24# ifdef DEAL_II_WITH_MPI
25# include <Epetra_MpiComm.h>
27# include <Epetra_Map.h>
28# include <Epetra_SerialComm.h>
29# ifdef DEAL_II_TRILINOS_WITH_TPETRA
30# include <Tpetra_Map.hpp>
38#ifdef DEAL_II_WITH_TRILINOS
40# ifdef DEAL_II_TRILINOS_WITH_TPETRA
43 Teuchos::RCP<
const Tpetra::Map<int, types::signed_global_dof_index>> map)
45 , index_space_size(1 + map->getMaxAllGlobalIndex())
46 , largest_range(
numbers::invalid_unsigned_int)
48 Assert(map->getMinAllGlobalIndex() == 0,
50 "The Tpetra::Map does not contain the global index 0, "
51 "which means some entries are not present on any processor."));
54 if (map->isContiguous())
59# if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
60 const size_type n_indices = map->getLocalNumElements();
62 const size_type n_indices = map->getNodeNumElements();
65 map->getMyGlobalIndices().data();
77# ifdef DEAL_II_WITH_64BIT_INDICES
81 , index_space_size(1 + map.MaxAllGID64())
82 , largest_range(
numbers::invalid_unsigned_int)
84 Assert(map.MinAllGID64() == 0,
86 "The Epetra_BlockMap does not contain the global index 0, "
87 "which means some entries are not present on any processor."));
94 const size_type n_indices = map.NumMyElements();
96 reinterpret_cast<size_type *
>(map.MyGlobalElements64());
107 : is_compressed(true)
108 , index_space_size(1 + map.MaxAllGID())
109 , largest_range(
numbers::invalid_unsigned_int)
111 Assert(map.MinAllGID() == 0,
113 "The Epetra_BlockMap does not contain the global index 0, "
114 "which means some entries are not present on any processor."));
121 const size_type n_indices = map.NumMyElements();
122 unsigned int *indices =
123 reinterpret_cast<unsigned int *
>(map.MyGlobalElements());
150 std::vector<Range>::iterator store =
ranges.begin();
151 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();)
153 std::vector<Range>::iterator next = i;
160 while (next !=
ranges.end() && (next->begin <= last_index))
162 last_index =
std::max(last_index, next->end);
168 *store =
Range(first_index, last_index);
172 if (store !=
ranges.end())
174 std::vector<Range> new_ranges(
ranges.begin(), store);
179 size_type next_index = 0, largest_range_size = 0;
180 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();
185 i->nth_index_in_set = next_index;
186 next_index += (i->end - i->begin);
187 if (i->end - i->begin > largest_range_size)
189 largest_range_size = i->end - i->begin;
205 for (
const auto &range :
ranges)
207 ExcMessage(
"In the process of creating the current IndexSet "
208 "object, you added indices beyond the size of the index "
209 "space. Specifically, you added elements that form the "
211 std::to_string(range.begin) +
"," +
212 std::to_string(range.end) +
213 "), but the size of the index space is only " +
229 std::vector<Range>::const_iterator r1 =
ranges.begin(),
237 if (r1->end <= r2->begin)
239 else if (r2->end <= r1->begin)
244 Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
245 ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
249 result.add_range(
std::max(r1->begin, r2->begin),
255 if (r1->end <= r2->end)
273 ExcMessage(
"End index needs to be larger or equal to begin index!"));
275 ExcMessage(
"You are asking for a view into an IndexSet object "
276 "that would cover the sub-range [" +
277 std::to_string(
begin) +
',' + std::to_string(
end) +
278 "). But this is not a subset of the range "
279 "of the current object, which is [0," +
280 std::to_string(
size()) +
")."));
283 std::vector<Range>::const_iterator r1 =
ranges.begin();
285 while (r1 !=
ranges.end())
287 if ((r1->end >
begin) && (r1->begin <
end))
292 else if (r1->begin >=
end)
308 ExcMessage(
"The mask must have the same size index space "
309 "as the index set it is applied to."));
321 if (mask.ranges.size() == 1)
322 return get_view(mask.ranges[0].begin, mask.ranges[0].end);
331 std::vector<Range> new_ranges;
333 std::vector<Range>::iterator own_it =
ranges.begin();
334 std::vector<Range>::iterator mask_it = mask.ranges.begin();
336 while ((own_it !=
ranges.end()) && (mask_it != mask.ranges.end()))
342 if (own_it->end <= mask_it->begin)
350 if (mask_it->end <= own_it->begin)
369 if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
371 new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
372 own_it->end - mask_it->nth_index_in_set);
376 if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
378 const size_type offset_within_mask_interval =
379 own_it->begin - mask_it->begin;
380 new_ranges.emplace_back(mask_it->nth_index_in_set +
381 offset_within_mask_interval,
382 mask_it->nth_index_in_set +
383 (mask_it->end - mask_it->begin));
387 if ((own_it->begin <= mask_it->begin) &&
388 (own_it->end >= mask_it->end))
390 new_ranges.emplace_back(mask_it->begin -
391 mask_it->nth_index_in_set,
392 mask_it->end - mask_it->nth_index_in_set);
396 if ((mask_it->begin <= own_it->begin) &&
397 (mask_it->end >= own_it->end))
399 const size_type offset_within_mask_interval =
400 own_it->begin - mask_it->begin;
401 new_ranges.emplace_back(mask_it->nth_index_in_set +
402 offset_within_mask_interval,
403 mask_it->nth_index_in_set +
404 offset_within_mask_interval +
405 (own_it->end - own_it->begin));
414 if (own_it->end < mask_it->end)
416 else if (mask_it->end < own_it->end)
431 for (
const auto &range : new_ranges)
432 result.
add_range(range.begin, range.end);
442 const std::vector<types::global_dof_index> &n_indices_per_block)
const
444 std::vector<IndexSet> partitioned;
445 const unsigned int n_blocks = n_indices_per_block.size();
447 partitioned.reserve(n_blocks);
449 for (
const auto n_block_indices : n_indices_per_block)
451 partitioned.push_back(this->
get_view(start, start + n_block_indices));
452 start += n_block_indices;
457 for (
const auto &partition : partitioned)
459 sum += partition.size();
479 std::vector<Range> new_ranges;
481 std::vector<Range>::iterator own_it =
ranges.begin();
482 std::vector<Range>::iterator other_it = other.
ranges.begin();
484 while (own_it !=
ranges.end() && other_it != other.
ranges.end())
487 if (own_it->end <= other_it->begin)
489 new_ranges.push_back(*own_it);
494 if (own_it->begin >= other_it->end)
502 if (own_it->begin < other_it->begin)
504 Range r(own_it->begin, other_it->begin);
506 new_ranges.push_back(r);
511 own_it->begin = other_it->end;
512 if (own_it->begin > own_it->end)
514 own_it->begin = own_it->end;
523 for (; own_it !=
ranges.end(); ++own_it)
524 new_ranges.push_back(*own_it);
529 const std::vector<Range>::iterator
end = new_ranges.end();
530 for (std::vector<Range>::iterator it = new_ranges.begin(); it !=
end; ++it)
542 for (
const auto el : *
this)
543 set.add_indices(other, el * other.
size());
555 "pop_back() failed, because this IndexSet contains no entries."));
573 "pop_front() failed, because this IndexSet contains no entries."));
596 const auto insert_position =
598 if (insert_position ==
ranges.end() ||
599 insert_position->begin > new_range.
begin ||
600 insert_position->end < new_range.
end)
601 ranges.insert(insert_position, new_range);
608 boost::container::small_vector<std::pair<size_type, size_type>, 200>
610 const bool ranges_are_sorted)
612 if (!ranges_are_sorted)
613 std::sort(tmp_ranges.begin(), tmp_ranges.end());
623 if (tmp_ranges.size() > 9)
626 tmp_set.
ranges.reserve(tmp_ranges.size());
627 for (
const auto &i : tmp_ranges)
632 if (this->
ranges.size() <= 1)
634 if (this->
ranges.size() == 1)
636 std::swap(*
this, tmp_set);
642 for (
const auto &i : tmp_ranges)
651 if ((
this == &other) && (offset == 0))
654 if (other.
ranges.size() != 0)
662 std::vector<Range>::const_iterator r1 =
ranges.begin(),
663 r2 = other.
ranges.begin();
665 std::vector<Range> new_ranges;
672 if (r2 == other.
ranges.end() ||
673 (r1 !=
ranges.end() && r1->end < (r2->begin + offset)))
675 new_ranges.push_back(*r1);
678 else if (r1 ==
ranges.end() || (r2->end + offset) < r1->begin)
680 new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
688 std::max(r1->end, r2->end + offset));
689 new_ranges.push_back(next);
706 ExcMessage(
"One index set can only be a subset of another if they "
707 "describe index spaces of the same size. The ones in "
708 "question here have sizes " +
709 std::to_string(
size()) +
" and " +
710 std::to_string(other.
size()) +
"."));
730 out <<
size() <<
" ";
731 out <<
ranges.size() << std::endl;
732 std::vector<Range>::const_iterator r =
ranges.begin();
733 for (; r !=
ranges.end(); ++r)
735 out << r->begin <<
" " << r->end << std::endl;
747 unsigned int n_ranges;
752 for (
unsigned int i = 0; i < n_ranges; ++i)
769 std::size_t n_ranges =
ranges.size();
770 out.write(
reinterpret_cast<const char *
>(&n_ranges),
sizeof(n_ranges));
771 if (
ranges.empty() ==
false)
772 out.write(
reinterpret_cast<const char *
>(&*
ranges.begin()),
781 std::size_t n_ranges;
782 in.read(
reinterpret_cast<char *
>(&
size),
sizeof(
size));
783 in.read(
reinterpret_cast<char *
>(&n_ranges),
sizeof(n_ranges));
789 in.read(
reinterpret_cast<char *
>(&*
ranges.begin()),
812 std::vector<Range>::const_iterator p = std::upper_bound(
820 return ((index >= p->begin) && (index < p->end));
828 return (p->end > index);
845 return p->begin + (n - p->nth_index_in_set);
857 std::vector<Range>::const_iterator p =
861 if (p ==
ranges.end() || p->end == n || p->begin > n)
867 return (n - p->begin) + p->nth_index_in_set;
881 std::vector<Range>::const_iterator main_range =
884 Range r(global_index, global_index + 1);
887 std::vector<Range>::const_iterator range_begin, range_end;
888 if (global_index < main_range->
begin)
890 range_begin =
ranges.begin();
891 range_end = main_range;
895 range_begin = main_range;
901 const std::vector<Range>::const_iterator p =
914 if (global_index < p->
begin)
917 return {
this,
static_cast<size_type>(p -
ranges.begin()), global_index};
922std::vector<IndexSet::size_type>
927 std::vector<size_type> indices;
930 for (
const auto &range :
ranges)
931 for (
size_type entry = range.begin; entry < range.end; ++entry)
932 indices.push_back(entry);
949#ifdef DEAL_II_WITH_TRILINOS
950# ifdef DEAL_II_TRILINOS_WITH_TPETRA
952Tpetra::Map<int, types::signed_global_dof_index>
961Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
974 ExcMessage(
"You are trying to create an Tpetra::Map object "
975 "that partitions elements of an index set "
976 "between processors. However, the union of the "
977 "index sets on different processors does not "
978 "contain all indices exactly once: the sum of "
979 "the number of entries the various processors "
980 "want to store locally is " +
981 std::to_string(n_global_elements) +
982 " whereas the total size of the object to be "
984 std::to_string(
size()) +
985 ". In other words, there are "
986 "either indices that are not spoken for "
987 "by any processor, or there are indices that are "
988 "claimed by multiple processors."));
998 Tpetra::Map<int, types::signed_global_dof_index>>(
1002# ifdef DEAL_II_WITH_MPI
1003 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1012 std::vector<types::signed_global_dof_index> int_indices(indices.size());
1013 std::copy(indices.begin(), indices.end(), int_indices.begin());
1014 const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
1018 Tpetra::Map<int, types::signed_global_dof_index>>(
1022# ifdef DEAL_II_WITH_MPI
1023 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1048 ExcMessage(
"You are trying to create an Epetra_Map object "
1049 "that partitions elements of an index set "
1050 "between processors. However, the union of the "
1051 "index sets on different processors does not "
1052 "contain all indices exactly once: the sum of "
1053 "the number of entries the various processors "
1054 "want to store locally is " +
1055 std::to_string(n_global_elements) +
1056 " whereas the total size of the object to be "
1058 std::to_string(
size()) +
1059 ". In other words, there are "
1060 "either indices that are not spoken for "
1061 "by any processor, or there are indices that are "
1062 "claimed by multiple processors."));
1075# ifdef DEAL_II_WITH_MPI
1076 Epetra_MpiComm(communicator)
1092# ifdef DEAL_II_WITH_MPI
1093 Epetra_MpiComm(communicator)
1103#ifdef DEAL_II_WITH_PETSC
1107 std::vector<size_type> indices;
1110 PetscInt n = indices.size();
1111 std::vector<PetscInt> pindices(indices.begin(), indices.end());
1114 PetscErrorCode ierr =
1115 ISCreateGeneral(communicator, n, pindices.data(), PETSC_COPY_VALUES, &is);
1131 if (n_global_elements !=
size())
1134 if (n_global_elements == 0)
1137#ifdef DEAL_II_WITH_MPI
1139 const bool all_contiguous =
1141 if (!all_contiguous)
1144 bool is_globally_ascending =
true;
1151 const std::vector<types::global_dof_index> global_dofs =
1161 for (; index < global_dofs.size(); ++index)
1166 if (new_dof <= old_dof)
1168 is_globally_ascending =
false;
1178 int is_ascending = is_globally_ascending ? 1 : 0;
1179 int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
1182 return (is_ascending == 1);
ElementIterator begin() const
bool is_subset_of(const IndexSet &other) const
bool is_element_binary_search(const size_type local_index) const
size_type index_within_set_binary_search(const size_type global_index) const
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
bool is_contiguous() const
ElementIterator at(const size_type global_index) const
void fill_index_vector(std::vector< size_type > &indices) const
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
size_type n_elements() const
void add_range_lower_bound(const Range &range)
ElementIterator begin() const
void set_size(const size_type size)
void read(std::istream &in)
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
IndexSet tensor_product(const IndexSet &other) const
void write(std::ostream &out) const
void block_read(std::istream &in)
void add_ranges_internal(boost::container::small_vector< std::pair< size_type, size_type >, 200 > &tmp_ranges, const bool ranges_are_sorted)
IntervalIterator begin_intervals() const
Tpetra::Map< int, types::signed_global_dof_index > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
std::vector< Range > ranges
void subtract_set(const IndexSet &other)
ElementIterator end() const
Threads::Mutex compress_mutex
size_type index_space_size
void block_write(std::ostream &out) const
IndexSet get_view(const size_type begin, const size_type end) const
void add_range(const size_type begin, const size_type end)
std::size_t memory_consumption() const
size_type nth_index_in_set_binary_search(const size_type local_index) 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
std::vector< size_type > get_index_vector() const
types::global_dof_index size_type
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
IndexSet operator&(const IndexSet &is) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
T sum(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
Teuchos::RCP< T > make_rcp(Args &&...args)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
const types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
size_type nth_index_in_set