deal.II version GIT relicensing-3696-g88d7b26363 2025-07-08 15:10: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
mapping_info.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2025 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_non_matching_mapping_info_h
17#define dealii_non_matching_mapping_info_h
18
19
20#include <deal.II/base/config.h>
21
26
27#include <deal.II/fe/fe_dgq.h>
29#include <deal.II/fe/mapping.h>
33
35
36#include <memory>
37
38
40
41namespace NonMatching
42{
43 namespace internal
44 {
45 template <int dim, int spacedim = dim>
47 {
50 spacedim>;
51
52 public:
53 static UpdateFlags
55 const ObserverPointer<const Mapping<dim, spacedim>> &mapping,
56 const UpdateFlags &update_flags)
57 {
58 return mapping->requires_update_flags(update_flags);
59 }
60
61 static void
63 const ObserverPointer<const Mapping<dim, spacedim>> &mapping,
64 const UpdateFlags &update_flags_mapping,
66 CellSimilarity::Similarity &cell_similarity,
67 const Quadrature<dim> &quadrature,
68 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
69 &internal_mapping_data,
70 MappingData &mapping_data)
71 {
72 mapping_data.initialize(quadrature.size(), update_flags_mapping);
73 internal_mapping_data->reinit(update_flags_mapping, quadrature);
74
75 cell_similarity = mapping->fill_fe_values(cell,
76 cell_similarity,
77 quadrature,
78 *internal_mapping_data,
79 mapping_data);
80 }
81
82
83
84 static void
86 const ObserverPointer<const Mapping<dim, spacedim>> &mapping,
87 const UpdateFlags &update_flags_mapping,
89 const ImmersedSurfaceQuadrature<dim> &quadrature,
90 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
91 &internal_mapping_data,
92 MappingData &mapping_data)
93 {
94 mapping_data.initialize(quadrature.size(), update_flags_mapping);
95
96 internal_mapping_data->reinit(update_flags_mapping, quadrature);
97
98 mapping->fill_fe_immersed_surface_values(cell,
99 quadrature,
100 *internal_mapping_data,
101 mapping_data);
102 }
103
104
105
106 static void
108 const ObserverPointer<const Mapping<dim, spacedim>> &mapping,
109 const UpdateFlags &update_flags_mapping,
111 const unsigned int face_no,
112 const Quadrature<dim - 1> &quadrature,
113 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
114 &internal_mapping_data,
115 MappingData &mapping_data)
116 {
117 mapping_data.initialize(quadrature.size(), update_flags_mapping);
118
119 // reuse internal_mapping_data for MappingQ to avoid memory allocations
120 if (const MappingQ<dim, spacedim> *mapping_q =
121 dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
122 {
123 auto &data =
124 dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
125 *internal_mapping_data);
126 data.initialize_face(update_flags_mapping,
128 cell->reference_cell(),
129 quadrature,
130 face_no,
131 cell->combined_face_orientation(face_no)),
132 quadrature.size());
133
134 mapping_q->fill_mapping_data_for_face_quadrature(
135 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
136 }
137 else
138 {
139 auto internal_mapping_data =
140 mapping->get_face_data(update_flags_mapping,
141 hp::QCollection<dim - 1>(quadrature));
142
143 mapping->fill_fe_face_values(cell,
144 face_no,
145 hp::QCollection<dim - 1>(quadrature),
146 *internal_mapping_data,
147 mapping_data);
148 }
149 }
150 };
151
152 template <int dim, int spacedim = dim>
155 const double diameter,
157 &inverse_jacobians)
158 {
159 const auto jac_0 = inverse_jacobians[0];
160 const double zero_tolerance_double =
161 1e4 / diameter * std::numeric_limits<double>::epsilon() * 1024.;
162 bool jacobian_constant = true;
163 for (unsigned int q = 1; q < inverse_jacobians.size(); ++q)
164 {
165 const DerivativeForm<1, spacedim, dim> &jac = inverse_jacobians[q];
166 for (unsigned int d = 0; d < dim; ++d)
167 for (unsigned int e = 0; e < spacedim; ++e)
168 if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
169 jacobian_constant = false;
170 if (!jacobian_constant)
171 break;
172 }
173
174 // check whether the Jacobian is diagonal to machine
175 // accuracy
176 bool cell_cartesian = jacobian_constant && dim == spacedim;
177 for (unsigned int d = 0; d < dim; ++d)
178 for (unsigned int e = 0; e < dim; ++e)
179 if (d != e)
180 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
181 {
182 cell_cartesian = false;
183 break;
184 }
185
186 // return cell type
187 if (cell_cartesian)
188 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
189 else if (jacobian_constant)
190 return ::internal::MatrixFreeFunctions::GeometryType::affine;
191 else
192 return ::internal::MatrixFreeFunctions::GeometryType::general;
193 }
194 } // namespace internal
195
218 template <int dim, int spacedim = dim, typename Number = double>
220 {
221 public:
225 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
226 Number>::vectorized_value_type;
227
233 {
238 const bool store_cells = false)
241 {}
242
249
258 };
259
274 const AdditionalData additional_data = AdditionalData());
275
279 MappingInfo(const MappingInfo &) = delete;
280
285 operator=(const MappingInfo &) = delete;
286
291 void
293 const std::vector<Point<dim>> &unit_points);
294
299 void
301 const ArrayView<const Point<dim>> &unit_points);
302
308 void
311
322 template <typename ContainerType>
323 void
325 const ContainerType &cell_iterator_range,
326 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
327 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
328
335 template <typename ContainerType>
336 void
338 const ContainerType &cell_iterator_range,
339 const std::vector<Quadrature<dim>> &quadrature_vector,
340 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
341
346 template <typename ContainerType>
347 void
349 const ContainerType &cell_iterator_range,
350 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
351 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
352
357 template <typename ContainerType>
358 void
360 const ContainerType &cell_iterator_range,
361 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
362 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
363
368 template <typename CellIteratorType>
369 void
370 reinit_faces(const std::vector<std::pair<CellIteratorType, unsigned int>>
371 &face_iterator_range_interior,
372 const std::vector<Quadrature<dim - 1>> &quadrature_vector);
373
378 bool
380
384 unsigned int
385 get_face_number(const unsigned int offset, const bool is_interior) const;
386
392 get_unit_point(const unsigned int offset) const;
393
398 const Point<dim - 1, VectorizedArrayType> *
399 get_unit_point_faces(const unsigned int offset) const;
400
406 get_jacobian(const unsigned int offset,
407 const bool is_interior = true) const;
408
414 get_inverse_jacobian(const unsigned int offset,
415 const bool is_interior = true) const;
416
422 get_normal_vector(const unsigned int offset) const;
423
428 const Number *
429 get_JxW(const unsigned int offset) const;
430
436 get_real_point(const unsigned int offset) const;
437
442 get_mapping() const;
443
449
455
459 template <bool is_face>
460 unsigned int
462 const unsigned int face_number) const;
463
467 unsigned int
468 compute_unit_point_index_offset(const unsigned int geometry_index) const;
469
473 unsigned int
474 compute_data_index_offset(const unsigned int geometry_index) const;
475
479 unsigned int
481 const unsigned int geometry_index) const;
482
486 unsigned int
487 get_n_q_points_unvectorized(const unsigned int geometry_index) const;
488
493 get_cell_type(const unsigned int geometry_index) const;
494
501 get_cell_iterator(const unsigned int cell_index) const;
502
506 std::size_t
508
509 private:
512 spacedim>;
513
517 template <typename NumberType>
518 unsigned int
520
524 void
526
530 void
531 resize_unit_points(const unsigned int n_unit_point_batches);
532
536 void
537 resize_unit_points_faces(const unsigned int n_unit_point_batches);
538
542 void
543 resize_data_fields(const unsigned int n_data_point_batches,
544 const bool is_face_centric = false);
545
549 void
550 store_unit_points(const unsigned int unit_points_index_offset,
551 const unsigned int n_q_points,
552 const unsigned int n_q_points_unvectorized,
553 const std::vector<Point<dim>> &points);
554
558 void
559 store_unit_points_faces(const unsigned int unit_points_index_offset,
560 const unsigned int n_q_points,
561 const unsigned int n_q_points_unvectorized,
562 const std::vector<Point<dim - 1>> &points);
563
567 void
568 store_mapping_data(const unsigned int unit_points_index_offset,
569 const unsigned int n_q_points,
570 const unsigned int n_q_points_unvectorized,
572 const std::vector<double> &weights,
573 const unsigned int compressed_unit_point_index_offset,
574 const bool affine_cell,
575 const bool is_interior = true);
576
580 unsigned int
581 compute_compressed_cell_index(const unsigned int cell_index) const;
582
586 template <typename ContainerType, typename QuadratureType>
587 void
589 const ContainerType &cell_iterator_range,
590 const std::vector<QuadratureType> &quadrature_vector,
591 const unsigned int n_unfiltered_cells,
592 const std::function<
593 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
594 const QuadratureType &quadrature,
595 MappingData &mapping_data)> &compute_mapping_data);
596
600 enum class State
601 {
602 invalid,
607 };
608
614
621
628
632 std::vector<unsigned int> unit_points_index;
633
637 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
639
644
651
656
661
666
671
678 std::vector<::internal::MatrixFreeFunctions::GeometryType> cell_type;
679
683 std::vector<unsigned int> data_index_offsets;
684
688 std::vector<unsigned int> compressed_data_index_offsets;
689
697
704
711 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
713
721 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
723
730
734 std::vector<unsigned int> n_q_points_unvectorized;
735
740 std::vector<unsigned int> cell_index_offset;
741
746 std::vector<unsigned int> cell_index_to_compressed_cell_index;
747
752
759
764 std::vector<std::pair<int, int>> cell_level_and_indices;
765
770 std::vector<std::pair<unsigned char, unsigned char>> face_number;
771 };
772
773 template <int dim, int spacedim, typename Number>
774 inline unsigned int
776 const unsigned int offset,
777 const bool is_interior) const
778 {
779 const auto &face_pair = face_number[offset];
780 return is_interior ? face_pair.first : face_pair.second;
781 }
782
783 // ----------------------- template functions ----------------------
784
785
786 template <int dim, int spacedim, typename Number>
788 const Mapping<dim, spacedim> &mapping,
789 const UpdateFlags update_flags,
790 const AdditionalData additional_data)
791 : mapping(&mapping)
792 , update_flags(update_flags)
793 , update_flags_mapping(update_default)
794 , additional_data(additional_data)
795 {
796 // translate update flags
806
807 // always save quadrature points for now
809
812 this->mapping, update_flags_mapping);
813
814 // construct internal_mapping_data for mappings for reuse in reinit()
815 // calls to avoid frequent memory allocations
817 }
818
819
820
821 template <int dim, int spacedim, typename Number>
822 void
824 {
825 n_q_points_unvectorized.clear();
826 unit_points_index.clear();
827 data_index_offsets.clear();
828 compressed_data_index_offsets.clear();
829 cell_type.clear();
830 }
831
832
833
834 template <int dim, int spacedim, typename Number>
835 void
838 const std::vector<Point<dim>> &unit_points_in)
839 {
840 reinit(cell, make_array_view(unit_points_in));
841 }
842
843
844
845 template <int dim, int spacedim, typename Number>
846 void
849 const ArrayView<const Point<dim>> &unit_points_in)
850 {
851 quadrature.initialize(unit_points_in);
852 reinit(cell, quadrature);
853 }
854
855
856
857 template <int dim, int spacedim, typename Number>
858 void
861 const Quadrature<dim> &quadrature)
862 {
863 n_q_points_unvectorized.resize(1);
864 n_q_points_unvectorized[0] = quadrature.size();
865
866 const unsigned int n_q_points =
867 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
868
869 const unsigned int n_q_points_data =
870 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
871
872 // resize data vectors
873 resize_unit_points(n_q_points);
874 resize_data_fields(n_q_points_data);
875
876 // store unit points
877 store_unit_points(0,
878 n_q_points,
879 n_q_points_unvectorized[0],
880 quadrature.get_points());
881
882 // compute mapping data
886 update_flags_mapping,
887 cell,
888 cell_similarity,
889 quadrature,
890 internal_mapping_data,
891 mapping_data);
892
893 // check for cartesian/affine cell
894 if (!quadrature.empty() &&
895 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
896 {
897 cell_type.push_back(
898 internal::compute_geometry_type(cell->diameter(),
899 mapping_data.inverse_jacobians));
900 }
901 else
902 cell_type.push_back(
904
905 // store mapping data
906 store_mapping_data(
907 0,
908 n_q_points_data,
909 n_q_points_unvectorized[0],
910 mapping_data,
911 quadrature.get_weights(),
912 0,
913 cell_type.back() <=
915
916 unit_points_index.push_back(0);
917 data_index_offsets.push_back(0);
918 compressed_data_index_offsets.push_back(0);
919
920 state = State::single_cell;
921 }
922
923
924
925 template <int dim, int spacedim, typename Number>
926 template <typename ContainerType>
927 void
929 const ContainerType &cell_iterator_range,
930 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
931 const unsigned int n_unfiltered_cells)
932 {
933 const unsigned int n_cells = unit_points_vector.size();
934 AssertDimension(n_cells,
935 std::distance(cell_iterator_range.begin(),
936 cell_iterator_range.end()));
937
938 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
939 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
940 quadrature_vector[cell_index] =
941 Quadrature<dim>(unit_points_vector[cell_index]);
942
943 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
944 }
945
946
947
948 template <int dim, int spacedim, typename Number>
949 template <typename ContainerType, typename QuadratureType>
950 void
952 const ContainerType &cell_iterator_range,
953 const std::vector<QuadratureType> &quadrature_vector,
954 const unsigned int n_unfiltered_cells,
955 const std::function<
956 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
957 const QuadratureType &quadrature,
958 MappingData &mapping_data)> &compute_mapping_data)
959 {
960 clear();
961
962 do_cell_index_compression =
963 n_unfiltered_cells != numbers::invalid_unsigned_int;
964
965 const unsigned int n_cells = quadrature_vector.size();
966 AssertDimension(n_cells,
967 std::distance(cell_iterator_range.begin(),
968 cell_iterator_range.end()));
969
970 n_q_points_unvectorized.reserve(n_cells);
971
972 cell_type.reserve(n_cells);
973
974 if (additional_data.store_cells)
975 cell_level_and_indices.resize(n_cells);
976
977 // fill unit points index offset vector
978 unit_points_index.reserve(n_cells + 1);
979 unit_points_index.push_back(0);
980 data_index_offsets.reserve(n_cells + 1);
981 data_index_offsets.push_back(0);
982 for (const auto &quadrature : quadrature_vector)
983 {
984 const unsigned int n_points = quadrature.size();
985 n_q_points_unvectorized.push_back(n_points);
986
987 const unsigned int n_q_points =
988 compute_n_q_points<VectorizedArrayType>(n_points);
989 unit_points_index.push_back(unit_points_index.back() + n_q_points);
990
991 const unsigned int n_q_points_data =
992 compute_n_q_points<Number>(n_points);
993 data_index_offsets.push_back(data_index_offsets.back() +
994 n_q_points_data);
995 }
996
997 const unsigned int n_unit_points = unit_points_index.back();
998 const unsigned int n_data_points = data_index_offsets.back();
999
1000 // resize data vectors
1001 resize_unit_points(n_unit_points);
1002 resize_data_fields(n_data_points);
1003
1004 if (do_cell_index_compression)
1005 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1007
1008 MappingData mapping_data_previous_cell;
1009 unsigned int size_compressed_data = 0;
1010 unsigned int cell_index = 0;
1011 for (const auto &cell : cell_iterator_range)
1012 {
1013 if (additional_data.store_cells)
1014 {
1015 this->triangulation = &cell->get_triangulation();
1016 cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
1017 }
1018
1019 const auto &quadrature = quadrature_vector[cell_index];
1020 const bool empty = quadrature.empty();
1021
1022 // store unit points
1023 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1024 n_q_points_unvectorized[cell_index]);
1025 store_unit_points(unit_points_index[cell_index],
1026 n_q_points,
1027 n_q_points_unvectorized[cell_index],
1028 quadrature.get_points());
1029
1030 // compute mapping data
1031 compute_mapping_data(cell, quadrature, mapping_data);
1032
1033 // store mapping data
1034 const unsigned int n_q_points_data =
1035 compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
1036
1037 // check for cartesian/affine cell
1038 if (!empty &&
1039 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1040 {
1041 cell_type.push_back(
1042 internal::compute_geometry_type(cell->diameter(),
1043 mapping_data.inverse_jacobians));
1044 }
1045 else
1046 cell_type.push_back(
1048
1049 if (cell_index > 0)
1050 {
1051 // check if current and previous cell are affine
1052 const bool affine_cells =
1053 cell_type[cell_index] <=
1055 cell_type[cell_index - 1] <=
1057
1058 // create a comparator to compare inverse Jacobian of current
1059 // and previous cell
1061 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1062 1024.);
1063
1064 // we can only compare if current and previous cell have at least
1065 // one quadrature point and both cells are at least affine
1066 const auto comparison_result =
1067 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1068 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1070 comparator.compare(
1071 mapping_data.inverse_jacobians[0],
1072 mapping_data_previous_cell.inverse_jacobians[0]);
1073
1074 // we can compress the Jacobians and inverse Jacobians if
1075 // inverse Jacobians are equal and cells are affine
1076 if (affine_cells &&
1077 comparison_result ==
1079 {
1080 compressed_data_index_offsets.push_back(
1081 compressed_data_index_offsets.back());
1082 }
1083 else
1084 {
1085 const unsigned int n_compressed_data_last_cell =
1086 cell_type[cell_index - 1] <=
1088 1 :
1089 compute_n_q_points<Number>(
1090 n_q_points_unvectorized[cell_index - 1]);
1091
1092 compressed_data_index_offsets.push_back(
1093 compressed_data_index_offsets.back() +
1094 n_compressed_data_last_cell);
1095 }
1096 }
1097 else
1098 compressed_data_index_offsets.push_back(0);
1099
1100 // cache mapping_data from previous cell
1101 mapping_data_previous_cell = mapping_data;
1102
1103 store_mapping_data(data_index_offsets[cell_index],
1104 n_q_points_data,
1105 n_q_points_unvectorized[cell_index],
1106 mapping_data,
1107 quadrature.get_weights(),
1108 compressed_data_index_offsets[cell_index],
1109 cell_type[cell_index] <=
1111
1112 // update size of compressed data depending on cell type and handle
1113 // empty quadratures
1114 if (cell_type[cell_index] <=
1116 size_compressed_data = compressed_data_index_offsets.back() + 1;
1117 else
1118 size_compressed_data =
1119 std::max(size_compressed_data,
1120 compressed_data_index_offsets.back() + n_q_points_data);
1121
1122 if (do_cell_index_compression)
1123 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1124 cell_index;
1125
1126 ++cell_index;
1127 }
1128
1129 if (update_flags_mapping & UpdateFlags::update_jacobians)
1130 {
1131 jacobians[0].resize(size_compressed_data);
1132 jacobians[0].shrink_to_fit();
1133 }
1134 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1135 {
1136 inverse_jacobians[0].resize(size_compressed_data);
1137 inverse_jacobians[0].shrink_to_fit();
1138 }
1139
1140 state = State::cell_vector;
1141 }
1142
1143
1144
1145 template <int dim, int spacedim, typename Number>
1146 template <typename ContainerType>
1147 void
1149 const ContainerType &cell_iterator_range,
1150 const std::vector<Quadrature<dim>> &quadrature_vector,
1151 const unsigned int n_unfiltered_cells)
1152 {
1153 auto compute_mapping_data_for_cells =
1154 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1155 const Quadrature<dim> &quadrature,
1156 MappingData &mapping_data) {
1157 CellSimilarity::Similarity cell_similarity =
1161 update_flags_mapping,
1162 cell,
1163 cell_similarity,
1164 quadrature,
1165 internal_mapping_data,
1166 mapping_data);
1167 };
1168
1169 do_reinit_cells<ContainerType, Quadrature<dim>>(
1170 cell_iterator_range,
1171 quadrature_vector,
1172 n_unfiltered_cells,
1173 compute_mapping_data_for_cells);
1174 }
1175
1176
1177
1178 template <int dim, int spacedim, typename Number>
1179 template <typename ContainerType>
1180 void
1182 const ContainerType &cell_iterator_range,
1183 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
1184 const unsigned int n_unfiltered_cells)
1185 {
1186 Assert(
1187 additional_data.use_global_weights == false,
1188 ExcMessage(
1189 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1190
1191 Assert(additional_data.store_cells == false, ExcNotImplemented());
1192
1193 if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
1194 update_flags_mapping |= update_covariant_transformation;
1195
1196 auto compute_mapping_data_for_surface =
1197 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1198 const ImmersedSurfaceQuadrature<dim> &quadrature,
1199 MappingData &mapping_data) {
1202 mapping,
1203 update_flags_mapping,
1204 cell,
1205 quadrature,
1206 internal_mapping_data,
1207 mapping_data);
1208 };
1209
1210 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1211 cell_iterator_range,
1212 quadrature_vector,
1213 n_unfiltered_cells,
1214 compute_mapping_data_for_surface);
1215 }
1216
1217
1218
1219 template <int dim, int spacedim, typename Number>
1220 template <typename ContainerType>
1221 void
1223 const ContainerType &cell_iterator_range,
1224 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
1225 const unsigned int n_unfiltered_cells)
1226 {
1227 clear();
1228
1229 Assert(additional_data.store_cells == false, ExcNotImplemented());
1230
1231 do_cell_index_compression =
1232 n_unfiltered_cells != numbers::invalid_unsigned_int;
1233
1234 const unsigned int n_cells = quadrature_vector.size();
1235 AssertDimension(n_cells,
1236 std::distance(cell_iterator_range.begin(),
1237 cell_iterator_range.end()));
1238
1239 // fill cell index offset vector
1240 cell_index_offset.resize(n_cells);
1241 unsigned int n_faces = 0;
1242 unsigned int cell_index = 0;
1243 for (const auto &cell : cell_iterator_range)
1244 {
1245 cell_index_offset[cell_index] = n_faces;
1246 n_faces += cell->n_faces();
1247 ++cell_index;
1248 }
1249
1250 n_q_points_unvectorized.reserve(n_faces);
1251
1252 cell_type.reserve(n_faces);
1253
1254 // fill unit points index offset vector
1255 unit_points_index.resize(n_faces + 1);
1256 data_index_offsets.resize(n_faces + 1);
1257 cell_index = 0;
1258 unsigned int n_unit_points = 0;
1259 unsigned int n_data_points = 0;
1260 for (const auto &cell : cell_iterator_range)
1261 {
1262 for (const auto &f : cell->face_indices())
1263 {
1264 const unsigned int current_face_index =
1265 cell_index_offset[cell_index] + f;
1266
1267 unit_points_index[current_face_index] = n_unit_points;
1268 data_index_offsets[current_face_index] = n_data_points;
1269
1270 const unsigned int n_points =
1271 quadrature_vector[cell_index][f].size();
1272 n_q_points_unvectorized.push_back(n_points);
1273
1274 const unsigned int n_q_points =
1275 compute_n_q_points<VectorizedArrayType>(n_points);
1276 n_unit_points += n_q_points;
1277
1278 const unsigned int n_q_points_data =
1279 compute_n_q_points<Number>(n_points);
1280 n_data_points += n_q_points_data;
1281 }
1282
1283 ++cell_index;
1284 }
1285 unit_points_index[n_faces] = n_unit_points;
1286 data_index_offsets[n_faces] = n_data_points;
1287
1288 // compress indices
1289 if (do_cell_index_compression)
1290 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1292
1293 // fill unit points and mapping data for every face of all cells
1294 // resize data vectors
1295 resize_unit_points_faces(n_unit_points);
1296 resize_data_fields(n_data_points);
1297
1298 MappingData mapping_data_previous_cell;
1299 MappingData mapping_data_first;
1300 bool first_set = false;
1301 unsigned int size_compressed_data = 0;
1302 cell_index = 0;
1303 for (const auto &cell : cell_iterator_range)
1304 {
1305 const auto &quadratures_on_faces = quadrature_vector[cell_index];
1306
1307 Assert(quadratures_on_faces.size() == cell->n_faces(),
1308 ExcDimensionMismatch(quadratures_on_faces.size(),
1309 cell->n_faces()));
1310
1311 for (const auto &f : cell->face_indices())
1312 {
1313 const auto &quadrature_on_face = quadratures_on_faces[f];
1314 const bool empty = quadrature_on_face.empty();
1315
1316 const unsigned int current_face_index =
1317 cell_index_offset[cell_index] + f;
1318
1319 // store unit points
1320 const unsigned int n_q_points =
1321 compute_n_q_points<VectorizedArrayType>(
1322 n_q_points_unvectorized[current_face_index]);
1323 store_unit_points_faces(unit_points_index[current_face_index],
1324 n_q_points,
1325 n_q_points_unvectorized[current_face_index],
1326 quadrature_on_face.get_points());
1327
1330 update_flags_mapping,
1331 cell,
1332 f,
1333 quadrature_on_face,
1334 internal_mapping_data,
1335 mapping_data);
1336
1337 // check for cartesian/affine cell
1338 if (!empty &&
1339 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1340 {
1341 cell_type.push_back(internal::compute_geometry_type(
1342 cell->diameter(), mapping_data.inverse_jacobians));
1343
1344 if (!first_set)
1345 {
1346 mapping_data_first = mapping_data;
1347 first_set = true;
1348 }
1349 }
1350 else
1351 cell_type.push_back(
1353
1354 if (current_face_index > 0)
1355 {
1356 // check if current and previous cell are affine
1357 const bool affine_cells =
1358 cell_type[current_face_index] <=
1360 cell_type[current_face_index - 1] <=
1362
1363 // create a comparator to compare inverse Jacobian of current
1364 // and previous cell
1366 1e4 / cell->diameter() *
1367 std::numeric_limits<double>::epsilon() * 1024.);
1368
1369 // we can only compare if current and previous cell have at
1370 // least one quadrature point and both cells are at least affine
1371 const auto comparison_result =
1372 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1373 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1375 comparator.compare(
1376 mapping_data.inverse_jacobians[0],
1377 mapping_data_previous_cell.inverse_jacobians[0]);
1378
1379 // we can compress the Jacobians and inverse Jacobians if
1380 // inverse Jacobians are equal and cells are affine
1381 if (affine_cells &&
1382 comparison_result ==
1384 {
1385 compressed_data_index_offsets.push_back(
1386 compressed_data_index_offsets.back());
1387 }
1388 else if (first_set &&
1389 (cell_type[current_face_index] <=
1391 (comparator.compare(
1392 mapping_data.inverse_jacobians[0],
1393 mapping_data_first.inverse_jacobians[0]) ==
1395 double>::ComparisonResult::equal))
1396 {
1397 compressed_data_index_offsets.push_back(0);
1398 }
1399 else
1400 {
1401 const unsigned int n_compressed_data_last_cell =
1402 cell_type[current_face_index - 1] <=
1404 1 :
1405 compute_n_q_points<Number>(
1406 n_q_points_unvectorized[current_face_index - 1]);
1407
1408 compressed_data_index_offsets.push_back(
1409 compressed_data_index_offsets.back() +
1410 n_compressed_data_last_cell);
1411 }
1412 }
1413 else
1414 compressed_data_index_offsets.push_back(0);
1415
1416 // cache mapping_data from previous cell
1417 mapping_data_previous_cell = mapping_data;
1418
1419 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1420 n_q_points_unvectorized[current_face_index]);
1421 store_mapping_data(data_index_offsets[current_face_index],
1422 n_q_points_data,
1423 n_q_points_unvectorized[current_face_index],
1424 mapping_data,
1425 quadrature_on_face.get_weights(),
1426 data_index_offsets[current_face_index],
1427 cell_type[current_face_index] <=
1429
1430 // update size of compressed data depending on cell type and handle
1431 // empty quadratures
1432 if (cell_type[current_face_index] <=
1434 size_compressed_data = compressed_data_index_offsets.back() + 1;
1435 else
1436 size_compressed_data =
1437 std::max(size_compressed_data,
1438 compressed_data_index_offsets.back() +
1439 n_q_points_data);
1440 }
1441 if (do_cell_index_compression)
1442 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1443 cell_index;
1444
1445 ++cell_index;
1446 }
1447
1448 if (update_flags_mapping & UpdateFlags::update_jacobians)
1449 {
1450 jacobians[0].resize(size_compressed_data);
1451 jacobians[0].shrink_to_fit();
1452 }
1453 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1454 {
1455 inverse_jacobians[0].resize(size_compressed_data);
1456 inverse_jacobians[0].shrink_to_fit();
1457 }
1458
1459 state = State::faces_on_cells_in_vector;
1460 }
1461
1462
1463
1464 template <int dim, int spacedim, typename Number>
1465 template <typename CellIteratorType>
1466 void
1468 const std::vector<std::pair<CellIteratorType, unsigned int>>
1469 &face_iterator_range_interior,
1470 const std::vector<Quadrature<dim - 1>> &quadrature_vector)
1471 {
1472 clear();
1473
1474 do_cell_index_compression = false;
1475
1476 Assert(additional_data.store_cells == false, ExcNotImplemented());
1477
1478
1479 const unsigned int n_faces = quadrature_vector.size();
1480 AssertDimension(n_faces,
1481 std::distance(face_iterator_range_interior.begin(),
1482 face_iterator_range_interior.end()));
1483
1484 n_q_points_unvectorized.reserve(n_faces);
1485
1486 cell_type.reserve(n_faces);
1487 face_number.reserve(n_faces);
1488
1489 // fill unit points index offset vector
1490 unit_points_index.reserve(n_faces + 1);
1491 unit_points_index.push_back(0);
1492 data_index_offsets.reserve(n_faces + 1);
1493 data_index_offsets.push_back(0);
1494 for (const auto &quadrature : quadrature_vector)
1495 {
1496 const unsigned int n_points = quadrature.size();
1497 n_q_points_unvectorized.push_back(n_points);
1498
1499 const unsigned int n_q_points =
1500 compute_n_q_points<VectorizedArrayType>(n_points);
1501 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1502
1503 const unsigned int n_q_points_data =
1504 compute_n_q_points<Number>(n_points);
1505 data_index_offsets.push_back(data_index_offsets.back() +
1506 n_q_points_data);
1507 }
1508
1509 const unsigned int n_unit_points = unit_points_index.back();
1510 const unsigned int n_data_points = data_index_offsets.back();
1511
1512 // resize data vectors
1513 resize_unit_points_faces(n_unit_points);
1514 resize_data_fields(n_data_points, true);
1515
1516 std::array<MappingData, 2> mapping_data;
1517 std::array<MappingData, 2> mapping_data_previous_cell;
1518 std::array<MappingData, 2> mapping_data_first;
1519 bool first_set = false;
1520 unsigned int size_compressed_data = 0;
1521 unsigned int face_index = 0;
1522 for (const auto &cell_and_f : face_iterator_range_interior)
1523 {
1524 const auto &quadrature_on_face = quadrature_vector[face_index];
1525 const bool empty = quadrature_on_face.empty();
1526
1527 // get interior cell and face number
1528 const auto &cell_m = cell_and_f.first;
1529 const auto f_m = cell_and_f.second;
1530
1531 // get exterior cell and face number
1532 const auto &cell_p =
1533 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1534 const auto f_p =
1535 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_face_no(f_m);
1536
1537 Assert(
1538 empty || (cell_m->level() == cell_p->level()),
1539 ExcMessage(
1540 "Intersected faces with quadrature points need to have the same "
1541 "refinement level!"));
1542
1543 face_number.emplace_back(f_m, f_p);
1544
1545 Assert(
1546 cell_m->combined_face_orientation(f_m) ==
1548 cell_p->combined_face_orientation(f_p) ==
1550 ExcMessage(
1551 "Non standard face orientation is currently not implemented."));
1552
1553 // store unit points
1554 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1555 n_q_points_unvectorized[face_index]);
1556 store_unit_points_faces(unit_points_index[face_index],
1557 n_q_points,
1558 n_q_points_unvectorized[face_index],
1559 quadrature_on_face.get_points());
1560
1561 // compute mapping for interior face
1564 update_flags_mapping,
1565 cell_m,
1566 f_m,
1567 quadrature_on_face,
1568 internal_mapping_data,
1569 mapping_data[0]);
1570
1571 // compute mapping for exterior face
1574 update_flags_mapping,
1575 cell_p,
1576 f_p,
1577 quadrature_on_face,
1578 internal_mapping_data,
1579 mapping_data[1]);
1580
1581 // check for cartesian/affine cell
1582 if (!empty &&
1583 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1584 {
1585 // select more general type of interior and exterior cell
1586 cell_type.push_back(std::max(
1588 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1590 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1591
1592 // cache mapping data of first cell pair with non-empty quadrature
1593 // on the face
1594 if (!first_set)
1595 {
1596 mapping_data_first = mapping_data;
1597 first_set = true;
1598 }
1599 }
1600 else
1601 cell_type.push_back(
1603
1604 if (face_index > 0)
1605 {
1606 // check if current and previous cell pairs are affine
1607 const bool affine_cells =
1608 cell_type[face_index] <=
1610 cell_type[face_index - 1] <=
1612
1613 // create a comparator to compare inverse Jacobian of current
1614 // and previous cell pair
1616 1e4 / cell_m->diameter() *
1617 std::numeric_limits<double>::epsilon() * 1024.);
1618
1619 // we can only compare if current and previous cell have at
1620 // least one quadrature point and both cells are at least affine
1621 const auto comparison_result_m =
1622 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1623 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1625 comparator.compare(
1626 mapping_data[0].inverse_jacobians[0],
1627 mapping_data_previous_cell[0].inverse_jacobians[0]);
1628
1629 const auto comparison_result_p =
1630 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1631 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1633 comparator.compare(
1634 mapping_data[1].inverse_jacobians[0],
1635 mapping_data_previous_cell[1].inverse_jacobians[0]);
1636
1637 // we can compress the Jacobians and inverse Jacobians if
1638 // inverse Jacobians are equal and cells are affine
1639 if (affine_cells &&
1640 comparison_result_m ==
1642 comparison_result_p ==
1644 {
1645 compressed_data_index_offsets.push_back(
1646 compressed_data_index_offsets.back());
1647 }
1648 else if (first_set &&
1649 (cell_type[face_index] <=
1651 (comparator.compare(
1652 mapping_data[0].inverse_jacobians[0],
1653 mapping_data_first[0].inverse_jacobians[0]) ==
1655 double>::ComparisonResult::equal) &&
1656 (comparator.compare(
1657 mapping_data[1].inverse_jacobians[0],
1658 mapping_data_first[1].inverse_jacobians[0]) ==
1660 {
1661 compressed_data_index_offsets.push_back(0);
1662 }
1663 else
1664 {
1665 const unsigned int n_compressed_data_last_cell =
1666 cell_type[face_index - 1] <=
1668 1 :
1669 compute_n_q_points<Number>(
1670 n_q_points_unvectorized[face_index - 1]);
1671
1672 compressed_data_index_offsets.push_back(
1673 compressed_data_index_offsets.back() +
1674 n_compressed_data_last_cell);
1675 }
1676 }
1677 else
1678 compressed_data_index_offsets.push_back(0);
1679
1680 // cache mapping_data from previous cell pair
1681 mapping_data_previous_cell = mapping_data;
1682
1683 const unsigned int n_q_points_data =
1684 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1685
1686 // store mapping data of interior face
1687 store_mapping_data(data_index_offsets[face_index],
1688 n_q_points_data,
1689 n_q_points_unvectorized[face_index],
1690 mapping_data[0],
1691 quadrature_on_face.get_weights(),
1692 data_index_offsets[face_index],
1693 cell_type[face_index] <=
1695 true);
1696
1697 // store only necessary mapping data for exterior face (Jacobians and
1698 // inverse Jacobians)
1699 store_mapping_data(data_index_offsets[face_index],
1700 n_q_points_data,
1701 n_q_points_unvectorized[face_index],
1702 mapping_data[1],
1703 quadrature_on_face.get_weights(),
1704 data_index_offsets[face_index],
1705 cell_type[face_index] <=
1707 false);
1708
1709 // update size of compressed data depending on cell type and handle
1710 // empty quadratures
1711 if (cell_type[face_index] <=
1713 size_compressed_data = compressed_data_index_offsets.back() + 1;
1714 else
1715 size_compressed_data =
1716 std::max(size_compressed_data,
1717 compressed_data_index_offsets.back() + n_q_points_data);
1718
1719 ++face_index;
1720 }
1721
1722 if (update_flags_mapping & UpdateFlags::update_jacobians)
1723 {
1724 jacobians[0].resize(size_compressed_data);
1725 jacobians[0].shrink_to_fit();
1726 jacobians[1].resize(size_compressed_data);
1727 jacobians[1].shrink_to_fit();
1728 }
1729 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1730 {
1731 inverse_jacobians[0].resize(size_compressed_data);
1732 inverse_jacobians[0].shrink_to_fit();
1733 inverse_jacobians[1].resize(size_compressed_data);
1734 inverse_jacobians[1].shrink_to_fit();
1735 }
1736
1737 state = State::face_vector;
1738 }
1739
1740
1741
1742 template <int dim, int spacedim, typename Number>
1743 bool
1745 {
1746 return state == State::faces_on_cells_in_vector;
1747 }
1748
1749
1750
1751 template <int dim, int spacedim, typename Number>
1752 unsigned int
1754 const unsigned int geometry_index) const
1755 {
1756 return n_q_points_unvectorized[geometry_index];
1757 }
1758
1759
1760 template <int dim, int spacedim, typename Number>
1763 const unsigned int geometry_index) const
1764 {
1765 return cell_type[geometry_index];
1766 }
1767
1768
1769
1770 template <int dim, int spacedim, typename Number>
1773 const unsigned int cell_index) const
1774 {
1775 Assert(
1776 additional_data.store_cells,
1777 ExcMessage(
1778 "Cells have been not stored. You can enable this by Additional::store_cells."));
1779 return {triangulation.get(),
1780 cell_level_and_indices[cell_index].first,
1781 cell_level_and_indices[cell_index].second};
1782 }
1783
1784
1785
1786 template <int dim, int spacedim, typename Number>
1787 template <typename NumberType>
1788 unsigned int
1790 const unsigned int n_q_points_unvectorized)
1791 {
1792 const unsigned int n_lanes =
1794 const unsigned int n_filled_lanes_last_batch =
1795 n_q_points_unvectorized % n_lanes;
1796 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1797 if (n_filled_lanes_last_batch > 0)
1798 ++n_q_points;
1799 return n_q_points;
1800 }
1801
1802
1803
1804 template <int dim, int spacedim, typename Number>
1805 template <bool is_face>
1806 unsigned int
1808 const unsigned int cell_index,
1809 const unsigned int face_number) const
1810 {
1812 return 0;
1813
1814 const unsigned int compressed_cell_index =
1815 compute_compressed_cell_index(cell_index);
1816 if (!is_face)
1817 {
1818 Assert(state == State::cell_vector,
1819 ExcMessage(
1820 "This mapping info is not reinitialized for a cell vector!"));
1821 return compressed_cell_index;
1822 }
1823 else
1824 {
1826 ExcMessage(
1827 "cell_index has to be set if face_number is specified!"));
1828 Assert(state == State::faces_on_cells_in_vector ||
1829 state == State::face_vector,
1830 ExcMessage("This mapping info is not reinitialized for faces"
1831 " on cells in a vector!"));
1832 if (state == State::faces_on_cells_in_vector)
1833 return cell_index_offset[compressed_cell_index] + face_number;
1834 else if (state == State::face_vector)
1835 return cell_index;
1836 }
1837 }
1838
1839
1840
1841 template <int dim, int spacedim, typename Number>
1842 unsigned int
1844 const unsigned int cell_index) const
1845 {
1846 if (do_cell_index_compression)
1847 {
1848 Assert(cell_index_to_compressed_cell_index[cell_index] !=
1850 ExcMessage("Mapping info object was not initialized for this"
1851 " active cell index!"));
1852 return cell_index_to_compressed_cell_index[cell_index];
1853 }
1854 else
1855 return cell_index;
1856 }
1857
1858
1859 template <int dim, int spacedim, typename Number>
1860 void
1862 const unsigned int unit_points_index_offset,
1863 const unsigned int n_q_points,
1864 const unsigned int n_q_points_unvectorized,
1865 const std::vector<Point<dim>> &points)
1866 {
1867 const unsigned int n_lanes =
1869
1870 for (unsigned int q = 0; q < n_q_points; ++q)
1871 {
1872 const unsigned int offset = unit_points_index_offset + q;
1873 for (unsigned int v = 0;
1874 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1875 ++v)
1876 for (unsigned int d = 0; d < dim; ++d)
1878 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1879 }
1880 }
1881
1882
1883
1884 template <int dim, int spacedim, typename Number>
1885 void
1887 const unsigned int unit_points_index_offset,
1888 const unsigned int n_q_points,
1889 const unsigned int n_q_points_unvectorized,
1890 const std::vector<Point<dim - 1>> &points)
1891 {
1892 const unsigned int n_lanes =
1894
1895 for (unsigned int q = 0; q < n_q_points; ++q)
1896 {
1897 const unsigned int offset = unit_points_index_offset + q;
1898 for (unsigned int v = 0;
1899 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1900 ++v)
1901 for (unsigned int d = 0; d < dim - 1; ++d)
1903 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1904 }
1905 }
1906
1907
1908
1909 template <int dim, int spacedim, typename Number>
1910 void
1912 const unsigned int unit_points_index_offset,
1913 const unsigned int n_q_points,
1914 const unsigned int n_q_points_unvectorized,
1915 const MappingInfo::MappingData &mapping_data,
1916 const std::vector<double> &weights,
1917 const unsigned int compressed_unit_point_index_offset,
1918 const bool affine_cell,
1919 const bool is_interior)
1920 {
1921 const unsigned int n_lanes =
1923
1924 for (unsigned int q = 0; q < n_q_points; ++q)
1925 {
1926 const unsigned int offset = unit_points_index_offset + q;
1927 const unsigned int compressed_offset =
1928 compressed_unit_point_index_offset + q;
1929 for (unsigned int v = 0;
1930 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1931 ++v)
1932 {
1933 if (q == 0 || !affine_cell)
1934 {
1935 if (update_flags_mapping & UpdateFlags::update_jacobians)
1936 for (unsigned int d = 0; d < dim; ++d)
1937 for (unsigned int s = 0; s < spacedim; ++s)
1939 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1940 v) = mapping_data.jacobians[q * n_lanes + v][d][s];
1941 if (update_flags_mapping &
1943 for (unsigned int d = 0; d < dim; ++d)
1944 for (unsigned int s = 0; s < spacedim; ++s)
1946 inverse_jacobians[is_interior ? 0 : 1]
1947 [compressed_offset][d][s],
1948 v) =
1949 mapping_data.inverse_jacobians[q * n_lanes + v][d][s];
1950 }
1951
1952 if (is_interior)
1953 {
1954 if (update_flags_mapping & UpdateFlags::update_JxW_values)
1955 {
1956 if (additional_data.use_global_weights)
1957 {
1959 JxW_values[offset], v) = weights[q * n_lanes + v];
1960 }
1961 else
1962 {
1964 JxW_values[offset], v) =
1965 mapping_data.JxW_values[q * n_lanes + v];
1966 }
1967 }
1968 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
1969 for (unsigned int s = 0; s < spacedim; ++s)
1971 normal_vectors[offset][s], v) =
1972 mapping_data.normal_vectors[q * n_lanes + v][s];
1973 if (update_flags_mapping &
1975 for (unsigned int s = 0; s < spacedim; ++s)
1977 real_points[offset][s], v) =
1978 mapping_data.quadrature_points[q * n_lanes + v][s];
1979 }
1980 }
1981 }
1982 }
1983
1984
1985
1986 template <int dim, int spacedim, typename Number>
1987 void
1989 const unsigned int n_unit_point_batches)
1990 {
1991 unit_points.resize(n_unit_point_batches);
1992 }
1993
1994
1995
1996 template <int dim, int spacedim, typename Number>
1997 void
1999 const unsigned int n_unit_point_batches)
2000 {
2001 unit_points_faces.resize(n_unit_point_batches);
2002 }
2003
2004
2005
2006 template <int dim, int spacedim, typename Number>
2007 void
2009 const unsigned int n_data_point_batches,
2010 const bool is_face_centric)
2011 {
2012 if (update_flags_mapping & UpdateFlags::update_jacobians)
2013 {
2014 jacobians[0].resize(n_data_point_batches);
2015 if (is_face_centric)
2016 jacobians[1].resize(n_data_point_batches);
2017 }
2018 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2019 {
2020 inverse_jacobians[0].resize(n_data_point_batches);
2021 if (is_face_centric)
2022 inverse_jacobians[1].resize(n_data_point_batches);
2023 }
2024 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2025 JxW_values.resize(n_data_point_batches);
2026 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2027 normal_vectors.resize(n_data_point_batches);
2028 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2029 real_points.resize(n_data_point_batches);
2030 }
2031
2032
2033
2034 template <int dim, int spacedim, typename Number>
2035 inline const Point<
2036 dim,
2039 const unsigned int offset) const
2040 {
2041 return unit_points.data() + offset;
2042 }
2043
2044
2045
2046 template <int dim, int spacedim, typename Number>
2047 inline const Point<
2048 dim - 1,
2051 const unsigned int offset) const
2052 {
2053 return unit_points_faces.data() + offset;
2054 }
2055
2056
2057
2058 template <int dim, int spacedim, typename Number>
2059 inline const Point<spacedim, Number> *
2061 const unsigned int offset) const
2062 {
2063 return real_points.data() + offset;
2064 }
2065
2066
2067
2068 template <int dim, int spacedim, typename Number>
2069 unsigned int
2071 const unsigned int geometry_index) const
2072 {
2073 return unit_points_index[geometry_index];
2074 }
2075
2076
2077
2078 template <int dim, int spacedim, typename Number>
2079 unsigned int
2081 const unsigned int geometry_index) const
2082 {
2083 return data_index_offsets[geometry_index];
2084 }
2085
2086
2087 template <int dim, int spacedim, typename Number>
2088 unsigned int
2090 const unsigned int geometry_index) const
2091 {
2092 return compressed_data_index_offsets[geometry_index];
2093 }
2094
2095
2096
2097 template <int dim, int spacedim, typename Number>
2100 const bool is_interior) const
2101 {
2102 return jacobians[is_interior ? 0 : 1].data() + offset;
2103 }
2104
2105
2106
2107 template <int dim, int spacedim, typename Number>
2110 const unsigned int offset,
2111 const bool is_interior) const
2112 {
2113 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2114 }
2115
2116
2117
2118 template <int dim, int spacedim, typename Number>
2119 inline const Tensor<1, spacedim, Number> *
2121 const unsigned int offset) const
2122 {
2123 return normal_vectors.data() + offset;
2124 }
2125
2126
2127
2128 template <int dim, int spacedim, typename Number>
2129 inline const Number *
2130 MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
2131 {
2132 return JxW_values.data() + offset;
2133 }
2134
2135
2136
2137 template <int dim, int spacedim, typename Number>
2140 {
2141 return *mapping;
2142 }
2143
2144
2145
2146 template <int dim, int spacedim, typename Number>
2149 {
2150 return update_flags;
2151 }
2152
2153
2154
2155 template <int dim, int spacedim, typename Number>
2158 {
2159 return update_flags_mapping;
2160 }
2161
2162
2163
2164 template <int dim, int spacedim, typename Number>
2165 std::size_t
2167 {
2168 std::size_t memory = MemoryConsumption::memory_consumption(unit_points);
2169 memory += MemoryConsumption::memory_consumption(unit_points_faces);
2170 memory += MemoryConsumption::memory_consumption(unit_points_index);
2171 memory += cell_type.capacity() *
2173 memory += MemoryConsumption::memory_consumption(data_index_offsets);
2174 memory +=
2175 MemoryConsumption::memory_consumption(compressed_data_index_offsets);
2176 memory += MemoryConsumption::memory_consumption(JxW_values);
2177 memory += MemoryConsumption::memory_consumption(normal_vectors);
2178 memory += MemoryConsumption::memory_consumption(jacobians);
2179 memory += MemoryConsumption::memory_consumption(inverse_jacobians);
2180 memory += MemoryConsumption::memory_consumption(real_points);
2181 memory += MemoryConsumption::memory_consumption(n_q_points_unvectorized);
2182 memory += MemoryConsumption::memory_consumption(cell_index_offset);
2184 cell_index_to_compressed_cell_index);
2185 memory += MemoryConsumption::memory_consumption(cell_level_and_indices);
2186 memory += sizeof(*this);
2187 return memory;
2188 }
2189} // namespace NonMatching
2190
2192
2193#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:159
Abstract base class for mapping classes.
Definition mapping.h:320
const UpdateFlags update_flags
unsigned int compute_compressed_data_index_offset(const unsigned int geometry_index) const
const DerivativeForm< 1, spacedim, dim, Number > * get_inverse_jacobian(const unsigned int offset, const bool is_interior=true) const
AlignedVector< Point< dim, VectorizedArrayType > > unit_points
std::array< AlignedVector< DerivativeForm< 1, dim, spacedim, Number > >, 2 > jacobians
AlignedVector< Number > JxW_values
::internal::MatrixFreeFunctions::GeometryType get_cell_type(const unsigned int geometry_index) const
unsigned int get_face_number(const unsigned int offset, const bool is_interior) const
const Mapping< dim, spacedim > & get_mapping() 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)
std::array< AlignedVector< DerivativeForm< 1, spacedim, dim, Number > >, 2 > inverse_jacobians
const Point< dim, VectorizedArrayType > * get_unit_point(const unsigned int offset) const
unsigned int get_n_q_points_unvectorized(const unsigned int geometry_index) const
const ObserverPointer< const Mapping< dim, spacedim > > mapping
unsigned int compute_compressed_cell_index(const unsigned int cell_index) const
const Point< dim - 1, VectorizedArrayType > * get_unit_point_faces(const unsigned int offset) const
Quadrature< dim > quadrature
std::vector< unsigned int > cell_index_offset
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Quadrature< dim > &quadrature)
void store_mapping_data(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const MappingData &mapping_data, const std::vector< double > &weights, const unsigned int compressed_unit_point_index_offset, const bool affine_cell, const bool is_interior=true)
AlignedVector< Tensor< 1, spacedim, Number > > normal_vectors
unsigned int compute_geometry_index_offset(const unsigned int cell_index, const unsigned int face_number) const
std::vector<::internal::MatrixFreeFunctions::GeometryType > cell_type
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim > > &unit_points)
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< Quadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
MappingInfo & operator=(const MappingInfo &)=delete
ObserverPointer< const Triangulation< dim, spacedim > > triangulation
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > internal_mapping_data
void reinit_surface(const ContainerType &cell_iterator_range, const std::vector< ImmersedSurfaceQuadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_unit_points_faces(const unsigned int n_unit_point_batches)
const Tensor< 1, spacedim, Number > * get_normal_vector(const unsigned int offset) const
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const std::vector< Point< dim > > &unit_points)
const DerivativeForm< 1, dim, spacedim, Number > * get_jacobian(const unsigned int offset, const bool is_interior=true) const
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< std::vector< Point< dim > > > &unit_points_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_data_fields(const unsigned int n_data_point_batches, const bool is_face_centric=false)
void resize_unit_points(const unsigned int n_unit_point_batches)
std::size_t memory_consumption() const
UpdateFlags get_update_flags() const
void store_unit_points_faces(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const std::vector< Point< dim - 1 > > &points)
unsigned int compute_unit_point_index_offset(const unsigned int geometry_index) const
std::vector< unsigned int > cell_index_to_compressed_cell_index
AlignedVector< Point< spacedim, Number > > real_points
unsigned int compute_data_index_offset(const unsigned int geometry_index) const
std::vector< std::pair< unsigned char, unsigned char > > face_number
MappingInfo(const Mapping< dim, spacedim > &mapping, const UpdateFlags update_flags, const AdditionalData additional_data=AdditionalData())
std::vector< unsigned int > unit_points_index
std::vector< std::pair< int, int > > cell_level_and_indices
std::vector< unsigned int > compressed_data_index_offsets
std::vector< unsigned int > data_index_offsets
std::vector< unsigned int > n_q_points_unvectorized
const Number * get_JxW(const unsigned int offset) const
const Point< spacedim, Number > * get_real_point(const unsigned int offset) const
void store_unit_points(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const std::vector< Point< dim > > &points)
const AdditionalData additional_data
AlignedVector< Point< dim - 1, VectorizedArrayType > > unit_points_faces
unsigned int compute_n_q_points(const unsigned int n_q_points_unvectorized)
MappingInfo(const MappingInfo &)=delete
void reinit_faces(const std::vector< std::pair< CellIteratorType, unsigned int > > &face_iterator_range_interior, const std::vector< Quadrature< dim - 1 > > &quadrature_vector)
UpdateFlags update_flags_mapping
void do_reinit_cells(const ContainerType &cell_iterator_range, const std::vector< QuadratureType > &quadrature_vector, const unsigned int n_unfiltered_cells, const std::function< void(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const QuadratureType &quadrature, MappingData &mapping_data)> &compute_mapping_data)
UpdateFlags get_update_flags_mapping() const
Triangulation< dim, spacedim >::cell_iterator get_cell_iterator(const unsigned int cell_index) const
static void compute_mapping_data_for_immersed_surface_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ImmersedSurfaceQuadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_face_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, CellSimilarity::Similarity &cell_similarity, const Quadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static UpdateFlags required_update_flags(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags)
Definition point.h:113
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
bool empty() const
unsigned int size() const
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
UpdateFlags
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
std::vector< index_type > data
Definition mpi.cc:750
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
::internal::MatrixFreeFunctions::GeometryType compute_geometry_type(const double diameter, const std::vector< DerivativeForm< 1, spacedim, dim, double > > &inverse_jacobians)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
ComparisonResult compare(const std::vector< T > &v1, const std::vector< T > &v2) const
AdditionalData(const bool use_global_weights=false, const bool store_cells=false)
static constexpr std::size_t width()
static value_type & get(value_type &value, unsigned int c)