Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+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
fe_remote_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 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
16#ifndef dealii_matrix_free_fe_remote_evaluation_h
17#define dealii_matrix_free_fe_remote_evaluation_h
18
20
23
27
29
30#include <algorithm>
31#include <variant>
32
34
35namespace internal
36{
41 template <int dim, int n_components, typename value_type_>
43 {
44 using value_type = typename internal::FEPointEvaluation::
45 EvaluatorTypeTraits<dim, dim, n_components, value_type_>::value_type;
46
47 using gradient_type = typename internal::FEPointEvaluation::
48 EvaluatorTypeTraits<dim, dim, n_components, value_type_>::
49 real_gradient_type;
50
55
60 };
61
71 {
75 unsigned int
76 get_shift(const unsigned int index) const;
77
81 unsigned int
82 get_shift(const unsigned int cell_index,
83 const unsigned int face_number) const;
84
88 unsigned int
89 size() const;
90
94 unsigned int start = 0;
95
99 std::vector<unsigned int> ptrs_ptrs;
100
104 std::vector<unsigned int> ptrs;
105 };
106
112 template <int dim, int n_components, typename value_type_>
114 {
121
122 public:
129
136 const value_type
137 get_value(const unsigned int q) const;
138
145 const gradient_type
146 get_gradient(const unsigned int q) const;
147
156 void
157 reinit(const unsigned int index);
158
165 void
166 reinit(const unsigned int index_0, const unsigned int index_1);
167
168 private:
174
179
183 unsigned int data_offset;
184 };
185
186} // namespace internal
187
188
189
199template <int dim>
201{
205 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe;
206
211 std::vector<std::pair<unsigned int, unsigned int>> batch_id_n_entities;
212
218 std::vector<std::pair<unsigned int, unsigned int>>
220};
221
228template <int dim>
230{
234 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe;
235
239 std::vector<unsigned int> indices;
240
246 std::vector<unsigned int>
248};
249
256template <int dim>
258{
262 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe;
263
268 std::vector<
269 std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
271
277 std::vector<
278 std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
280};
281
285template <int dim>
287{
288public:
294 void
296 &comm_objects,
297 const std::pair<unsigned int, unsigned int> &face_batch_range,
298 const std::vector<unsigned int> &quadrature_sizes);
299
305 void
307 const std::vector<FERemoteCommunicationObject<dim>> &comm_objects,
308 const std::pair<unsigned int, unsigned int> &face_range,
309 const std::vector<unsigned int> &quadrature_sizes);
310
317 template <typename Iterator>
318 void
320 const std::vector<FERemoteCommunicationObjectTwoLevel<dim>> &comm_objects,
321 const IteratorRange<Iterator> &cell_iterator_range,
322 const std::vector<std::vector<unsigned int>> &quadrature_sizes);
323
324
325
329 template <int n_components,
330 typename PrecomputedEvaluationDataType,
331 typename MeshType,
332 typename VectorType>
333 void
335 PrecomputedEvaluationDataType &dst,
336 const MeshType &mesh,
337 const VectorType &src,
338 const EvaluationFlags::EvaluationFlags eval_flags,
339 const unsigned int first_selected_component,
341
346 get_view() const;
347
348private:
354
358 std::vector<std::variant<FERemoteCommunicationObjectEntityBatches<dim>,
362
368 {
369 public:
374 template <typename T1, typename T2>
375 static void
378 const std::vector<T2> &src,
379 const std::vector<unsigned int> &indices);
380
385 template <typename T1, typename T2>
386 static void
387 copy_data(
390 const std::vector<T2> &src,
391 const std::vector<std::pair<typename Triangulation<dim>::cell_iterator,
392 unsigned int>> &cell_face_nos);
393
398 template <typename T1, typename T2>
399 static void
402 const std::vector<T2> &src,
403 const std::vector<std::pair<unsigned int, unsigned int>>
404 &batch_id_n_entities);
405
406 private:
410 template <typename T1, std::size_t n_lanes>
411 static void
413 const unsigned int v,
414 const T1 &src);
415
419 template <typename T1, int rank_, std::size_t n_lanes, int dim_>
420 static void
422 const unsigned int v,
423 const Tensor<rank_, dim_, T1> &src);
424
428 template <typename T1,
429 int rank_,
430 std::size_t n_lanes,
431 int n_components_,
432 int dim_>
433 static void
435 Tensor<rank_,
436 n_components_,
437 Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>>> &dst,
438 const unsigned int v,
439 const Tensor<rank_, n_components_, Tensor<rank_, dim_, T1>> &src);
440
445 template <typename T1, typename T2>
446 static void
447 copy_data_entries(T1 &, const unsigned int, const T2 &);
448 };
449};
450
451
452
457namespace Utilities
458{
476 template <int dim,
477 typename Number,
478 typename VectorizedArrayType = VectorizedArray<Number>>
482 const std::vector<
483 std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
484 &non_matching_faces_marked_vertices,
485 const unsigned int quad_no = 0,
486 const unsigned int dof_no = 0,
487 const double tolerance = 1e-9);
488
489
490
509 template <int dim,
510 typename Number,
511 typename VectorizedArrayType = VectorizedArray<Number>>
515 const std::vector<
516 std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
517 &non_matching_faces_marked_vertices,
518 const unsigned int n_q_pnts_1D,
519 const unsigned int dof_no = 0,
520 NonMatching::MappingInfo<dim, dim, Number> *nm_mapping_info = nullptr,
521 const double tolerance = 1e-9);
522} // namespace Utilities
523
524
525
537template <int dim, int n_components, typename value_type>
634
635
636
637namespace internal
638{
639 unsigned int
641 {
642 Assert(ptrs_ptrs.size() == 0, ExcMessage("Two level CRS set up"));
643
645 ExcMessage("Index has to be valid!"));
646
648 AssertIndexRange(index - start, ptrs.size());
649 return ptrs[index - start];
650 }
651
652 unsigned int
654 const unsigned int face_number) const
655 {
656 Assert(ptrs_ptrs.size() > 0, ExcMessage("No two level CRS set up"));
657
659 ExcMessage("Cell index has to be valid!"));
661 ExcMessage("Face number has to be valid!"));
662
664
666 const unsigned int face_index = ptrs_ptrs[cell_index - start] + face_number;
667 AssertIndexRange(face_index, ptrs.size());
668 return ptrs[face_index];
669 }
670
671 unsigned int
673 {
674 Assert(ptrs.size() > 0, ExcInternalError());
675 return ptrs.back();
676 }
677
678 template <int dim, int n_components, typename value_type_>
687
688 template <int dim, int n_components, typename value_type_>
689 const typename PrecomputedEvaluationData<dim,
690 n_components,
691 value_type_>::value_type
693 const unsigned int q) const
694 {
696 ExcMessage("reinit() not called."));
697 AssertIndexRange(data_offset + q, data.values.size());
698 return data.values[data_offset + q];
699 }
700
701 template <int dim, int n_components, typename value_type_>
703 gradient_type
705 get_gradient(const unsigned int q) const
706 {
708 ExcMessage("reinit() not called."));
709 AssertIndexRange(data_offset + q, data.gradients.size());
710 return data.gradients[data_offset + q];
711 }
712
713 template <int dim, int n_components, typename value_type_>
714 void
716 const unsigned int index)
717 {
718 data_offset = view.get_shift(index);
719 }
720
721 template <int dim, int n_components, typename value_type_>
722 void
724 const unsigned int index_0,
725 const unsigned int index_1)
726 {
727 data_offset = view.get_shift(index_0, index_1);
728 }
729
730} // namespace internal
731
732
733
734template <int dim>
735std::vector<std::pair<unsigned int, unsigned int>>
737 const
738{
739 return batch_id_n_entities;
740}
741
742template <int dim>
743std::vector<unsigned int>
748
749template <int dim>
750std::vector<std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
755
756
757
758template <int dim>
759void
762 &comm_objects,
763 const std::pair<unsigned int, unsigned int> &face_batch_range,
764 const std::vector<unsigned int> &quadrature_sizes)
765{
766 // erase type by converting to the base object
767 communication_objects.clear();
768 for (const auto &co : comm_objects)
769 communication_objects.push_back(co);
770
771 // fetch points and update communication patterns
772 const unsigned int n_cells = quadrature_sizes.size();
773 AssertDimension(n_cells, face_batch_range.second - face_batch_range.first);
774
775 // construct view:
776 view.start = face_batch_range.first;
777
778 view.ptrs.resize(n_cells + 1);
779
780 view.ptrs[0] = 0;
781 for (unsigned int face = 0; face < n_cells; ++face)
782 {
783 view.ptrs[face + 1] = view.ptrs[face] + quadrature_sizes[face];
784 }
785}
786
787template <int dim>
788void
790 const std::vector<FERemoteCommunicationObject<dim>> &comm_objects,
791 const std::pair<unsigned int, unsigned int> &face_range,
792 const std::vector<unsigned int> &quadrature_sizes)
793{
794 // erase type
795 communication_objects.clear();
796 for (const auto &co : comm_objects)
797 communication_objects.push_back(co);
798
799 const unsigned int n_faces = quadrature_sizes.size();
800 AssertDimension(n_faces, face_range.second - face_range.first);
801
802 // construct view:
803 view.start = face_range.first;
804
805 view.ptrs.resize(n_faces + 1);
806
807 view.ptrs[0] = 0;
808 for (unsigned int face = 0; face < n_faces; ++face)
809 view.ptrs[face + 1] = view.ptrs[face] + quadrature_sizes[face];
810}
811
812template <int dim>
813template <typename Iterator>
814void
816 const std::vector<FERemoteCommunicationObjectTwoLevel<dim>> &comm_objects,
817 const IteratorRange<Iterator> &cell_iterator_range,
818 const std::vector<std::vector<unsigned int>> &quadrature_sizes)
819{
820 // erase type
821 communication_objects.clear();
822 for (const auto &co : comm_objects)
823 communication_objects.push_back(co);
824
825 const unsigned int n_cells = quadrature_sizes.size();
826 AssertDimension(n_cells,
827 std::distance(cell_iterator_range.begin(),
828 cell_iterator_range.end()));
829
830 // construct view:
831 auto &cell_ptrs = view.ptrs_ptrs;
832 auto &face_ptrs = view.ptrs;
833
834 view.start = 0;
835 cell_ptrs.resize(n_cells);
836 unsigned int n_faces = 0;
837 for (const auto &cell : cell_iterator_range)
838 {
839 cell_ptrs[cell->active_cell_index()] = n_faces;
840 n_faces += cell->n_faces();
841 }
842
843 face_ptrs.resize(n_faces + 1);
844 face_ptrs[0] = 0;
845 for (const auto &cell : cell_iterator_range)
846 {
847 for (const auto &f : cell->face_indices())
848 {
849 const unsigned int face_index =
850 cell_ptrs[cell->active_cell_index()] + f;
851
852 face_ptrs[face_index + 1] =
853 face_ptrs[face_index] +
854 quadrature_sizes[cell->active_cell_index()][f];
855 }
856 }
857}
858
859template <int dim>
860template <int n_components,
861 typename PrecomputedEvaluationDataType,
862 typename MeshType,
863 typename VectorType>
864void
866 PrecomputedEvaluationDataType &dst,
867 const MeshType &mesh,
868 const VectorType &src,
869 const EvaluationFlags::EvaluationFlags eval_flags,
870 const unsigned int first_selected_component,
872{
873 const bool has_ghost_elements = src.has_ghost_elements();
874
875 if (has_ghost_elements == false)
876 src.update_ghost_values();
877
878
879 for (const auto &communication_object : communication_objects)
880 {
881 if (eval_flags & EvaluationFlags::values)
882 {
883 std::visit(
884 [&](const auto &obj) {
885 CopyInstructions::copy_data(
886 view,
887 dst.values,
888 VectorTools::point_values<n_components>(
889 *obj.rpe, mesh, src, vec_flags, first_selected_component),
890 obj.get_communication_object_pntrs());
891 },
892 communication_object);
893 }
894
895 if (eval_flags & EvaluationFlags::gradients)
896 {
897 std::visit(
898 [&](const auto &obj) {
899 CopyInstructions::copy_data(
900 view,
901 dst.gradients,
902 VectorTools::point_gradients<n_components>(
903 *obj.rpe, mesh, src, vec_flags, first_selected_component),
904 obj.get_communication_object_pntrs());
905 },
906 communication_object);
907 }
908
910 }
911
912 if (has_ghost_elements == false)
913 src.zero_out_ghost_values();
914}
915
916template <int dim>
919{
920 return view;
921}
922
923template <int dim>
924template <typename T1, typename T2>
925void
929 const std::vector<T2> &src,
930 const std::vector<unsigned int> &indices)
931{
932 dst.resize(view.size());
933
934 unsigned int c = 0;
935 for (const auto idx : indices)
936 {
937 for (unsigned int j = view.get_shift(idx); j < view.get_shift(idx + 1);
938 ++j, ++c)
939 {
940 AssertIndexRange(j, dst.size());
941 AssertIndexRange(c, src.size());
942 dst[j] = src[c];
943 }
944 }
945}
946
947template <int dim>
948template <typename T1, typename T2>
949void
953 const std::vector<T2> &src,
954 const std::vector<std::pair<typename Triangulation<dim>::cell_iterator,
955 unsigned int>> &cell_face_nos)
956{
957 dst.resize(view.size());
958
959 unsigned int c = 0;
960 for (const auto &[cell, f] : cell_face_nos)
961 {
962 for (unsigned int j = view.get_shift(cell->active_cell_index(), f);
963 j < view.get_shift(cell->active_cell_index(), f + 1);
964 ++j, ++c)
965 {
966 AssertIndexRange(j, dst.size());
967 AssertIndexRange(c, src.size());
968
969 dst[j] = src[c];
970 }
971 }
972}
973
974template <int dim>
975template <typename T1, typename T2>
976void
980 const std::vector<T2> &src,
981 const std::vector<std::pair<unsigned int, unsigned int>> &batch_id_n_entities)
982{
983 dst.resize(view.size());
984
985 unsigned int c = 0;
986 for (const auto &[batch_id, n_entries] : batch_id_n_entities)
987 {
988 for (unsigned int v = 0; v < n_entries; ++v)
989 for (unsigned int j = view.get_shift(batch_id);
990 j < view.get_shift(batch_id + 1);
991 ++j, ++c)
992 {
993 AssertIndexRange(j, dst.size());
994 AssertIndexRange(c, src.size());
995
996 copy_data_entries(dst[j], v, src[c]);
997 }
998 }
999}
1000
1001template <int dim>
1002template <typename T1, std::size_t n_lanes>
1003void
1006 const unsigned int v,
1007 const T1 &src)
1008{
1009 AssertIndexRange(v, n_lanes);
1010
1011 dst[v] = src;
1012}
1013
1014template <int dim>
1015template <typename T1, int rank_, std::size_t n_lanes, int dim_>
1016void
1018 Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>> &dst,
1019 const unsigned int v,
1020 const Tensor<rank_, dim_, T1> &src)
1021{
1022 AssertIndexRange(v, n_lanes);
1023
1024 if constexpr (rank_ == 1)
1025 {
1026 for (unsigned int i = 0; i < dim_; ++i)
1027 dst[i][v] = src[i];
1028 }
1029 else
1030 {
1031 for (unsigned int i = 0; i < rank_; ++i)
1032 for (unsigned int j = 0; j < dim_; ++j)
1033 dst[i][j][v] = src[i][j];
1034 }
1035}
1036
1037template <int dim>
1038template <typename T1,
1039 int rank_,
1040 std::size_t n_lanes,
1041 int n_components_,
1042 int dim_>
1043void
1045 Tensor<rank_,
1046 n_components_,
1047 Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>>> &dst,
1048 const unsigned int v,
1049 const Tensor<rank_, n_components_, Tensor<rank_, dim_, T1>> &src)
1050{
1051 if constexpr (rank_ == 1)
1052 {
1053 for (unsigned int i = 0; i < n_components_; ++i)
1054 copy_data(dst[i], v, src[i]);
1055 }
1056 else
1057 {
1058 for (unsigned int i = 0; i < rank_; ++i)
1059 for (unsigned int j = 0; j < n_components_; ++j)
1060 dst[i][j][v] = src[i][j];
1061 }
1062}
1063
1064template <int dim>
1065template <typename T1, typename T2>
1066void
1068 T1 &,
1069 const unsigned int,
1070 const T2 &)
1071{
1072 Assert(false,
1073 ExcMessage(
1074 "copy_data_entries() not implemented for given arguments."));
1075}
1076
1077
1078
1079namespace Utilities
1080{
1081 template <int dim, typename Number, typename VectorizedArrayType>
1085 const std::vector<
1086 std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
1087 &non_matching_faces_marked_vertices,
1088 const unsigned int quad_no,
1089 const unsigned int dof_no,
1090 const double tolerance)
1091 {
1092 const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
1093 const auto &tria = dof_handler.get_triangulation();
1094 const auto &mapping = *matrix_free.get_mapping_info().mapping;
1095
1096 // Communication objects know about the communication pattern. I.e.,
1097 // they know about the cells and quadrature points that have to be
1098 // evaluated at remote faces. This information is given via
1099 // RemotePointEvaluation. Additionally, the communication objects
1100 // have to be able to match the quadrature points of the remote
1101 // points (that provide exterior information) to the quadrature points
1102 // defined at the interior cell. In case of point-to-point interpolation
1103 // a vector of pairs with face batch Ids and the number of faces in the
1104 // batch is needed. @c FERemoteCommunicationObjectEntityBatches
1105 // is a container to store this information.
1106 //
1107 // We need multiple communication objects (one for each non-matching face
1108 // ID).
1109 std::vector<FERemoteCommunicationObjectEntityBatches<dim>> comm_objects;
1110
1111 // Additionally to the communication objects we need a vector
1112 // that stores quadrature rule sizes for every face batch.
1113 // The quadrature can have size zero in case of non non-matching faces,
1114 // i.e. boundary faces. Internally this information is needed to correctly
1115 // access values over multiple communication objects.
1116 std::vector<unsigned int> global_quadrature_sizes(
1118
1119 // Get the range of face batches we have to look at during construction of
1120 // the communication objects. We only have to look at boundary faces.
1121 const auto face_batch_range =
1122 std::make_pair(matrix_free.n_inner_face_batches(),
1123 matrix_free.n_inner_face_batches() +
1124 matrix_free.n_boundary_face_batches());
1125
1126 // Iterate over all non-matching face IDs.
1127 for (const auto &[nm_face, marked_vertices] :
1128 non_matching_faces_marked_vertices)
1129 {
1130 // Construct the communication object for every face ID:
1131 // 1) RemotePointEvaluation with user specified function for marked
1132 // vertices.
1133 auto rpe = std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(
1134 tolerance, false, 0, marked_vertices);
1135
1136 // 2) Face batch IDs and number of faces in batch.
1137 std::vector<std::pair<unsigned int, unsigned int>>
1138 face_batch_id_n_faces;
1139
1140 // Points that are searched by rpe.
1141 std::vector<Point<dim>> points;
1142
1143 // Temporarily set up FEFaceEvaluation to access the quadrature points
1144 // at the faces on the non-matching interface.
1145 FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free,
1146 true,
1147 dof_no,
1148 quad_no);
1149
1150 // Iterate over the boundary faces.
1151 for (unsigned int bface = 0;
1152 bface < face_batch_range.second - face_batch_range.first;
1153 ++bface)
1154 {
1155 const unsigned int face = face_batch_range.first + bface;
1156
1157 if (matrix_free.get_boundary_id(face) == nm_face)
1158 {
1159 phi.reinit(face);
1160
1161 // If @c face is on the current side of the non-matching
1162 // interface. Add the face batch ID and the number of faces in
1163 // the batch to the corresponding data structure.
1164 const unsigned int n_faces =
1165 matrix_free.n_active_entries_per_face_batch(face);
1166 face_batch_id_n_faces.emplace_back(face, n_faces);
1167
1168 // Append the quadrature points to the points we need to search
1169 // for.
1170 for (unsigned int v = 0; v < n_faces; ++v)
1171 {
1172 for (unsigned int q : phi.quadrature_point_indices())
1173 {
1174 const auto point = phi.quadrature_point(q);
1175 Point<dim> temp;
1176 for (unsigned int i = 0; i < dim; ++i)
1177 temp[i] = point[i][v];
1178
1179 points.push_back(temp);
1180 }
1181 }
1182
1183 // Insert the quadrature size into the global vector.
1184 // First check that each face is only considered once.
1185 Assert(global_quadrature_sizes[bface] ==
1187 ExcMessage(
1188 "Quadrature for given face already provided."));
1189
1190 global_quadrature_sizes[bface] = phi.n_q_points;
1191 }
1192 }
1193
1194 // Reinit RPE and ensure all points are found.
1195 rpe->reinit(points, tria, mapping);
1196 Assert(rpe->all_points_found(),
1197 ExcMessage("Not all remote points found."));
1198
1199 // Add communication object to the list of objects.
1201 co.batch_id_n_entities = face_batch_id_n_faces;
1202 co.rpe = rpe;
1203 comm_objects.push_back(co);
1204 }
1205
1206 // Reinit the communicator `FERemoteEvaluationCommunicator`
1207 // with the communication objects.
1208 FERemoteEvaluationCommunicator<dim> remote_communicator;
1209
1210 // if no quadrature size is set, an empty quadrature is considered
1211 std::replace(global_quadrature_sizes.begin(),
1212 global_quadrature_sizes.end(),
1214 0u);
1215
1216 remote_communicator.reinit_faces(comm_objects,
1217 face_batch_range,
1218 global_quadrature_sizes);
1219
1220 return remote_communicator;
1221 }
1222
1223
1224
1225 template <int dim, typename Number, typename VectorizedArrayType>
1229 const std::vector<
1230 std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
1231 &non_matching_faces_marked_vertices,
1232 const unsigned int n_q_pnts_1D,
1233 const unsigned int dof_no,
1235 const double tolerance)
1236 {
1237 const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
1238 const auto &tria = dof_handler.get_triangulation();
1239 const auto &mapping = *matrix_free.get_mapping_info().mapping;
1240
1241 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
1242
1243 std::pair<unsigned int, unsigned int> face_range =
1244 std::make_pair(matrix_free.n_inner_face_batches(),
1245 matrix_free.n_inner_face_batches() +
1246 matrix_free.n_boundary_face_batches());
1247
1248 std::vector<Quadrature<dim - 1>> global_quadrature_vector(
1249 (matrix_free.n_inner_face_batches() +
1250 matrix_free.n_boundary_face_batches()) *
1251 n_lanes);
1252
1253 // In case of Nitsche-type mortaring a vector of face indices is
1254 // needed as communication object.
1255 // @c FERemoteCommunicationObjectFaces is a container to store this
1256 // information.
1257 //
1258 // We need multiple communication objects (one for each non-matching face
1259 // ID).
1260 std::vector<FERemoteCommunicationObject<dim>> comm_objects;
1261
1262 // Create bounding boxes and GridTools::Cache which is needed in
1263 // the following loop.
1264 std::vector<BoundingBox<dim>> local_boxes;
1265 for (const auto &cell : tria.active_cell_iterators())
1266 if (cell->is_locally_owned())
1267 local_boxes.emplace_back(mapping.get_bounding_box(cell));
1268
1269 // Create r-tree of bounding boxes
1270 const auto local_tree = pack_rtree(local_boxes);
1271
1272 // Compress r-tree to a minimal set of bounding boxes
1273 std::vector<std::vector<BoundingBox<dim>>> global_bboxes(1);
1274 global_bboxes[0] = extract_rtree_level(local_tree, 0);
1275
1276 const GridTools::Cache<dim, dim> cache(tria, mapping);
1277
1278 // Iterate over all sides of the non-matching interface.
1279 for (const auto &[nm_face, marked_vertices] :
1280 non_matching_faces_marked_vertices)
1281 {
1282 // 1) compute cell face pairs
1283 std::vector<
1284 std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
1285 cell_face_pairs;
1286
1287 std::vector<unsigned int> indices;
1288
1289 for (unsigned int face = face_range.first; face < face_range.second;
1290 ++face)
1291 {
1292 if (matrix_free.get_boundary_id(face) == nm_face)
1293 {
1294 for (unsigned int v = 0;
1295 v < matrix_free.n_active_entries_per_face_batch(face);
1296 ++v)
1297 {
1298 const auto &[c, f] = matrix_free.get_face_iterator(face, v);
1299
1300 cell_face_pairs.emplace_back(std::make_pair(c, f));
1301 indices.push_back(face * n_lanes + v);
1302 }
1303 }
1304 }
1305
1306 // 2) Create RPE.
1307 // In the Nitsche-type case we do not collect points for the setup
1308 // of RemotePointEvaluation. Instead we compute intersections between
1309 // the faces and set up RemotePointEvaluation with the computed
1310 // intersections.
1311
1312 // Build intersection requests. Intersection requests
1313 // correspond to vertices at faces.
1314 std::vector<std::vector<Point<dim>>> intersection_requests;
1315 for (const auto &[cell, f] : cell_face_pairs)
1316 {
1317 std::vector<Point<dim>> vertices(cell->face(f)->n_vertices());
1318 std::copy_n(mapping.get_vertices(cell, f).begin(),
1319 cell->face(f)->n_vertices(),
1320 vertices.begin());
1321 intersection_requests.emplace_back(vertices);
1322 }
1323
1324 // Compute intersection data with user specified function for marked
1325 // vertices.
1326 auto intersection_data =
1328 dim - 1>(cache,
1329 intersection_requests,
1330 global_bboxes,
1331 marked_vertices(),
1332 tolerance);
1333
1334 // Convert to RPE.
1335 std::vector<Quadrature<dim>> mapped_quadratures_recv_comp;
1336
1337 auto rpe =
1338 std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>();
1339 rpe->reinit(
1340 intersection_data
1341 .template convert_to_distributed_compute_point_locations_internal<
1342 dim>(n_q_pnts_1D, tria, mapping, &mapped_quadratures_recv_comp),
1343 tria,
1344 mapping);
1345
1346 // 3) Fill global quadrature vector.
1347 for (unsigned int i = 0; i < intersection_requests.size(); ++i)
1348 {
1349 const auto idx = indices[i];
1350
1351 // We do not use a structural binding here, since with
1352 // C++17 capturing structural bindings in lambdas leads
1353 // to an ill formed program.
1354 const auto &cell = std::get<0>(cell_face_pairs[i]);
1355 const auto &f = std::get<1>(cell_face_pairs[i]);
1356
1357 std::vector<Point<dim - 1>> q_points;
1358 std::vector<double> weights;
1359 for (unsigned int ptr = intersection_data.recv_ptrs[i];
1360 ptr < intersection_data.recv_ptrs[i + 1];
1361 ++ptr)
1362 {
1363 const auto &quad = mapped_quadratures_recv_comp[ptr];
1364
1365 const auto &ps = quad.get_points();
1366 std::transform(
1367 ps.begin(),
1368 ps.end(),
1369 std::back_inserter(q_points),
1370 [&](const Point<dim> &p) {
1371 return mapping.project_real_point_to_unit_point_on_face(
1372 cell, f, p);
1373 });
1374
1375 const auto &ws = quad.get_weights();
1376 weights.insert(weights.end(), ws.begin(), ws.end());
1377 }
1378 Quadrature<dim - 1> quad(q_points, weights);
1379
1380 Assert(global_quadrature_vector[idx].size() == 0,
1381 ExcMessage("Quadrature for given face already provided."));
1382
1383 global_quadrature_vector[idx] = quad;
1384 }
1385
1386 // Add communication object.
1388 co.indices = indices;
1389 co.rpe = rpe;
1390 comm_objects.push_back(co);
1391 }
1392
1393 // Reinit the communicator with the communication objects.
1394 FERemoteEvaluationCommunicator<dim> remote_communicator;
1395
1396 std::vector<unsigned int> global_quadrature_sizes(
1397 global_quadrature_vector.size());
1398 std::transform(global_quadrature_vector.cbegin(),
1399 global_quadrature_vector.cend(),
1400 global_quadrature_sizes.begin(),
1401 [](const auto &q) { return q.size(); });
1402
1403 remote_communicator.reinit_faces(
1404 comm_objects,
1405 std::make_pair(0, global_quadrature_vector.size()),
1406 global_quadrature_sizes);
1407
1408 if (nm_mapping_info != nullptr)
1409 {
1410 std::vector<
1411 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
1412 vector_face_accessors;
1413 vector_face_accessors.reserve((matrix_free.n_inner_face_batches() +
1414 matrix_free.n_boundary_face_batches()) *
1415 n_lanes);
1416
1417 // fill container for inner face batches
1418 unsigned int face_batch = 0;
1419 for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch)
1420 {
1421 for (unsigned int v = 0; v < n_lanes; ++v)
1422 {
1423 if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
1424 vector_face_accessors.push_back(
1425 matrix_free.get_face_iterator(face_batch, v));
1426 else
1427 vector_face_accessors.push_back(
1428 matrix_free.get_face_iterator(face_batch, 0));
1429 }
1430 }
1431 // and boundary face batches
1432 for (; face_batch < (matrix_free.n_inner_face_batches() +
1433 matrix_free.n_boundary_face_batches());
1434 ++face_batch)
1435 {
1436 for (unsigned int v = 0; v < n_lanes; ++v)
1437 {
1438 if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
1439 vector_face_accessors.push_back(
1440 matrix_free.get_face_iterator(face_batch, v));
1441 else
1442 vector_face_accessors.push_back(
1443 matrix_free.get_face_iterator(face_batch, 0));
1444 }
1445 }
1446
1447 nm_mapping_info->reinit_faces(vector_face_accessors,
1448 global_quadrature_vector);
1449 }
1450
1451 return remote_communicator;
1452 }
1453} // namespace Utilities
1454
1455
1456
1457template <int dim, int n_components, typename value_type>
1458template <typename MeshType>
1461 const MeshType &mesh,
1462 const unsigned int first_selected_component,
1463 const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags)
1464 : comm(&comm)
1465 , first_selected_component(first_selected_component)
1466 , evaluation_flags(evaluation_flags)
1467{
1468 set_mesh(mesh);
1469}
1470
1471template <int dim, int n_components, typename value_type>
1472template <typename VectorType>
1473void
1475 const VectorType &src,
1477{
1478 if (tria)
1479 {
1480 Assert(n_components == 1, ExcNotImplemented());
1481 comm->template update_ghost_values<n_components>(this->data,
1482 *tria,
1483 src,
1484 flags,
1485 first_selected_component,
1486 evaluation_flags);
1487 }
1488 else if (dof_handler)
1489 {
1490 comm->template update_ghost_values<n_components>(this->data,
1491 *dof_handler,
1492 src,
1493 flags,
1494 first_selected_component,
1495 evaluation_flags);
1496 }
1497 else
1499}
1500
1501template <int dim, int n_components, typename value_type>
1504{
1506 comm->get_view());
1507 return data_accessor;
1508}
1509
1510template <int dim, int n_components, typename value_type>
1511void
1513 const Triangulation<dim> &tria)
1514{
1515 this->tria = &tria;
1516}
1517
1518template <int dim, int n_components, typename value_type>
1519void
1521 const DoFHandler<dim> &dof_handler)
1522{
1523 this->dof_handler = &dof_handler;
1524}
1525
1527
1528#endif
size_type size() const
void resize(const size_type new_size)
static void copy_data(const internal::PrecomputedEvaluationDataView &view, AlignedVector< T1 > &dst, const std::vector< T2 > &src, const std::vector< unsigned int > &indices)
static void copy_data_entries(VectorizedArray< T1, n_lanes > &dst, const unsigned int v, const T1 &src)
const internal::PrecomputedEvaluationDataView & get_view() const
void reinit_faces(const std::vector< FERemoteCommunicationObjectEntityBatches< dim > > &comm_objects, const std::pair< unsigned int, unsigned int > &face_batch_range, const std::vector< unsigned int > &quadrature_sizes)
void update_ghost_values(PrecomputedEvaluationDataType &dst, const MeshType &mesh, const VectorType &src, const EvaluationFlags::EvaluationFlags eval_flags, const unsigned int first_selected_component, const VectorTools::EvaluationFlags::EvaluationFlags vec_flags) const
internal::PrecomputedEvaluationDataView view
std::vector< std::variant< FERemoteCommunicationObjectEntityBatches< dim >, FERemoteCommunicationObject< dim >, FERemoteCommunicationObjectTwoLevel< dim > > > communication_objects
SmartPointer< const DoFHandler< dim > > dof_handler
FERemoteEvaluation(const FERemoteEvaluationCommunicator< dim > &comm, const MeshType &mesh, const unsigned int first_selected_component=0, const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags=VectorTools::EvaluationFlags::avg)
SmartPointer< const Triangulation< dim > > tria
void set_mesh(const Triangulation< dim > &tria)
SmartPointer< const FERemoteEvaluationCommunicator< dim > > comm
internal::PrecomputedEvaluationData< dim, n_components, value_type > data
const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags
const unsigned int first_selected_component
internal::PrecomputedEvaluationDataAccessor< dim, n_components, value_type > get_data_accessor() const
void gather_evaluate(const VectorType &src, const EvaluationFlags::EvaluationFlags flags)
IteratorOverIterators end() const
IteratorOverIterators begin()
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
unsigned int n_inner_face_batches() const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
unsigned int n_boundary_face_batches() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
void reinit_faces(const ContainerType &cell_iterator_range, const std::vector< std::vector< Quadrature< dim - 1 > > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
Definition point.h:111
const gradient_type get_gradient(const unsigned int q) const
PrecomputedEvaluationDataAccessor(const PrecomputedEvaluationData< dim, n_components, value_type_ > &data, const PrecomputedEvaluationDataView &view)
const PrecomputedEvaluationDataView & view
typename PrecomputedEvaluationData< dim, n_components, value_type_ >::value_type value_type
const PrecomputedEvaluationData< dim, n_components, value_type_ > & data
const value_type get_value(const unsigned int q) const
typename PrecomputedEvaluationData< dim, n_components, value_type_ >::gradient_type gradient_type
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
EvaluationFlags
The EvaluationFlags enum.
DistributedComputeIntersectionLocationsInternal< structdim, spacedim > distributed_compute_intersection_locations(const Cache< dim, spacedim > &cache, const std::vector< std::vector< Point< spacedim > > > &intersection_requests, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance)
FERemoteEvaluationCommunicator< dim > compute_remote_communicator_faces_nitsche_type_mortaring(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::vector< std::pair< types::boundary_id, std::function< std::vector< bool >()> > > &non_matching_faces_marked_vertices, const unsigned int n_q_pnts_1D, const unsigned int dof_no=0, NonMatching::MappingInfo< dim, dim, Number > *nm_mapping_info=nullptr, const double tolerance=1e-9)
FERemoteEvaluationCommunicator< dim > compute_remote_communicator_faces_point_to_point_interpolation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::vector< std::pair< types::boundary_id, std::function< std::vector< bool >()> > > &non_matching_faces_marked_vertices, const unsigned int quad_no=0, const unsigned int dof_no=0, const double tolerance=1e-9)
static const unsigned int invalid_unsigned_int
Definition types.h:220
*braid_SplitCommworld & comm
std::vector< BoundingBox< boost::geometry::dimension< typename Rtree::indexable_type >::value > > extract_rtree_level(const Rtree &tree, const unsigned int level)
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
std::vector< std::pair< unsigned int, unsigned int > > batch_id_n_entities
std::vector< std::pair< unsigned int, unsigned int > > get_communication_object_pntrs() const
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::vector< std::pair< typename Triangulation< dim >::cell_iterator, unsigned int > > cell_face_nos
std::vector< std::pair< typename Triangulation< dim >::cell_iterator, unsigned int > > get_communication_object_pntrs() const
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::vector< unsigned int > indices
std::vector< unsigned int > get_communication_object_pntrs() const
unsigned int get_shift(const unsigned int index) const
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, dim, n_components, value_type_ >::value_type value_type
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, dim, n_components, value_type_ >::real_gradient_type gradient_type
AlignedVector< gradient_type > gradients