deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mpi_remote_point_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_mpi_mpi_remote_point_evaluation_h
16#define dealii_mpi_mpi_remote_point_evaluation_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/mpi.h>
21#include <deal.II/base/mpi_tags.h>
22
23#include <deal.II/dofs/dof_handler.h>
24
25DEAL_II_NAMESPACE_OPEN
26
27// Forward declarations
28namespace GridTools
29{
30 template <int dim, int spacedim>
31 class Cache;
32
33 namespace internal
34 {
35 template <int dim, int spacedim>
36 struct DistributedComputePointLocationsInternal;
37 }
38} // namespace GridTools
39
40namespace Utilities
41{
42 namespace MPI
43 {
54 template <int dim, int spacedim = dim>
55 class RemotePointEvaluation
56 {
57 public:
62 struct AdditionalData
63 {
64 public:
68 AdditionalData(
69 const double tolerance = 1e-6,
70 const bool enforce_unique_mapping = false,
71 const unsigned int rtree_level = 0,
72 const std::function<std::vector<bool>()> &marked_vertices = {});
73
82 double tolerance;
83
88 bool enforce_unique_mapping;
89
93 unsigned int rtree_level;
94
99 std::function<std::vector<bool>()> marked_vertices;
100 };
101
107 RemotePointEvaluation(
108 const AdditionalData &additional_data = AdditionalData());
109
128 DEAL_II_DEPRECATED_WITH_COMMENT(
129 "Use the constructor with AdditionalData struct.")
130 RemotePointEvaluation(
131 const double tolerance,
132 const bool enforce_unique_mapping = false,
133 const unsigned int rtree_level = 0,
134 const std::function<std::vector<bool>()> &marked_vertices = {});
135
139 ~RemotePointEvaluation();
140
152 void
153 reinit(const std::vector<Point<spacedim>> &points,
154 const Triangulation<dim, spacedim> &tria,
155 const Mapping<dim, spacedim> &mapping);
156
167 void
168 reinit(const GridTools::Cache<dim, spacedim> &cache,
169 const std::vector<Point<spacedim>> &points);
170
180 void
181 reinit(const GridTools::internal::
182 DistributedComputePointLocationsInternal<dim, spacedim> &data,
183 const Triangulation<dim, spacedim> &tria,
184 const Mapping<dim, spacedim> &mapping);
185
190 class CellData
191 {
192 public:
198 CellData(const Triangulation<dim, spacedim> &triangulation);
199
225 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
226 cell_indices() const;
227
231 typename Triangulation<dim, spacedim>::active_cell_iterator
232 get_active_cell_iterator(const unsigned int cell) const;
233
237 ArrayView<const Point<dim>>
238 get_unit_points(const unsigned int cell) const;
239
243 template <typename T>
244 ArrayView<T>
245 get_data_view(const unsigned int cell,
246 const ArrayView<T> &values) const;
247
251 std::vector<std::pair<int, int>> cells;
252
257 std::vector<unsigned int> reference_point_ptrs;
258
262 std::vector<Point<dim>> reference_point_values;
263
264 private:
269 const Triangulation<dim, spacedim> &triangulation;
270 };
271
276 const CellData &
277 get_cell_data() const;
278
299 template <typename T, unsigned int n_components = 1>
300 void
301 evaluate_and_process(
302 std::vector<T> &output,
303 std::vector<T> &buffer,
304 const std::function<void(const ArrayView<T> &, const CellData &)>
305 &evaluation_function,
306 const bool sort_data = true) const;
307
312 template <typename T, unsigned int n_components = 1>
313 std::vector<T>
314 evaluate_and_process(
315 const std::function<void(const ArrayView<T> &, const CellData &)>
316 &evaluation_function,
317 const bool sort_data = true) const;
318
333 template <typename T, unsigned int n_components = 1>
334 void
335 process_and_evaluate(
336 const std::vector<T> &input,
337 std::vector<T> &buffer,
338 const std::function<void(const ArrayView<const T> &, const CellData &)>
339 &evaluation_function,
340 const bool sort_data = true) const;
341
346 template <typename T, unsigned int n_components = 1>
347 void
348 process_and_evaluate(
349 const std::vector<T> &input,
350 const std::function<void(const ArrayView<const T> &, const CellData &)>
351 &evaluation_function,
352 const bool sort_data = true) const;
353
358 const std::vector<unsigned int> &
359 get_point_ptrs() const;
360
367 bool
368 is_map_unique() const;
369
373 bool
374 all_points_found() const;
375
379 bool
380 point_found(const unsigned int i) const;
381
385 const Triangulation<dim, spacedim> &
386 get_triangulation() const;
387
391 const Mapping<dim, spacedim> &
392 get_mapping() const;
393
399 bool
400 is_ready() const;
401
407 const std::vector<unsigned int> &
408 get_send_permutation() const;
409
416 const std::vector<unsigned int> &
417 get_inverse_recv_permutation() const;
418
419
420 private:
424 const AdditionalData additional_data;
425
429 boost::signals2::connection tria_signal;
430
436 bool ready_flag;
437
441 ObserverPointer<const Triangulation<dim, spacedim>> tria;
442
446 ObserverPointer<const Mapping<dim, spacedim>> mapping;
447
451 bool unique_mapping;
452
456 bool all_points_found_flag;
457
463 std::vector<unsigned int> point_ptrs;
464
468 std::vector<unsigned int> recv_permutation;
469
473 std::vector<unsigned int> recv_permutation_inv;
474
479 std::vector<unsigned int> recv_ptrs;
480
484 std::vector<unsigned int> recv_ranks;
485
490 std::unique_ptr<CellData> cell_data;
491
495 std::vector<unsigned int> send_permutation;
496
500 std::vector<unsigned int> send_permutation_inv;
501
505 std::vector<unsigned int> send_ranks;
506
511 std::vector<unsigned int> send_ptrs;
512
526 unsigned int buffer_size_with_sorting;
527
532 unsigned int buffer_size_without_sorting;
533 };
534
535 namespace internal
536 {
537#ifdef DEAL_II_WITH_MPI
541 template <typename T>
542 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == false, void>
543 pack_and_isend(const ArrayView<const T> &data,
544 const unsigned int rank,
545 const unsigned int tag,
546 const MPI_Comm comm,
547 std::vector<std::vector<char>> &buffers,
548 std::vector<MPI_Request> &requests)
549 {
550 requests.emplace_back(MPI_Request());
551
552 buffers.emplace_back(Utilities::pack(
553 std::vector<T>(data.data(), data.data() + data.size()), false));
554
555 const int ierr = MPI_Isend(buffers.back().data(),
556 buffers.back().size(),
557 MPI_CHAR,
558 rank,
559 tag,
560 comm,
561 &requests.back());
562 AssertThrowMPI(ierr);
563 }
564
565
566
571 template <typename T>
572 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == true, void>
573 pack_and_isend(const ArrayView<const T> &data,
574 const unsigned int rank,
575 const unsigned int tag,
576 const MPI_Comm comm,
577 std::vector<std::vector<char>> & /*buffers*/,
578 std::vector<MPI_Request> &requests)
579 {
580 requests.emplace_back(MPI_Request());
581
582 const int ierr = MPI_Isend(data.data(),
583 data.size(),
584 Utilities::MPI::mpi_type_id_for_type<T>,
585 rank,
586 tag,
587 comm,
588 &requests.back());
589 AssertThrowMPI(ierr);
590 }
591
592
593
598 template <int rank_, int dim, typename T>
599 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == true, void>
600 pack_and_isend(const ArrayView<const Tensor<rank_, dim, T>> &data,
601 const unsigned int rank,
602 const unsigned int tag,
603 const MPI_Comm comm,
604 std::vector<std::vector<char>> &buffers,
605 std::vector<MPI_Request> &requests)
606 {
607 ArrayView<const T> data_(reinterpret_cast<const T *>(data.data()),
608 data.size() * Utilities::pow(dim, rank_));
609
610 pack_and_isend(data_, rank, tag, comm, buffers, requests);
611 }
612
613
617 template <typename T>
618 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == false, void>
619 recv_and_unpack(const ArrayView<T> &data,
620 const MPI_Comm comm,
621 const MPI_Status &status,
622 std::vector<char> &buffer)
623 {
624 int message_length;
625 int ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
626 AssertThrowMPI(ierr);
627
628 buffer.resize(message_length);
629
630 ierr = MPI_Recv(buffer.data(),
631 buffer.size(),
632 MPI_CHAR,
633 status.MPI_SOURCE,
634 internal::Tags::remote_point_evaluation,
635 comm,
636 MPI_STATUS_IGNORE);
637 AssertThrowMPI(ierr);
638
639 // unpack data
640 const auto temp = Utilities::unpack<std::vector<T>>(buffer, false);
641
642 for (unsigned int i = 0; i < data.size(); ++i)
643 data[i] = temp[i];
644 }
645
646
647
652 template <typename T>
653 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == true, void>
654 recv_and_unpack(const ArrayView<T> &data,
655 const MPI_Comm comm,
656 const MPI_Status &status,
657 std::vector<char> & /*buffer*/)
658 {
659 const auto ierr = MPI_Recv(data.data(),
660 data.size(),
661 Utilities::MPI::mpi_type_id_for_type<T>,
662 status.MPI_SOURCE,
663 internal::Tags::remote_point_evaluation,
664 comm,
665 MPI_STATUS_IGNORE);
666 AssertThrowMPI(ierr);
667 }
668
669
670
675 template <int rank_, int dim, typename T>
676 std::enable_if_t<Utilities::MPI::is_mpi_type<T> == true, void>
677 recv_and_unpack(const ArrayView<Tensor<rank_, dim, T>> &data,
678 const MPI_Comm comm,
679 const MPI_Status &status,
680 std::vector<char> &buffer)
681 {
682 const ArrayView<T> data_(reinterpret_cast<T *>(data.data()),
683 data.size() * Utilities::pow(dim, rank_));
684
685 recv_and_unpack(data_, comm, status, buffer);
686 }
687#endif
688 } // namespace internal
689
690
691
692 template <int dim, int spacedim>
693 template <typename T>
694 ArrayView<T>
695 RemotePointEvaluation<dim, spacedim>::CellData::get_data_view(
696 const unsigned int cell,
697 const ArrayView<T> &values) const
698 {
699 AssertIndexRange(cell, cells.size());
700 return {values.data() + reference_point_ptrs[cell],
701 reference_point_ptrs[cell + 1] - reference_point_ptrs[cell]};
702 }
703
704
705
706 template <int dim, int spacedim>
707 template <typename T, unsigned int n_components>
708 void
709 RemotePointEvaluation<dim, spacedim>::evaluate_and_process(
710 std::vector<T> &output,
711 std::vector<T> &buffer,
712 const std::function<void(const ArrayView<T> &, const CellData &)>
713 &evaluation_function,
714 const bool sort_data) const
715 {
716#ifndef DEAL_II_WITH_MPI
717 Assert(false, ExcNeedsMPI());
718 (void)output;
719 (void)buffer;
720 (void)evaluation_function;
721 (void)sort_data;
722#else
723 static CollectiveMutex mutex;
724 CollectiveMutex::ScopedLock lock(mutex, tria->get_mpi_communicator());
725
726 const unsigned int my_rank =
727 Utilities::MPI::this_mpi_process(tria->get_mpi_communicator());
728
729 // allocate memory for output and buffer
730 output.resize(point_ptrs.back() * n_components);
731
732 buffer.resize(
733 (sort_data ? buffer_size_with_sorting : buffer_size_without_sorting) *
734 n_components);
735
736 // ... for evaluation
737 ArrayView<T> buffer_eval(buffer.data(),
738 send_permutation.size() * n_components);
739
740 // ... for communication (send)
741 ArrayView<T> buffer_send(sort_data ?
742 (buffer.data() +
743 send_permutation.size() * n_components) :
744 buffer.data(),
745 send_permutation.size() * n_components);
746
747 // more arrays
748 std::vector<MPI_Request> send_requests;
749 std::vector<std::vector<char>> send_buffers_packed;
750 std::vector<char> recv_buffer_packed;
751
752 // evaluate functions at points
753 evaluation_function(buffer_eval, *cell_data);
754
755 // sort for communication (optional)
756 if (sort_data)
757 {
758 const auto my_rank_local_recv_ptr =
759 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
760
761 if (my_rank_local_recv_ptr != recv_ranks.end())
762 {
763 const unsigned int my_rank_local_recv =
764 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
765 const unsigned int my_rank_local_send = std::distance(
766 send_ranks.begin(),
767 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
768 const unsigned int start = send_ptrs[my_rank_local_send];
769 const unsigned int end = send_ptrs[my_rank_local_send + 1];
770 const unsigned int *recv_ptr =
771 recv_permutation.data() + recv_ptrs[my_rank_local_recv];
772 for (unsigned int i = 0; i < send_permutation.size(); ++i)
773 {
774 const unsigned int send_index = send_permutation[i];
775
776 if (start <= send_index && send_index < end)
777 // local data -> can be copied to output directly
778 for (unsigned int c = 0; c < n_components; ++c)
779 output[recv_ptr[send_index - start] * n_components + c] =
780 buffer_eval[i * n_components + c];
781 else
782 // data to be sent
783 for (unsigned int c = 0; c < n_components; ++c)
784 buffer_send[send_index * n_components + c] =
785 buffer_eval[i * n_components + c];
786 }
787 }
788 else
789 {
790 for (unsigned int i = 0; i < send_permutation.size(); ++i)
791 for (unsigned int c = 0; c < n_components; ++c)
792 buffer_send[send_permutation[i] * n_components + c] =
793 buffer_eval[i * n_components + c];
794 }
795 }
796
797 // send data
798 send_buffers_packed.reserve(send_ranks.size());
799 send_requests.reserve(send_ranks.size());
800
801 for (unsigned int i = 0; i < send_ranks.size(); ++i)
802 {
803 if (send_ranks[i] == my_rank)
804 continue;
805
806 internal::pack_and_isend(
807 ArrayView<const T>(
808 buffer_send.begin() + send_ptrs[i] * n_components,
809 (send_ptrs[i + 1] - send_ptrs[i]) * n_components),
810 send_ranks[i],
811 internal::Tags::remote_point_evaluation,
812 tria->get_mpi_communicator(),
813 send_buffers_packed,
814 send_requests);
815 }
816
817 // copy local data directly to output if no sorting is requested
818 if (!sort_data)
819 {
820 const auto my_rank_local_recv_ptr =
821 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
822
823 if (my_rank_local_recv_ptr != recv_ranks.end())
824 {
825 const unsigned int my_rank_local_recv =
826 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
827 const unsigned int my_rank_local_send = std::distance(
828 send_ranks.begin(),
829 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
830
831 for (unsigned int j = recv_ptrs[my_rank_local_recv],
832 k = send_ptrs[my_rank_local_send];
833 j < recv_ptrs[my_rank_local_recv + 1];
834 ++j, ++k)
835 for (unsigned int c = 0; c < n_components; ++c)
836 output[j * n_components + c] =
837 buffer_eval[k * n_components + c];
838 }
839 }
840
841 // receive data
842 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
843 {
844 if (recv_ranks[i] == my_rank)
845 continue;
846
847 // receive remote data
848 MPI_Status status;
849 int ierr = MPI_Probe(MPI_ANY_SOURCE,
850 internal::Tags::remote_point_evaluation,
851 tria->get_mpi_communicator(),
852 &status);
853 AssertThrowMPI(ierr);
854
855 const auto ptr =
856 std::find(recv_ranks.begin(), recv_ranks.end(), status.MPI_SOURCE);
857
858 Assert(ptr != recv_ranks.end(), ExcNotImplemented());
859
860 const unsigned int j = std::distance(recv_ranks.begin(), ptr);
861
862 // ... for communication (recv)
863 ArrayView<T> recv_buffer(
864 sort_data ?
865 (buffer.data() + send_permutation.size() * 2 * n_components) :
866 (output.data() + recv_ptrs[j] * n_components),
867 (recv_ptrs[j + 1] - recv_ptrs[j]) * n_components);
868
869 internal::recv_and_unpack(recv_buffer,
870 tria->get_mpi_communicator(),
871 status,
872 recv_buffer_packed);
873
874 if (sort_data)
875 {
876 // sort data into output vector (optional)
877 for (unsigned int i = recv_ptrs[j], k = 0; i < recv_ptrs[j + 1];
878 ++i, ++k)
879 for (unsigned int c = 0; c < n_components; ++c)
880 output[recv_permutation[i] * n_components + c] =
881 recv_buffer[k * n_components + c];
882 }
883 }
884
885 // make sure all messages have been sent
886 if (!send_requests.empty())
887 {
888 const int ierr = MPI_Waitall(send_requests.size(),
889 send_requests.data(),
890 MPI_STATUSES_IGNORE);
891 AssertThrowMPI(ierr);
892 }
893#endif
894 }
895
896
897 template <int dim, int spacedim>
898 template <typename T, unsigned int n_components>
899 std::vector<T>
900 RemotePointEvaluation<dim, spacedim>::evaluate_and_process(
901 const std::function<void(const ArrayView<T> &, const CellData &)>
902 &evaluation_function,
903 const bool sort_data) const
904 {
905 std::vector<T> output;
906 std::vector<T> buffer;
907
908 this->evaluate_and_process<T, n_components>(output,
909 buffer,
910 evaluation_function,
911 sort_data);
912
913 return output;
914 }
915
916
917
918 template <int dim, int spacedim>
919 template <typename T, unsigned int n_components>
920 void
921 RemotePointEvaluation<dim, spacedim>::process_and_evaluate(
922 const std::vector<T> &input,
923 std::vector<T> &buffer,
924 const std::function<void(const ArrayView<const T> &, const CellData &)>
925 &evaluation_function,
926 const bool sort_data) const
927 {
928#ifndef DEAL_II_WITH_MPI
929 Assert(false, ExcNeedsMPI());
930 (void)input;
931 (void)buffer;
932 (void)evaluation_function;
933 (void)sort_data;
934#else
935 static CollectiveMutex mutex;
936 CollectiveMutex::ScopedLock lock(mutex, tria->get_mpi_communicator());
937
938 const unsigned int my_rank =
939 Utilities::MPI::this_mpi_process(tria->get_mpi_communicator());
940
941 // allocate memory for buffer
942 const auto &point_ptrs = this->get_point_ptrs();
943
944 AssertDimension(input.size(),
945 sort_data ? ((point_ptrs.size() - 1) * n_components) :
946 (point_ptrs.back() * n_components));
947 buffer.resize(
948 (sort_data ? buffer_size_with_sorting : buffer_size_without_sorting) *
949 n_components);
950
951 // ... for evaluation
952 ArrayView<T> buffer_eval(buffer.data(),
953 send_permutation.size() * n_components);
954
955 // ... for communication (send)
956 ArrayView<T> buffer_send(sort_data ?
957 (buffer.data() +
958 send_permutation.size() * n_components) :
959 const_cast<T *>(input.data()),
960 point_ptrs.back() * n_components);
961
962 // more arrays
963 std::vector<MPI_Request> send_requests;
964 std::vector<std::vector<char>> send_buffers_packed;
965 std::vector<char> recv_buffer_packed;
966
967 // sort for communication (and duplicate data if necessary)
968
969 if (sort_data)
970 {
971 const auto my_rank_local_recv_ptr =
972 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
973
974 if (my_rank_local_recv_ptr != recv_ranks.end())
975 {
976 // optimize the case if we have our own rank among the possible
977 // list
978 const unsigned int my_rank_local_recv =
979 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
980 const unsigned int my_rank_local_send = std::distance(
981 send_ranks.begin(),
982 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
983
984 const unsigned int start = recv_ptrs[my_rank_local_recv];
985 const unsigned int end = recv_ptrs[my_rank_local_recv + 1];
986 const unsigned int *send_ptr =
987 send_permutation_inv.data() + send_ptrs[my_rank_local_send];
988 for (unsigned int i = 0, k = 0; i < point_ptrs.size() - 1; ++i)
989 {
990 const unsigned int next = point_ptrs[i + 1];
991 for (unsigned int j = point_ptrs[i]; j < next; ++j, ++k)
992 {
993 const unsigned int recv_index = recv_permutation_inv[k];
994
995 // local data -> can be copied to final buffer directly
996 if (start <= recv_index && recv_index < end)
997 for (unsigned int c = 0; c < n_components; ++c)
998 buffer_eval[send_ptr[recv_index - start] *
999 n_components +
1000 c] = input[i * n_components + c];
1001 else // data to be sent
1002 for (unsigned int c = 0; c < n_components; ++c)
1003 buffer_send[recv_index * n_components + c] =
1004 input[i * n_components + c];
1005 }
1006 }
1007 }
1008 else
1009 {
1010 for (unsigned int i = 0, k = 0; i < point_ptrs.size() - 1; ++i)
1011 for (unsigned int j = point_ptrs[i]; j < point_ptrs[i + 1];
1012 ++j, ++k)
1013 for (unsigned int c = 0; c < n_components; ++c)
1014 buffer_send[recv_permutation_inv[k] * n_components + c] =
1015 input[i * n_components + c];
1016 }
1017 }
1018
1019 // send data
1020 send_buffers_packed.reserve(recv_ranks.size());
1021 send_requests.reserve(recv_ranks.size());
1022
1023 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
1024 {
1025 if (recv_ranks[i] == my_rank)
1026 continue;
1027
1028 internal::pack_and_isend(
1029 ArrayView<const T>(
1030 buffer_send.begin() + recv_ptrs[i] * n_components,
1031 (recv_ptrs[i + 1] - recv_ptrs[i]) * n_components),
1032 recv_ranks[i],
1033 internal::Tags::remote_point_evaluation,
1034 tria->get_mpi_communicator(),
1035 send_buffers_packed,
1036 send_requests);
1037 }
1038
1039 // copy local data directly to output if no sorting is requested
1040 if (!sort_data)
1041 {
1042 const auto my_rank_local_recv_ptr =
1043 std::find(recv_ranks.begin(), recv_ranks.end(), my_rank);
1044
1045 if (my_rank_local_recv_ptr != recv_ranks.end())
1046 {
1047 const unsigned int my_rank_local_recv =
1048 std::distance(recv_ranks.begin(), my_rank_local_recv_ptr);
1049 const unsigned int my_rank_local_send = std::distance(
1050 send_ranks.begin(),
1051 std::find(send_ranks.begin(), send_ranks.end(), my_rank));
1052
1053 for (unsigned int j = recv_ptrs[my_rank_local_recv],
1054 k = send_ptrs[my_rank_local_send];
1055 j < recv_ptrs[my_rank_local_recv + 1];
1056 ++j, ++k)
1057 for (unsigned int c = 0; c < n_components; ++c)
1058 buffer_eval[k * n_components + c] =
1059 input[j * n_components + c];
1060 }
1061 }
1062
1063 for (unsigned int i = 0; i < send_ranks.size(); ++i)
1064 {
1065 if (send_ranks[i] == my_rank)
1066 continue;
1067
1068 // receive remote data
1069 MPI_Status status;
1070 int ierr = MPI_Probe(MPI_ANY_SOURCE,
1071 internal::Tags::remote_point_evaluation,
1072 tria->get_mpi_communicator(),
1073 &status);
1074 AssertThrowMPI(ierr);
1075
1076 // write data into buffer vector
1077 const auto ptr =
1078 std::find(send_ranks.begin(), send_ranks.end(), status.MPI_SOURCE);
1079
1080 Assert(ptr != send_ranks.end(), ExcNotImplemented());
1081
1082 const unsigned int j = std::distance(send_ranks.begin(), ptr);
1083
1084 ArrayView<T> recv_buffer(
1085 sort_data ?
1086 (buffer.data() +
1087 (point_ptrs.back() + send_permutation.size()) * n_components) :
1088 (buffer.data() + send_ptrs[j] * n_components),
1089 (send_ptrs[j + 1] - send_ptrs[j]) * n_components);
1090
1091 internal::recv_and_unpack(recv_buffer,
1092 tria->get_mpi_communicator(),
1093 status,
1094 recv_buffer_packed);
1095
1096 if (sort_data)
1097 {
1098 for (unsigned int i = send_ptrs[j], k = 0; i < send_ptrs[j + 1];
1099 ++i, ++k)
1100 for (unsigned int c = 0; c < n_components; ++c)
1101 buffer_eval[send_permutation_inv[i] * n_components + c] =
1102 recv_buffer[k * n_components + c];
1103 }
1104 }
1105
1106 if (!send_requests.empty())
1107 {
1108 const int ierr = MPI_Waitall(send_requests.size(),
1109 send_requests.data(),
1110 MPI_STATUSES_IGNORE);
1111 AssertThrowMPI(ierr);
1112 }
1113
1114 // evaluate function at points
1115 evaluation_function(buffer_eval, *cell_data);
1116#endif
1117 }
1118
1119
1120
1121 template <int dim, int spacedim>
1122 template <typename T, unsigned int n_components>
1123 void
1124 RemotePointEvaluation<dim, spacedim>::process_and_evaluate(
1125 const std::vector<T> &input,
1126 const std::function<void(const ArrayView<const T> &, const CellData &)>
1127 &evaluation_function,
1128 const bool sort_data) const
1129 {
1130 std::vector<T> buffer;
1131 this->process_and_evaluate<T, n_components>(input,
1132 buffer,
1133 evaluation_function,
1134 sort_data);
1135 }
1136
1137 } // end of namespace MPI
1138} // end of namespace Utilities
1139
1140
1141DEAL_II_NAMESPACE_CLOSE
1142
1143#endif
ArrayView::begin
iterator begin() const
Definition array_view.h:702
ArrayView::empty
bool empty() const
Definition array_view.h:693
ArrayView::data
value_type * data() const noexcept
Definition array_view.h:661
ArrayView::size
std::size_t size() const
Definition array_view.h:684
MPI_Comm
Mapping
Abstract base class for mapping classes.
Definition mapping.h:318
Point
Definition point.h:111
Triangulation::get_mpi_communicator
virtual MPI_Comm get_mpi_communicator() const
Utilities::MPI::RemotePointEvaluation::CellData::get_data_view
ArrayView< T > get_data_view(const unsigned int cell, const ArrayView< T > &values) const
Definition mpi_remote_point_evaluation.h:695
Utilities::MPI::RemotePointEvaluation::CellData::cell_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > cell_indices() const
Definition mpi_remote_point_evaluation.cc:279
Utilities::MPI::RemotePointEvaluation::CellData::get_unit_points
ArrayView< const Point< dim > > get_unit_points(const unsigned int cell) const
Definition mpi_remote_point_evaluation.cc:299
Utilities::MPI::RemotePointEvaluation::CellData::get_active_cell_iterator
Triangulation< dim, spacedim >::active_cell_iterator get_active_cell_iterator(const unsigned int cell) const
Definition mpi_remote_point_evaluation.cc:288
Utilities::MPI::RemotePointEvaluation::get_send_permutation
const std::vector< unsigned int > & get_send_permutation() const
Definition mpi_remote_point_evaluation.cc:389
Utilities::MPI::RemotePointEvaluation::mapping
ObserverPointer< const Mapping< dim, spacedim > > mapping
Definition mpi_remote_point_evaluation.h:446
Utilities::MPI::RemotePointEvaluation::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
Definition mpi_remote_point_evaluation.cc:362
Utilities::MPI::RemotePointEvaluation::process_and_evaluate
void process_and_evaluate(const std::vector< T > &input, std::vector< T > &buffer, const std::function< void(const ArrayView< const T > &, const CellData &)> &evaluation_function, const bool sort_data=true) const
Definition mpi_remote_point_evaluation.h:921
Utilities::MPI::RemotePointEvaluation::tria
ObserverPointer< const Triangulation< dim, spacedim > > tria
Definition mpi_remote_point_evaluation.h:441
Utilities::MPI::RemotePointEvaluation::get_point_ptrs
const std::vector< unsigned int > & get_point_ptrs() const
Definition mpi_remote_point_evaluation.cc:320
Utilities::MPI::RemotePointEvaluation::evaluate_and_process
std::vector< T > evaluate_and_process(const std::function< void(const ArrayView< T > &, const CellData &)> &evaluation_function, const bool sort_data=true) const
Definition mpi_remote_point_evaluation.h:900
Utilities::MPI::RemotePointEvaluation::get_inverse_recv_permutation
const std::vector< unsigned int > & get_inverse_recv_permutation() const
Definition mpi_remote_point_evaluation.cc:398
Utilities::MPI::RemotePointEvaluation::evaluate_and_process
void evaluate_and_process(std::vector< T > &output, std::vector< T > &buffer, const std::function< void(const ArrayView< T > &, const CellData &)> &evaluation_function, const bool sort_data=true) const
Definition mpi_remote_point_evaluation.h:709
Utilities::MPI::RemotePointEvaluation::get_mapping
const Mapping< dim, spacedim > & get_mapping() const
Definition mpi_remote_point_evaluation.cc:371
Utilities::MPI::RemotePointEvaluation::process_and_evaluate
void process_and_evaluate(const std::vector< T > &input, const std::function< void(const ArrayView< const T > &, const CellData &)> &evaluation_function, const bool sort_data=true) const
Definition mpi_remote_point_evaluation.h:1124
Utilities::MPI::RemotePointEvaluation::reinit
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
Definition mpi_remote_point_evaluation.cc:86
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_DEPRECATED_WITH_COMMENT
#define DEAL_II_DEPRECATED_WITH_COMMENT(comment)
Definition config.h:206
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
dof_handler.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
StandardExceptions::ExcNeedsMPI
static ::ExceptionBase & ExcNeedsMPI()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exceptions.h:2174
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
InitializeLibrary::MPI
@ MPI
my_rank
const unsigned int my_rank
Definition mpi.cc:918
data
std::vector< index_type > data
Definition mpi.cc:735
comm
const MPI_Comm comm
Definition mpi.cc:913
mpi_tags.h
Utilities::MPI::internal::recv_and_unpack
std::enable_if_t< Utilities::MPI::is_mpi_type< T >==false, void > recv_and_unpack(const ArrayView< T > &data, const MPI_Comm comm, const MPI_Status &status, std::vector< char > &buffer)
Definition mpi_remote_point_evaluation.h:619
Utilities::MPI::internal::pack_and_isend
std::enable_if_t< Utilities::MPI::is_mpi_type< T >==false, void > pack_and_isend(const ArrayView< const T > &data, const unsigned int rank, const unsigned int tag, const MPI_Comm comm, std::vector< std::vector< char > > &buffers, std::vector< MPI_Request > &requests)
Definition mpi_remote_point_evaluation.h:543
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
Utilities::pack
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1381
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std_cxx20::ranges::iota_view
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
std
STL namespace.
Generated by doxygen 1.9.8