Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00: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<Quadrature<dim>> &quadrature_vector);
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<Quadrature<dim - 1>> &quadrature_vector);
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<Quadrature<dim - 1>>> &quadrature_vector);
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<Quadrature<dim>> &quadrature_vector)
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_vector.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_vector[face].size();
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<Quadrature<dim - 1>> &quadrature_vector)
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_vector.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_vector[face].size();
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<Quadrature<dim - 1>>> &quadrature_vector)
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_vector.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_vector[cell->active_cell_index()][f].size();
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 rules for every face batch.
1113 // The quadrature can be empty 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<Quadrature<dim>> global_quadrature_vector(
1117 matrix_free.n_boundary_face_batches());
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(
1167 std::make_pair(face, n_faces));
1168
1169 // Append the quadrature points to the points we need to search
1170 // for.
1171 for (unsigned int v = 0; v < n_faces; ++v)
1172 {
1173 for (unsigned int q : phi.quadrature_point_indices())
1174 {
1175 const auto point = phi.quadrature_point(q);
1176 Point<dim> temp;
1177 for (unsigned int i = 0; i < dim; ++i)
1178 temp[i] = point[i][v];
1179
1180 points.push_back(temp);
1181 }
1182 }
1183
1184 // Insert a quadrature rule of correct size into the global
1185 // quadrature vector. First check that each face is only
1186 // considered once.
1187 Assert(global_quadrature_vector[bface].size() == 0,
1188 ExcMessage(
1189 "Quadrature for given face already provided."));
1190
1191 global_quadrature_vector[bface] =
1192 Quadrature<dim>(phi.n_q_points);
1193 }
1194 }
1195
1196 // Reinit RPE and ensure all points are found.
1197 rpe->reinit(points, tria, mapping);
1198 Assert(rpe->all_points_found(),
1199 ExcMessage("Not all remote points found."));
1200
1201 // Add communication object to the list of objects.
1203 co.batch_id_n_entities = face_batch_id_n_faces;
1204 co.rpe = rpe;
1205 comm_objects.push_back(co);
1206 }
1207
1208 // Reinit the communicator `FERemoteEvaluationCommunicator`
1209 // with the communication objects.
1210 FERemoteEvaluationCommunicator<dim> remote_communicator;
1211
1212 remote_communicator.reinit_faces(comm_objects,
1213 face_batch_range,
1214 global_quadrature_vector);
1215
1216 return remote_communicator;
1217 }
1218
1219
1220
1221 template <int dim, typename Number, typename VectorizedArrayType>
1225 const std::vector<
1226 std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
1227 &non_matching_faces_marked_vertices,
1228 const unsigned int n_q_pnts_1D,
1229 const unsigned int dof_no,
1231 const double tolerance)
1232 {
1233 const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
1234 const auto &tria = dof_handler.get_triangulation();
1235 const auto &mapping = *matrix_free.get_mapping_info().mapping;
1236
1237 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
1238
1239 std::pair<unsigned int, unsigned int> face_range =
1240 std::make_pair(matrix_free.n_inner_face_batches(),
1241 matrix_free.n_inner_face_batches() +
1242 matrix_free.n_boundary_face_batches());
1243
1244 std::vector<Quadrature<dim - 1>> global_quadrature_vector(
1245 (matrix_free.n_inner_face_batches() +
1246 matrix_free.n_boundary_face_batches()) *
1247 n_lanes);
1248
1249 // In case of Nitsche-type mortaring a vector of face indices is
1250 // needed as communication object.
1251 // @c FERemoteCommunicationObjectFaces is a container to store this
1252 // information.
1253 //
1254 // We need multiple communication objects (one for each non-matching face
1255 // ID).
1256 std::vector<FERemoteCommunicationObject<dim>> comm_objects;
1257
1258 // Create bounding boxes and GridTools::Cache which is needed in
1259 // the following loop.
1260 std::vector<BoundingBox<dim>> local_boxes;
1261 for (const auto &cell : tria.active_cell_iterators())
1262 if (cell->is_locally_owned())
1263 local_boxes.emplace_back(mapping.get_bounding_box(cell));
1264
1265 // Create r-tree of bounding boxes
1266 const auto local_tree = pack_rtree(local_boxes);
1267
1268 // Compress r-tree to a minimal set of bounding boxes
1269 std::vector<std::vector<BoundingBox<dim>>> global_bboxes(1);
1270 global_bboxes[0] = extract_rtree_level(local_tree, 0);
1271
1272 const GridTools::Cache<dim, dim> cache(tria, mapping);
1273
1274 // Iterate over all sides of the non-matching interface.
1275 for (const auto &[nm_face, marked_vertices] :
1276 non_matching_faces_marked_vertices)
1277 {
1278 // 1) compute cell face pairs
1279 std::vector<
1280 std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
1281 cell_face_pairs;
1282
1283 std::vector<unsigned int> indices;
1284
1285 for (unsigned int face = face_range.first; face < face_range.second;
1286 ++face)
1287 {
1288 if (matrix_free.get_boundary_id(face) == nm_face)
1289 {
1290 for (unsigned int v = 0;
1291 v < matrix_free.n_active_entries_per_face_batch(face);
1292 ++v)
1293 {
1294 const auto &[c, f] = matrix_free.get_face_iterator(face, v);
1295
1296 cell_face_pairs.emplace_back(std::make_pair(c, f));
1297 indices.push_back(face * n_lanes + v);
1298 }
1299 }
1300 }
1301
1302 // 2) Create RPE.
1303 // In the Nitsche-type case we do not collect points for the setup
1304 // of RemotePointEvaluation. Instead we compute intersections between
1305 // the faces and set up RemotePointEvaluation with the computed
1306 // intersections.
1307
1308 // Build intersection requests. Intersection requests
1309 // correspond to vertices at faces.
1310 std::vector<std::vector<Point<dim>>> intersection_requests;
1311 for (const auto &[cell, f] : cell_face_pairs)
1312 {
1313 std::vector<Point<dim>> vertices(cell->face(f)->n_vertices());
1314 std::copy_n(mapping.get_vertices(cell, f).begin(),
1315 cell->face(f)->n_vertices(),
1316 vertices.begin());
1317 intersection_requests.emplace_back(vertices);
1318 }
1319
1320 // Compute intersection data with user specified function for marked
1321 // vertices.
1322 auto intersection_data =
1324 dim - 1>(cache,
1325 intersection_requests,
1326 global_bboxes,
1327 marked_vertices(),
1328 tolerance);
1329
1330 // Convert to RPE.
1331 std::vector<Quadrature<dim>> mapped_quadratures_recv_comp;
1332
1333 auto rpe =
1334 std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>();
1335 rpe->reinit(
1336 intersection_data
1337 .template convert_to_distributed_compute_point_locations_internal<
1338 dim>(n_q_pnts_1D, tria, mapping, &mapped_quadratures_recv_comp),
1339 tria,
1340 mapping);
1341
1342 // 3) Fill global quadrature vector.
1343 for (unsigned int i = 0; i < intersection_requests.size(); ++i)
1344 {
1345 const auto idx = indices[i];
1346
1347 // We do not use a structural binding here, since with
1348 // C++17 caputuring structural bindings in lambdas leads
1349 // to an ill formed program.
1350 const auto &cell = std::get<0>(cell_face_pairs[i]);
1351 const auto &f = std::get<1>(cell_face_pairs[i]);
1352
1353 std::vector<Point<dim - 1>> q_points;
1354 std::vector<double> weights;
1355 for (unsigned int ptr = intersection_data.recv_ptrs[i];
1356 ptr < intersection_data.recv_ptrs[i + 1];
1357 ++ptr)
1358 {
1359 const auto &quad = mapped_quadratures_recv_comp[ptr];
1360
1361 const auto &ps = quad.get_points();
1362 std::transform(
1363 ps.begin(),
1364 ps.end(),
1365 std::back_inserter(q_points),
1366 [&](const Point<dim> &p) {
1367 return mapping.project_real_point_to_unit_point_on_face(
1368 cell, f, p);
1369 });
1370
1371 const auto &ws = quad.get_weights();
1372 weights.insert(weights.end(), ws.begin(), ws.end());
1373 }
1374 Quadrature<dim - 1> quad(q_points, weights);
1375
1376 Assert(global_quadrature_vector[idx].size() == 0,
1377 ExcMessage("Quadrature for given face already provided."));
1378
1379 global_quadrature_vector[idx] = quad;
1380 }
1381
1382 // Add communication object.
1384 co.indices = indices;
1385 co.rpe = rpe;
1386 comm_objects.push_back(co);
1387 }
1388
1389 // Reinit the communicator with the communication objects.
1390 FERemoteEvaluationCommunicator<dim> remote_communicator;
1391
1392 remote_communicator.reinit_faces(
1393 comm_objects,
1394 std::make_pair(0, global_quadrature_vector.size()),
1395 global_quadrature_vector);
1396
1397 if (nm_mapping_info != nullptr)
1398 {
1399 std::vector<
1400 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
1401 vector_face_accessors;
1402 vector_face_accessors.reserve((matrix_free.n_inner_face_batches() +
1403 matrix_free.n_boundary_face_batches()) *
1404 n_lanes);
1405
1406 // fill container for inner face batches
1407 unsigned int face_batch = 0;
1408 for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch)
1409 {
1410 for (unsigned int v = 0; v < n_lanes; ++v)
1411 {
1412 if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
1413 vector_face_accessors.push_back(
1414 matrix_free.get_face_iterator(face_batch, v));
1415 else
1416 vector_face_accessors.push_back(
1417 matrix_free.get_face_iterator(face_batch, 0));
1418 }
1419 }
1420 // and boundary face batches
1421 for (; face_batch < (matrix_free.n_inner_face_batches() +
1422 matrix_free.n_boundary_face_batches());
1423 ++face_batch)
1424 {
1425 for (unsigned int v = 0; v < n_lanes; ++v)
1426 {
1427 if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
1428 vector_face_accessors.push_back(
1429 matrix_free.get_face_iterator(face_batch, v));
1430 else
1431 vector_face_accessors.push_back(
1432 matrix_free.get_face_iterator(face_batch, 0));
1433 }
1434 }
1435
1436 nm_mapping_info->reinit_faces(vector_face_accessors,
1437 global_quadrature_vector);
1438 }
1439
1440 return remote_communicator;
1441 }
1442} // namespace Utilities
1443
1444
1445
1446template <int dim, int n_components, typename value_type>
1447template <typename MeshType>
1450 const MeshType &mesh,
1451 const unsigned int first_selected_component,
1452 const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags)
1453 : comm(&comm)
1454 , first_selected_component(first_selected_component)
1455 , evaluation_flags(evaluation_flags)
1456{
1457 set_mesh(mesh);
1458}
1459
1460template <int dim, int n_components, typename value_type>
1461template <typename VectorType>
1462void
1464 const VectorType &src,
1466{
1467 if (tria)
1468 {
1469 Assert(n_components == 1, ExcNotImplemented());
1470 comm->template update_ghost_values<n_components>(this->data,
1471 *tria,
1472 src,
1473 flags,
1474 first_selected_component,
1475 evaluation_flags);
1476 }
1477 else if (dof_handler)
1478 {
1479 comm->template update_ghost_values<n_components>(this->data,
1480 *dof_handler,
1481 src,
1482 flags,
1483 first_selected_component,
1484 evaluation_flags);
1485 }
1486 else
1488}
1489
1490template <int dim, int n_components, typename value_type>
1493{
1495 comm->get_view());
1496 return data_accessor;
1497}
1498
1499template <int dim, int n_components, typename value_type>
1500void
1502 const Triangulation<dim> &tria)
1503{
1504 this->tria = &tria;
1505}
1506
1507template <int dim, int n_components, typename value_type>
1508void
1510 const DoFHandler<dim> &dof_handler)
1511{
1512 this->dof_handler = &dof_handler;
1513}
1514
1516
1517#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)
void reinit_faces(const std::vector< FERemoteCommunicationObjectEntityBatches< dim > > &comm_objects, const std::pair< unsigned int, unsigned int > &face_batch_range, const std::vector< Quadrature< dim > > &quadrature_vector)
const internal::PrecomputedEvaluationDataView & get_view() const
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
Triangulation< dim, spacedim > & get_triangulation()
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
IteratorRange< active_cell_iterator > active_cell_iterators() const
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