24 #ifdef DEAL_II_WITH_TRILINOS
25 # ifdef DEAL_II_WITH_MPI
26 # include <Epetra_MpiComm.h>
28 # include <Epetra_Map.h>
29 # include <Epetra_SerialComm.h>
30 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
31 # include <Tpetra_Map.hpp>
39 #ifdef DEAL_II_WITH_TRILINOS
41 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
44 Teuchos::RCP<
const Tpetra::Map<int, types::signed_global_dof_index>> map)
46 , index_space_size(1 + map->getMaxAllGlobalIndex())
49 Assert(map->getMinAllGlobalIndex() == 0,
51 "The Tpetra::Map does not contain the global index 0, "
52 "which means some entries are not present on any processor."));
55 if (map->isContiguous())
60 # if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
61 const size_type n_indices = map->getLocalNumElements();
63 const size_type n_indices = map->getNodeNumElements();
66 map->getMyGlobalIndices().data();
78 # ifdef DEAL_II_WITH_64BIT_INDICES
82 , index_space_size(1 + map.MaxAllGID64())
85 Assert(map.MinAllGID64() == 0,
87 "The Epetra_BlockMap does not contain the global index 0, "
88 "which means some entries are not present on any processor."));
95 const size_type n_indices = map.NumMyElements();
97 reinterpret_cast<size_type *
>(map.MyGlobalElements64());
108 : is_compressed(true)
109 , index_space_size(1 + map.MaxAllGID())
112 Assert(map.MinAllGID() == 0,
114 "The Epetra_BlockMap does not contain the global index 0, "
115 "which means some entries are not present on any processor."));
122 const size_type n_indices = map.NumMyElements();
123 unsigned int *indices =
124 reinterpret_cast<unsigned int *
>(map.MyGlobalElements());
151 std::vector<Range>::iterator store =
ranges.begin();
152 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();)
154 std::vector<Range>::iterator next = i;
161 while (next !=
ranges.end() && (next->begin <= last_index))
163 last_index =
std::max(last_index, next->end);
169 *store =
Range(first_index, last_index);
173 if (store !=
ranges.end())
175 std::vector<Range> new_ranges(
ranges.begin(), store);
180 size_type next_index = 0, largest_range_size = 0;
181 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();
186 i->nth_index_in_set = next_index;
187 next_index += (i->end - i->begin);
188 if (i->end - i->begin > largest_range_size)
190 largest_range_size = i->end - i->begin;
206 for (
const auto &range :
ranges)
208 ExcMessage(
"In the process of creating the current IndexSet "
209 "object, you added indices beyond the size of the index "
210 "space. Specifically, you added elements that form the "
214 "), but the size of the index space is only " +
230 std::vector<Range>::const_iterator r1 =
ranges.begin(),
238 if (r1->end <= r2->begin)
240 else if (r2->end <= r1->begin)
245 Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
246 ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
250 result.add_range(
std::max(r1->begin, r2->begin),
256 if (r1->end <= r2->end)
274 ExcMessage(
"End index needs to be larger or equal to begin index!"));
276 ExcMessage(
"You are asking for a view into an IndexSet object "
277 "that would cover the sub-range [" +
279 "). But this is not a subset of the range "
280 "of the current object, which is [0," +
284 std::vector<Range>::const_iterator r1 =
ranges.begin();
286 while (r1 !=
ranges.end())
288 if ((r1->end >
begin) && (r1->begin <
end))
293 else if (r1->begin >=
end)
309 ExcMessage(
"The mask must have the same size index space "
310 "as the index set it is applied to."));
322 if (mask.ranges.size() == 1)
323 return get_view(mask.ranges[0].begin, mask.ranges[0].end);
332 std::vector<Range> new_ranges;
334 std::vector<Range>::iterator own_it =
ranges.begin();
335 std::vector<Range>::iterator mask_it = mask.ranges.begin();
337 while ((own_it !=
ranges.end()) && (mask_it != mask.ranges.end()))
343 if (own_it->end <= mask_it->begin)
351 if (mask_it->end <= own_it->begin)
370 if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
372 new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
373 own_it->end - mask_it->nth_index_in_set);
377 if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
379 const size_type offset_within_mask_interval =
380 own_it->begin - mask_it->begin;
381 new_ranges.emplace_back(mask_it->nth_index_in_set +
382 offset_within_mask_interval,
383 mask_it->nth_index_in_set +
384 (mask_it->end - mask_it->begin));
388 if ((own_it->begin <= mask_it->begin) &&
389 (own_it->end >= mask_it->end))
391 new_ranges.emplace_back(mask_it->begin -
392 mask_it->nth_index_in_set,
393 mask_it->end - mask_it->nth_index_in_set);
397 if ((mask_it->begin <= own_it->begin) &&
398 (mask_it->end >= own_it->end))
400 const size_type offset_within_mask_interval =
401 own_it->begin - mask_it->begin;
402 new_ranges.emplace_back(mask_it->nth_index_in_set +
403 offset_within_mask_interval,
404 mask_it->nth_index_in_set +
405 offset_within_mask_interval +
406 (own_it->end - own_it->begin));
415 if (own_it->end < mask_it->end)
417 else if (mask_it->end < own_it->end)
432 for (
const auto &range : new_ranges)
433 result.
add_range(range.begin, range.end);
441 std::vector<IndexSet>
443 const std::vector<types::global_dof_index> &n_indices_per_block)
const
445 std::vector<IndexSet> partitioned;
446 const unsigned int n_blocks = n_indices_per_block.size();
450 for (
const auto n_block_indices : n_indices_per_block)
452 partitioned.push_back(this->
get_view(start, start + n_block_indices));
453 start += n_block_indices;
458 for (
const auto &
partition : partitioned)
480 std::vector<Range> new_ranges;
482 std::vector<Range>::iterator own_it =
ranges.begin();
483 std::vector<Range>::iterator other_it = other.
ranges.begin();
485 while (own_it !=
ranges.end() && other_it != other.
ranges.end())
488 if (own_it->end <= other_it->begin)
490 new_ranges.push_back(*own_it);
495 if (own_it->begin >= other_it->end)
503 if (own_it->begin < other_it->begin)
505 Range r(own_it->begin, other_it->begin);
507 new_ranges.push_back(r);
512 own_it->begin = other_it->end;
513 if (own_it->begin > own_it->end)
515 own_it->begin = own_it->end;
524 for (; own_it !=
ranges.end(); ++own_it)
525 new_ranges.push_back(*own_it);
530 const std::vector<Range>::iterator
end = new_ranges.end();
531 for (std::vector<Range>::iterator it = new_ranges.begin(); it !=
end; ++it)
543 for (
const auto el : *
this)
544 set.add_indices(other, el * other.
size());
556 "pop_back() failed, because this IndexSet contains no entries."));
574 "pop_front() failed, because this IndexSet contains no entries."));
597 const auto insert_position =
599 if (insert_position ==
ranges.end() ||
600 insert_position->begin > new_range.
begin ||
601 insert_position->end < new_range.
end)
602 ranges.insert(insert_position, new_range);
609 boost::container::small_vector<std::pair<size_type, size_type>, 200>
611 const bool ranges_are_sorted)
613 if (!ranges_are_sorted)
614 std::sort(tmp_ranges.begin(), tmp_ranges.end());
624 if (tmp_ranges.size() > 9)
627 tmp_set.
ranges.reserve(tmp_ranges.size());
628 for (
const auto &i : tmp_ranges)
633 if (this->
ranges.size() <= 1)
635 if (this->
ranges.size() == 1)
643 for (
const auto &i : tmp_ranges)
652 if ((
this == &other) && (offset == 0))
655 if (other.
ranges.size() != 0)
663 std::vector<Range>::const_iterator r1 =
ranges.begin(),
664 r2 = other.
ranges.begin();
666 std::vector<Range> new_ranges;
673 if (r2 == other.
ranges.end() ||
674 (r1 !=
ranges.end() && r1->end < (r2->begin + offset)))
676 new_ranges.push_back(*r1);
679 else if (r1 ==
ranges.end() || (r2->end + offset) < r1->begin)
681 new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
689 std::max(r1->end, r2->end + offset));
690 new_ranges.push_back(next);
707 ExcMessage(
"One index set can only be a subset of another if they "
708 "describe index spaces of the same size. The ones in "
709 "question here have sizes " +
731 out <<
size() <<
" ";
732 out <<
ranges.size() << std::endl;
733 std::vector<Range>::const_iterator r =
ranges.begin();
734 for (; r !=
ranges.end(); ++r)
736 out << r->begin <<
" " << r->end << std::endl;
748 unsigned int n_ranges;
753 for (
unsigned int i = 0; i < n_ranges; ++i)
770 std::size_t n_ranges =
ranges.size();
771 out.write(
reinterpret_cast<const char *
>(&n_ranges),
sizeof(n_ranges));
772 if (
ranges.empty() ==
false)
773 out.write(
reinterpret_cast<const char *
>(&*
ranges.begin()),
782 std::size_t n_ranges;
783 in.read(
reinterpret_cast<char *
>(&
size),
sizeof(
size));
784 in.read(
reinterpret_cast<char *
>(&n_ranges),
sizeof(n_ranges));
790 in.read(
reinterpret_cast<char *
>(&*
ranges.begin()),
813 std::vector<Range>::const_iterator p = std::upper_bound(
821 return ((index >= p->begin) && (index < p->
end));
829 return (p->end > index);
846 return p->begin + (n - p->nth_index_in_set);
858 std::vector<Range>::const_iterator p =
862 if (p ==
ranges.end() || p->end == n || p->begin > n)
868 return (n - p->begin) + p->nth_index_in_set;
882 std::vector<Range>::const_iterator main_range =
888 std::vector<Range>::const_iterator range_begin, range_end;
889 if (global_index < main_range->
begin)
891 range_begin =
ranges.begin();
892 range_end = main_range;
896 range_begin = main_range;
902 const std::vector<Range>::const_iterator p =
915 if (global_index < p->
begin)
923 std::vector<IndexSet::size_type>
928 std::vector<size_type> indices;
931 for (
const auto &range :
ranges)
932 for (
size_type entry = range.begin; entry < range.end; ++entry)
933 indices.push_back(entry);
950 #ifdef DEAL_II_WITH_TRILINOS
951 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
953 Tpetra::Map<int, types::signed_global_dof_index>
955 const bool overlapping)
const
962 Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
964 const bool overlapping)
const
975 ExcMessage(
"You are trying to create an Tpetra::Map object "
976 "that partitions elements of an index set "
977 "between processors. However, the union of the "
978 "index sets on different processors does not "
979 "contain all indices exactly once: the sum of "
980 "the number of entries the various processors "
981 "want to store locally is " +
983 " whereas the total size of the object to be "
986 ". In other words, there are "
987 "either indices that are not spoken for "
988 "by any processor, or there are indices that are "
989 "claimed by multiple processors."));
999 Tpetra::Map<int, types::signed_global_dof_index>>(
1003 # ifdef DEAL_II_WITH_MPI
1004 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1013 std::vector<types::signed_global_dof_index> int_indices(indices.size());
1014 std::copy(indices.begin(), indices.end(), int_indices.begin());
1015 const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
1019 Tpetra::Map<int, types::signed_global_dof_index>>(
1023 # ifdef DEAL_II_WITH_MPI
1024 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1038 const bool overlapping)
const
1049 ExcMessage(
"You are trying to create an Epetra_Map object "
1050 "that partitions elements of an index set "
1051 "between processors. However, the union of the "
1052 "index sets on different processors does not "
1053 "contain all indices exactly once: the sum of "
1054 "the number of entries the various processors "
1055 "want to store locally is " +
1057 " whereas the total size of the object to be "
1060 ". In other words, there are "
1061 "either indices that are not spoken for "
1062 "by any processor, or there are indices that are "
1063 "claimed by multiple processors."));
1076 # ifdef DEAL_II_WITH_MPI
1077 Epetra_MpiComm(communicator)
1093 # ifdef DEAL_II_WITH_MPI
1094 Epetra_MpiComm(communicator)
1104 #ifdef DEAL_II_WITH_PETSC
1108 std::vector<size_type> indices;
1111 PetscInt n = indices.size();
1112 std::vector<PetscInt> pindices(indices.begin(), indices.end());
1115 PetscErrorCode ierr =
1116 ISCreateGeneral(communicator, n, pindices.data(), PETSC_COPY_VALUES, &is);
1138 #ifdef DEAL_II_WITH_MPI
1140 const bool all_contiguous =
1142 if (!all_contiguous)
1145 bool is_globally_ascending =
true;
1152 const std::vector<types::global_dof_index> global_dofs =
1162 for (; index < global_dofs.size(); ++index)
1167 if (new_dof <= old_dof)
1169 is_globally_ascending =
false;
1179 int is_ascending = is_globally_ascending ? 1 : 0;
1180 int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
1183 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
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Teuchos::RCP< T > make_rcp(Args &&...args)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
const types::global_dof_index invalid_dof_index
int signed_global_dof_index
unsigned int global_dof_index
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