Reference documentation for deal.II version GIT relicensing-453-g9b4d7b6a5a 2024-04-20 16: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
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 - 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_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 SmartPointer<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 SmartPointer<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 SmartPointer<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 SmartPointer<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 ReferenceCells::get_hypercube<dim>(),
129 quadrature,
130 face_no,
131 cell->face_orientation(face_no),
132 cell->face_flip(face_no),
133 cell->face_rotation(face_no)),
134 quadrature.size());
135
136 mapping_q->fill_mapping_data_for_face_quadrature(
137 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
138 }
139 else
140 {
141 auto internal_mapping_data =
142 mapping->get_face_data(update_flags_mapping,
143 hp::QCollection<dim - 1>(quadrature));
144
145 mapping->fill_fe_face_values(cell,
146 face_no,
147 hp::QCollection<dim - 1>(quadrature),
148 *internal_mapping_data,
149 mapping_data);
150 }
151 }
152 };
153
154 template <int dim, int spacedim = dim>
157 const double diameter,
159 &inverse_jacobians)
160 {
161 const auto jac_0 = inverse_jacobians[0];
162 const double zero_tolerance_double =
163 1e4 / diameter * std::numeric_limits<double>::epsilon() * 1024.;
164 bool jacobian_constant = true;
165 for (unsigned int q = 1; q < inverse_jacobians.size(); ++q)
166 {
167 const DerivativeForm<1, spacedim, dim> &jac = inverse_jacobians[q];
168 for (unsigned int d = 0; d < dim; ++d)
169 for (unsigned int e = 0; e < spacedim; ++e)
170 if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
171 jacobian_constant = false;
172 if (!jacobian_constant)
173 break;
174 }
175
176 // check whether the Jacobian is diagonal to machine
177 // accuracy
178 bool cell_cartesian = jacobian_constant && dim == spacedim;
179 for (unsigned int d = 0; d < dim; ++d)
180 for (unsigned int e = 0; e < dim; ++e)
181 if (d != e)
182 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
183 {
184 cell_cartesian = false;
185 break;
186 }
187
188 // return cell type
189 if (cell_cartesian)
190 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
191 else if (jacobian_constant)
192 return ::internal::MatrixFreeFunctions::GeometryType::affine;
193 else
194 return ::internal::MatrixFreeFunctions::GeometryType::general;
195 }
196 } // namespace internal
197
220 template <int dim, int spacedim = dim, typename Number = double>
222 {
223 public:
227 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
228 Number>::vectorized_value_type;
229
235 {
240 const bool store_cells = false)
243 {}
244
251
260 };
261
276 const AdditionalData additional_data = AdditionalData());
277
281 MappingInfo(const MappingInfo &) = delete;
282
287 operator=(const MappingInfo &) = delete;
288
293 void
295 const std::vector<Point<dim>> &unit_points);
296
301 void
303 const ArrayView<const Point<dim>> &unit_points);
304
310 void
313
324 template <typename ContainerType>
325 void
327 const ContainerType &cell_iterator_range,
328 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
329 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
330
337 template <typename ContainerType>
338 void
340 const ContainerType &cell_iterator_range,
341 const std::vector<Quadrature<dim>> &quadrature_vector,
342 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
343
348 template <typename ContainerType>
349 void
351 const ContainerType &cell_iterator_range,
352 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
353 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
354
359 template <typename ContainerType>
360 void
362 const ContainerType &cell_iterator_range,
363 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
364 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
365
370 template <typename CellIteratorType>
371 void
372 reinit_faces(const std::vector<std::pair<CellIteratorType, unsigned int>>
373 &face_iterator_range_interior,
374 const std::vector<Quadrature<dim - 1>> &quadrature_vector);
375
380 bool
382
386 unsigned int
387 get_face_number(const unsigned int offset, const bool is_interior) const;
388
394 get_unit_point(const unsigned int offset) const;
395
400 const Point<dim - 1, VectorizedArrayType> *
401 get_unit_point_faces(const unsigned int offset) const;
402
408 get_jacobian(const unsigned int offset,
409 const bool is_interior = true) const;
410
416 get_inverse_jacobian(const unsigned int offset,
417 const bool is_interior = true) const;
418
424 get_normal_vector(const unsigned int offset) const;
425
430 const Number *
431 get_JxW(const unsigned int offset) const;
432
438 get_real_point(const unsigned int offset) const;
439
444 get_mapping() const;
445
451
457
461 boost::signals2::connection
462 connect_is_reinitialized(const std::function<void()> &set_is_reinitialized);
463
467 template <bool is_face>
468 unsigned int
470 const unsigned int face_number) const;
471
475 unsigned int
476 compute_unit_point_index_offset(const unsigned int geometry_index) const;
477
481 unsigned int
482 compute_data_index_offset(const unsigned int geometry_index) const;
483
487 unsigned int
489 const unsigned int geometry_index) const;
490
494 unsigned int
495 get_n_q_points_unvectorized(const unsigned int geometry_index) const;
496
501 get_cell_type(const unsigned int geometry_index) const;
502
509 get_cell_iterator(const unsigned int cell_index) const;
510
514 std::size_t
516
517 private:
520 spacedim>;
521
525 template <typename NumberType>
526 unsigned int
528
532 void
534
538 void
539 resize_unit_points(const unsigned int n_unit_point_batches);
540
544 void
545 resize_unit_points_faces(const unsigned int n_unit_point_batches);
546
550 void
551 resize_data_fields(const unsigned int n_data_point_batches,
552 const bool is_face_centric = false);
553
557 void
558 store_unit_points(const unsigned int unit_points_index_offset,
559 const unsigned int n_q_points,
560 const unsigned int n_q_points_unvectorized,
561 const std::vector<Point<dim>> &points);
562
566 void
567 store_unit_points_faces(const unsigned int unit_points_index_offset,
568 const unsigned int n_q_points,
569 const unsigned int n_q_points_unvectorized,
570 const std::vector<Point<dim - 1>> &points);
571
575 void
576 store_mapping_data(const unsigned int unit_points_index_offset,
577 const unsigned int n_q_points,
578 const unsigned int n_q_points_unvectorized,
580 const std::vector<double> &weights,
581 const unsigned int compressed_unit_point_index_offset,
582 const bool affine_cell,
583 const bool is_interior = true);
584
588 unsigned int
589 compute_compressed_cell_index(const unsigned int cell_index) const;
590
594 template <typename ContainerType, typename QuadratureType>
595 void
597 const ContainerType &cell_iterator_range,
598 const std::vector<QuadratureType> &quadrature_vector,
599 const unsigned int n_unfiltered_cells,
600 const std::function<
601 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
602 const QuadratureType &quadrature,
603 MappingData &mapping_data)> &compute_mapping_data);
604
608 enum class State
609 {
610 invalid,
615 };
616
622
629
636
640 std::vector<unsigned int> unit_points_index;
641
645 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
647
652
659
664
669
674
679
686 std::vector<::internal::MatrixFreeFunctions::GeometryType> cell_type;
687
691 std::vector<unsigned int> data_index_offsets;
692
696 std::vector<unsigned int> compressed_data_index_offsets;
697
705
712
719 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
721
729 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
731
738
742 std::vector<unsigned int> n_q_points_unvectorized;
743
748 std::vector<unsigned int> cell_index_offset;
749
754 std::vector<unsigned int> cell_index_to_compressed_cell_index;
755
760
765 boost::signals2::signal<void()> is_reinitialized;
766
773
778 std::vector<std::pair<int, int>> cell_level_and_indices;
779
784 std::vector<std::pair<unsigned char, unsigned char>> face_number;
785 };
786
787 template <int dim, int spacedim, typename Number>
788 inline unsigned int
790 const unsigned int offset,
791 const bool is_interior) const
792 {
793 const auto &face_pair = face_number[offset];
794 return is_interior ? face_pair.first : face_pair.second;
795 }
796
797 // ----------------------- template functions ----------------------
798
799
800 template <int dim, int spacedim, typename Number>
802 const Mapping<dim, spacedim> &mapping,
803 const UpdateFlags update_flags,
804 const AdditionalData additional_data)
805 : mapping(&mapping)
806 , update_flags(update_flags)
807 , update_flags_mapping(update_default)
808 , additional_data(additional_data)
809 {
810 // translate update flags
820
821 // always save quadrature points for now
823
826 this->mapping, update_flags_mapping);
827
828 // construct internal_mapping_data for mappings for reuse in reinit()
829 // calls to avoid frequent memory allocations
831 }
832
833
834
835 template <int dim, int spacedim, typename Number>
836 void
838 {
839 n_q_points_unvectorized.clear();
840 unit_points_index.clear();
841 data_index_offsets.clear();
842 compressed_data_index_offsets.clear();
843 cell_type.clear();
844 }
845
846
847
848 template <int dim, int spacedim, typename Number>
849 void
852 const std::vector<Point<dim>> &unit_points_in)
853 {
854 reinit(cell, make_array_view(unit_points_in));
855 }
856
857
858
859 template <int dim, int spacedim, typename Number>
860 void
863 const ArrayView<const Point<dim>> &unit_points_in)
864 {
865 quadrature.initialize(unit_points_in);
866 reinit(cell, quadrature);
867 }
868
869
870
871 template <int dim, int spacedim, typename Number>
872 void
875 const Quadrature<dim> &quadrature)
876 {
877 n_q_points_unvectorized.resize(1);
878 n_q_points_unvectorized[0] = quadrature.size();
879
880 const unsigned int n_q_points =
881 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
882
883 const unsigned int n_q_points_data =
884 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
885
886 // resize data vectors
887 resize_unit_points(n_q_points);
888 resize_data_fields(n_q_points_data);
889
890 // store unit points
891 store_unit_points(0,
892 n_q_points,
893 n_q_points_unvectorized[0],
894 quadrature.get_points());
895
896 // compute mapping data
900 update_flags_mapping,
901 cell,
902 cell_similarity,
903 quadrature,
904 internal_mapping_data,
905 mapping_data);
906
907 // check for cartesian/affine cell
908 if (!quadrature.empty() &&
909 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
910 {
911 cell_type.push_back(
912 internal::compute_geometry_type(cell->diameter(),
913 mapping_data.inverse_jacobians));
914 }
915 else
916 cell_type.push_back(
918
919 // store mapping data
920 store_mapping_data(
921 0,
922 n_q_points_data,
923 n_q_points_unvectorized[0],
924 mapping_data,
925 quadrature.get_weights(),
926 0,
927 cell_type.back() <=
929
930 unit_points_index.push_back(0);
931 data_index_offsets.push_back(0);
932 compressed_data_index_offsets.push_back(0);
933
934 state = State::single_cell;
935 is_reinitialized();
936 }
937
938
939
940 template <int dim, int spacedim, typename Number>
941 template <typename ContainerType>
942 void
944 const ContainerType &cell_iterator_range,
945 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
946 const unsigned int n_unfiltered_cells)
947 {
948 const unsigned int n_cells = unit_points_vector.size();
949 AssertDimension(n_cells,
950 std::distance(cell_iterator_range.begin(),
951 cell_iterator_range.end()));
952
953 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
954 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
955 quadrature_vector[cell_index] =
956 Quadrature<dim>(unit_points_vector[cell_index]);
957
958 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
959 }
960
961
962
963 template <int dim, int spacedim, typename Number>
964 template <typename ContainerType, typename QuadratureType>
965 void
967 const ContainerType &cell_iterator_range,
968 const std::vector<QuadratureType> &quadrature_vector,
969 const unsigned int n_unfiltered_cells,
970 const std::function<
971 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
972 const QuadratureType &quadrature,
973 MappingData &mapping_data)> &compute_mapping_data)
974 {
975 clear();
976
977 do_cell_index_compression =
978 n_unfiltered_cells != numbers::invalid_unsigned_int;
979
980 const unsigned int n_cells = quadrature_vector.size();
981 AssertDimension(n_cells,
982 std::distance(cell_iterator_range.begin(),
983 cell_iterator_range.end()));
984
985 n_q_points_unvectorized.reserve(n_cells);
986
987 cell_type.reserve(n_cells);
988
989 if (additional_data.store_cells)
990 cell_level_and_indices.resize(n_cells);
991
992 // fill unit points index offset vector
993 unit_points_index.reserve(n_cells + 1);
994 unit_points_index.push_back(0);
995 data_index_offsets.reserve(n_cells + 1);
996 data_index_offsets.push_back(0);
997 for (const auto &quadrature : quadrature_vector)
998 {
999 const unsigned int n_points = quadrature.size();
1000 n_q_points_unvectorized.push_back(n_points);
1001
1002 const unsigned int n_q_points =
1003 compute_n_q_points<VectorizedArrayType>(n_points);
1004 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1005
1006 const unsigned int n_q_points_data =
1007 compute_n_q_points<Number>(n_points);
1008 data_index_offsets.push_back(data_index_offsets.back() +
1009 n_q_points_data);
1010 }
1011
1012 const unsigned int n_unit_points = unit_points_index.back();
1013 const unsigned int n_data_points = data_index_offsets.back();
1014
1015 // resize data vectors
1016 resize_unit_points(n_unit_points);
1017 resize_data_fields(n_data_points);
1018
1019 if (do_cell_index_compression)
1020 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1022
1023 MappingData mapping_data_previous_cell;
1024 unsigned int size_compressed_data = 0;
1025 unsigned int cell_index = 0;
1026 for (const auto &cell : cell_iterator_range)
1027 {
1028 if (additional_data.store_cells)
1029 {
1030 this->triangulation = &cell->get_triangulation();
1031 cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
1032 }
1033
1034 const auto &quadrature = quadrature_vector[cell_index];
1035 const bool empty = quadrature.empty();
1036
1037 // store unit points
1038 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1039 n_q_points_unvectorized[cell_index]);
1040 store_unit_points(unit_points_index[cell_index],
1041 n_q_points,
1042 n_q_points_unvectorized[cell_index],
1043 quadrature.get_points());
1044
1045 // compute mapping data
1046 compute_mapping_data(cell, quadrature, mapping_data);
1047
1048 // store mapping data
1049 const unsigned int n_q_points_data =
1050 compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
1051
1052 // check for cartesian/affine cell
1053 if (!empty &&
1054 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1055 {
1056 cell_type.push_back(
1057 internal::compute_geometry_type(cell->diameter(),
1058 mapping_data.inverse_jacobians));
1059 }
1060 else
1061 cell_type.push_back(
1063
1064 if (cell_index > 0)
1065 {
1066 // check if current and previous cell are affine
1067 const bool affine_cells =
1068 cell_type[cell_index] <=
1070 cell_type[cell_index - 1] <=
1072
1073 // create a comparator to compare inverse Jacobian of current
1074 // and previous cell
1076 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1077 1024.);
1078
1079 // we can only compare if current and previous cell have at least
1080 // one quadrature point and both cells are at least affine
1081 const auto comparison_result =
1082 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1083 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1085 comparator.compare(
1086 mapping_data.inverse_jacobians[0],
1087 mapping_data_previous_cell.inverse_jacobians[0]);
1088
1089 // we can compress the Jacobians and inverse Jacobians if
1090 // inverse Jacobians are equal and cells are affine
1091 if (affine_cells &&
1092 comparison_result ==
1094 {
1095 compressed_data_index_offsets.push_back(
1096 compressed_data_index_offsets.back());
1097 }
1098 else
1099 {
1100 const unsigned int n_compressed_data_last_cell =
1101 cell_type[cell_index - 1] <=
1103 1 :
1104 compute_n_q_points<Number>(
1105 n_q_points_unvectorized[cell_index - 1]);
1106
1107 compressed_data_index_offsets.push_back(
1108 compressed_data_index_offsets.back() +
1109 n_compressed_data_last_cell);
1110 }
1111 }
1112 else
1113 compressed_data_index_offsets.push_back(0);
1114
1115 // cache mapping_data from previous cell
1116 mapping_data_previous_cell = mapping_data;
1117
1118 store_mapping_data(data_index_offsets[cell_index],
1119 n_q_points_data,
1120 n_q_points_unvectorized[cell_index],
1121 mapping_data,
1122 quadrature.get_weights(),
1123 compressed_data_index_offsets[cell_index],
1124 cell_type[cell_index] <=
1126
1127 // update size of compressed data depending on cell type and handle
1128 // empty quadratures
1129 if (cell_type[cell_index] <=
1131 size_compressed_data = compressed_data_index_offsets.back() + 1;
1132 else
1133 size_compressed_data =
1134 std::max(size_compressed_data,
1135 compressed_data_index_offsets.back() + n_q_points_data);
1136
1137 if (do_cell_index_compression)
1138 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1139 cell_index;
1140
1141 ++cell_index;
1142 }
1143
1144 if (update_flags_mapping & UpdateFlags::update_jacobians)
1145 {
1146 jacobians[0].resize(size_compressed_data);
1147 jacobians[0].shrink_to_fit();
1148 }
1149 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1150 {
1151 inverse_jacobians[0].resize(size_compressed_data);
1152 inverse_jacobians[0].shrink_to_fit();
1153 }
1154
1155 state = State::cell_vector;
1156 is_reinitialized();
1157 }
1158
1159
1160
1161 template <int dim, int spacedim, typename Number>
1162 template <typename ContainerType>
1163 void
1165 const ContainerType &cell_iterator_range,
1166 const std::vector<Quadrature<dim>> &quadrature_vector,
1167 const unsigned int n_unfiltered_cells)
1168 {
1169 auto compute_mapping_data_for_cells =
1170 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1171 const Quadrature<dim> &quadrature,
1172 MappingData &mapping_data) {
1173 CellSimilarity::Similarity cell_similarity =
1177 update_flags_mapping,
1178 cell,
1179 cell_similarity,
1180 quadrature,
1181 internal_mapping_data,
1182 mapping_data);
1183 };
1184
1185 do_reinit_cells<ContainerType, Quadrature<dim>>(
1186 cell_iterator_range,
1187 quadrature_vector,
1188 n_unfiltered_cells,
1189 compute_mapping_data_for_cells);
1190 }
1191
1192
1193
1194 template <int dim, int spacedim, typename Number>
1195 template <typename ContainerType>
1196 void
1198 const ContainerType &cell_iterator_range,
1199 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
1200 const unsigned int n_unfiltered_cells)
1201 {
1202 Assert(
1203 additional_data.use_global_weights == false,
1204 ExcMessage(
1205 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1206
1207 Assert(additional_data.store_cells == false, ExcNotImplemented());
1208
1209 if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
1210 update_flags_mapping |= update_covariant_transformation;
1211
1212 auto compute_mapping_data_for_surface =
1213 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1214 const ImmersedSurfaceQuadrature<dim> &quadrature,
1215 MappingData &mapping_data) {
1218 mapping,
1219 update_flags_mapping,
1220 cell,
1221 quadrature,
1222 internal_mapping_data,
1223 mapping_data);
1224 };
1225
1226 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1227 cell_iterator_range,
1228 quadrature_vector,
1229 n_unfiltered_cells,
1230 compute_mapping_data_for_surface);
1231 }
1232
1233
1234
1235 template <int dim, int spacedim, typename Number>
1236 template <typename ContainerType>
1237 void
1239 const ContainerType &cell_iterator_range,
1240 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
1241 const unsigned int n_unfiltered_cells)
1242 {
1243 clear();
1244
1245 Assert(additional_data.store_cells == false, ExcNotImplemented());
1246
1247 do_cell_index_compression =
1248 n_unfiltered_cells != numbers::invalid_unsigned_int;
1249
1250 const unsigned int n_cells = quadrature_vector.size();
1251 AssertDimension(n_cells,
1252 std::distance(cell_iterator_range.begin(),
1253 cell_iterator_range.end()));
1254
1255 // fill cell index offset vector
1256 cell_index_offset.resize(n_cells);
1257 unsigned int n_faces = 0;
1258 unsigned int cell_index = 0;
1259 for (const auto &cell : cell_iterator_range)
1260 {
1261 cell_index_offset[cell_index] = n_faces;
1262 n_faces += cell->n_faces();
1263 ++cell_index;
1264 }
1265
1266 n_q_points_unvectorized.reserve(n_faces);
1267
1268 cell_type.reserve(n_faces);
1269
1270 // fill unit points index offset vector
1271 unit_points_index.resize(n_faces + 1);
1272 data_index_offsets.resize(n_faces + 1);
1273 cell_index = 0;
1274 unsigned int n_unit_points = 0;
1275 unsigned int n_data_points = 0;
1276 for (const auto &cell : cell_iterator_range)
1277 {
1278 for (const auto &f : cell->face_indices())
1279 {
1280 const unsigned int current_face_index =
1281 cell_index_offset[cell_index] + f;
1282
1283 unit_points_index[current_face_index] = n_unit_points;
1284 data_index_offsets[current_face_index] = n_data_points;
1285
1286 const unsigned int n_points =
1287 quadrature_vector[cell_index][f].size();
1288 n_q_points_unvectorized.push_back(n_points);
1289
1290 const unsigned int n_q_points =
1291 compute_n_q_points<VectorizedArrayType>(n_points);
1292 n_unit_points += n_q_points;
1293
1294 const unsigned int n_q_points_data =
1295 compute_n_q_points<Number>(n_points);
1296 n_data_points += n_q_points_data;
1297 }
1298
1299 ++cell_index;
1300 }
1301 unit_points_index[n_faces] = n_unit_points;
1302 data_index_offsets[n_faces] = n_data_points;
1303
1304 // compress indices
1305 if (do_cell_index_compression)
1306 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1308
1309 // fill unit points and mapping data for every face of all cells
1310 // resize data vectors
1311 resize_unit_points_faces(n_unit_points);
1312 resize_data_fields(n_data_points);
1313
1314 MappingData mapping_data_previous_cell;
1315 MappingData mapping_data_first;
1316 bool first_set = false;
1317 unsigned int size_compressed_data = 0;
1318 cell_index = 0;
1319 for (const auto &cell : cell_iterator_range)
1320 {
1321 const auto &quadratures_on_faces = quadrature_vector[cell_index];
1322
1323 Assert(quadratures_on_faces.size() == cell->n_faces(),
1324 ExcDimensionMismatch(quadratures_on_faces.size(),
1325 cell->n_faces()));
1326
1327 for (const auto &f : cell->face_indices())
1328 {
1329 const auto &quadrature_on_face = quadratures_on_faces[f];
1330 const bool empty = quadrature_on_face.empty();
1331
1332 const unsigned int current_face_index =
1333 cell_index_offset[cell_index] + f;
1334
1335 // store unit points
1336 const unsigned int n_q_points =
1337 compute_n_q_points<VectorizedArrayType>(
1338 n_q_points_unvectorized[current_face_index]);
1339 store_unit_points_faces(unit_points_index[current_face_index],
1340 n_q_points,
1341 n_q_points_unvectorized[current_face_index],
1342 quadrature_on_face.get_points());
1343
1346 update_flags_mapping,
1347 cell,
1348 f,
1349 quadrature_on_face,
1350 internal_mapping_data,
1351 mapping_data);
1352
1353 // check for cartesian/affine cell
1354 if (!empty &&
1355 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1356 {
1357 cell_type.push_back(internal::compute_geometry_type(
1358 cell->diameter(), mapping_data.inverse_jacobians));
1359
1360 if (!first_set)
1361 {
1362 mapping_data_first = mapping_data;
1363 first_set = true;
1364 }
1365 }
1366 else
1367 cell_type.push_back(
1369
1370 if (current_face_index > 0)
1371 {
1372 // check if current and previous cell are affine
1373 const bool affine_cells =
1374 cell_type[current_face_index] <=
1376 cell_type[current_face_index - 1] <=
1378
1379 // create a comparator to compare inverse Jacobian of current
1380 // and previous cell
1382 1e4 / cell->diameter() *
1383 std::numeric_limits<double>::epsilon() * 1024.);
1384
1385 // we can only compare if current and previous cell have at
1386 // least one quadrature point and both cells are at least affine
1387 const auto comparison_result =
1388 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1389 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1391 comparator.compare(
1392 mapping_data.inverse_jacobians[0],
1393 mapping_data_previous_cell.inverse_jacobians[0]);
1394
1395 // we can compress the Jacobians and inverse Jacobians if
1396 // inverse Jacobians are equal and cells are affine
1397 if (affine_cells &&
1398 comparison_result ==
1400 {
1401 compressed_data_index_offsets.push_back(
1402 compressed_data_index_offsets.back());
1403 }
1404 else if (first_set &&
1405 (cell_type[current_face_index] <=
1407 (comparator.compare(
1408 mapping_data.inverse_jacobians[0],
1409 mapping_data_first.inverse_jacobians[0]) ==
1411 double>::ComparisonResult::equal))
1412 {
1413 compressed_data_index_offsets.push_back(0);
1414 }
1415 else
1416 {
1417 const unsigned int n_compressed_data_last_cell =
1418 cell_type[current_face_index - 1] <=
1420 1 :
1421 compute_n_q_points<Number>(
1422 n_q_points_unvectorized[current_face_index - 1]);
1423
1424 compressed_data_index_offsets.push_back(
1425 compressed_data_index_offsets.back() +
1426 n_compressed_data_last_cell);
1427 }
1428 }
1429 else
1430 compressed_data_index_offsets.push_back(0);
1431
1432 // cache mapping_data from previous cell
1433 mapping_data_previous_cell = mapping_data;
1434
1435 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1436 n_q_points_unvectorized[current_face_index]);
1437 store_mapping_data(data_index_offsets[current_face_index],
1438 n_q_points_data,
1439 n_q_points_unvectorized[current_face_index],
1440 mapping_data,
1441 quadrature_on_face.get_weights(),
1442 data_index_offsets[current_face_index],
1443 cell_type[current_face_index] <=
1445
1446 // update size of compressed data depending on cell type and handle
1447 // empty quadratures
1448 if (cell_type[current_face_index] <=
1450 size_compressed_data = compressed_data_index_offsets.back() + 1;
1451 else
1452 size_compressed_data =
1453 std::max(size_compressed_data,
1454 compressed_data_index_offsets.back() +
1455 n_q_points_data);
1456 }
1457 if (do_cell_index_compression)
1458 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1459 cell_index;
1460
1461 ++cell_index;
1462 }
1463
1464 if (update_flags_mapping & UpdateFlags::update_jacobians)
1465 {
1466 jacobians[0].resize(size_compressed_data);
1467 jacobians[0].shrink_to_fit();
1468 }
1469 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1470 {
1471 inverse_jacobians[0].resize(size_compressed_data);
1472 inverse_jacobians[0].shrink_to_fit();
1473 }
1474
1475 state = State::faces_on_cells_in_vector;
1476 is_reinitialized();
1477 }
1478
1479
1480
1481 template <int dim, int spacedim, typename Number>
1482 template <typename CellIteratorType>
1483 void
1485 const std::vector<std::pair<CellIteratorType, unsigned int>>
1486 &face_iterator_range_interior,
1487 const std::vector<Quadrature<dim - 1>> &quadrature_vector)
1488 {
1489 clear();
1490
1491 do_cell_index_compression = false;
1492
1493 Assert(additional_data.store_cells == false, ExcNotImplemented());
1494
1495
1496 const unsigned int n_faces = quadrature_vector.size();
1497 AssertDimension(n_faces,
1498 std::distance(face_iterator_range_interior.begin(),
1499 face_iterator_range_interior.end()));
1500
1501 n_q_points_unvectorized.reserve(n_faces);
1502
1503 cell_type.reserve(n_faces);
1504 face_number.reserve(n_faces);
1505
1506 // fill unit points index offset vector
1507 unit_points_index.reserve(n_faces + 1);
1508 unit_points_index.push_back(0);
1509 data_index_offsets.reserve(n_faces + 1);
1510 data_index_offsets.push_back(0);
1511 for (const auto &quadrature : quadrature_vector)
1512 {
1513 const unsigned int n_points = quadrature.size();
1514 n_q_points_unvectorized.push_back(n_points);
1515
1516 const unsigned int n_q_points =
1517 compute_n_q_points<VectorizedArrayType>(n_points);
1518 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1519
1520 const unsigned int n_q_points_data =
1521 compute_n_q_points<Number>(n_points);
1522 data_index_offsets.push_back(data_index_offsets.back() +
1523 n_q_points_data);
1524 }
1525
1526 const unsigned int n_unit_points = unit_points_index.back();
1527 const unsigned int n_data_points = data_index_offsets.back();
1528
1529 // resize data vectors
1530 resize_unit_points_faces(n_unit_points);
1531 resize_data_fields(n_data_points, true);
1532
1533 std::array<MappingData, 2> mapping_data;
1534 std::array<MappingData, 2> mapping_data_previous_cell;
1535 std::array<MappingData, 2> mapping_data_first;
1536 bool first_set = false;
1537 unsigned int size_compressed_data = 0;
1538 unsigned int face_index = 0;
1539 for (const auto &cell_and_f : face_iterator_range_interior)
1540 {
1541 const auto &quadrature_on_face = quadrature_vector[face_index];
1542 const bool empty = quadrature_on_face.empty();
1543
1544 // get interior cell and face number
1545 const auto &cell_m = cell_and_f.first;
1546 const auto f_m = cell_and_f.second;
1547
1548 // get exterior cell and face number
1549 const auto &cell_p =
1550 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1551 const auto f_p =
1552 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1553
1554 face_number.emplace_back(f_m, f_p);
1555
1556 Assert(
1557 cell_m->combined_face_orientation(f_m) ==
1559 cell_p->combined_face_orientation(f_p) ==
1561 ExcMessage(
1562 "Non standard face orientation is currently not implemented."));
1563
1564 // store unit points
1565 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1566 n_q_points_unvectorized[face_index]);
1567 store_unit_points_faces(unit_points_index[face_index],
1568 n_q_points,
1569 n_q_points_unvectorized[face_index],
1570 quadrature_on_face.get_points());
1571
1572 // compute mapping for interior face
1575 update_flags_mapping,
1576 cell_m,
1577 f_m,
1578 quadrature_on_face,
1579 internal_mapping_data,
1580 mapping_data[0]);
1581
1582 // compute mapping for exterior face
1585 update_flags_mapping,
1586 cell_p,
1587 f_p,
1588 quadrature_on_face,
1589 internal_mapping_data,
1590 mapping_data[1]);
1591
1592 // check for cartesian/affine cell
1593 if (!empty &&
1594 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1595 {
1596 // select more general type of interior and exterior cell
1597 cell_type.push_back(std::max(
1599 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1601 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1602
1603 // cache mapping data of first cell pair with non-empty quadrature
1604 // on the face
1605 if (!first_set)
1606 {
1607 mapping_data_first = mapping_data;
1608 first_set = true;
1609 }
1610 }
1611 else
1612 cell_type.push_back(
1614
1615 if (face_index > 0)
1616 {
1617 // check if current and previous cell pairs are affine
1618 const bool affine_cells =
1619 cell_type[face_index] <=
1621 cell_type[face_index - 1] <=
1623
1624 // create a comparator to compare inverse Jacobian of current
1625 // and previous cell pair
1627 1e4 / cell_m->diameter() *
1628 std::numeric_limits<double>::epsilon() * 1024.);
1629
1630 // we can only compare if current and previous cell have at
1631 // least one quadrature point and both cells are at least affine
1632 const auto comparison_result_m =
1633 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1634 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1636 comparator.compare(
1637 mapping_data[0].inverse_jacobians[0],
1638 mapping_data_previous_cell[0].inverse_jacobians[0]);
1639
1640 const auto comparison_result_p =
1641 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1642 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1644 comparator.compare(
1645 mapping_data[1].inverse_jacobians[0],
1646 mapping_data_previous_cell[1].inverse_jacobians[0]);
1647
1648 // we can compress the Jacobians and inverse Jacobians if
1649 // inverse Jacobians are equal and cells are affine
1650 if (affine_cells &&
1651 comparison_result_m ==
1653 comparison_result_p ==
1655 {
1656 compressed_data_index_offsets.push_back(
1657 compressed_data_index_offsets.back());
1658 }
1659 else if (first_set &&
1660 (cell_type[face_index] <=
1662 (comparator.compare(
1663 mapping_data[0].inverse_jacobians[0],
1664 mapping_data_first[0].inverse_jacobians[0]) ==
1666 double>::ComparisonResult::equal) &&
1667 (comparator.compare(
1668 mapping_data[1].inverse_jacobians[0],
1669 mapping_data_first[1].inverse_jacobians[0]) ==
1671 {
1672 compressed_data_index_offsets.push_back(0);
1673 }
1674 else
1675 {
1676 const unsigned int n_compressed_data_last_cell =
1677 cell_type[face_index - 1] <=
1679 1 :
1680 compute_n_q_points<Number>(
1681 n_q_points_unvectorized[face_index - 1]);
1682
1683 compressed_data_index_offsets.push_back(
1684 compressed_data_index_offsets.back() +
1685 n_compressed_data_last_cell);
1686 }
1687 }
1688 else
1689 compressed_data_index_offsets.push_back(0);
1690
1691 // cache mapping_data from previous cell pair
1692 mapping_data_previous_cell = mapping_data;
1693
1694 const unsigned int n_q_points_data =
1695 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1696
1697 // store mapping data of interior face
1698 store_mapping_data(data_index_offsets[face_index],
1699 n_q_points_data,
1700 n_q_points_unvectorized[face_index],
1701 mapping_data[0],
1702 quadrature_on_face.get_weights(),
1703 data_index_offsets[face_index],
1704 cell_type[face_index] <=
1706 true);
1707
1708 // store only necessary mapping data for exterior face (Jacobians and
1709 // inverse Jacobians)
1710 store_mapping_data(data_index_offsets[face_index],
1711 n_q_points_data,
1712 n_q_points_unvectorized[face_index],
1713 mapping_data[1],
1714 quadrature_on_face.get_weights(),
1715 data_index_offsets[face_index],
1716 cell_type[face_index] <=
1718 false);
1719
1720 // update size of compressed data depending on cell type and handle
1721 // empty quadratures
1722 if (cell_type[face_index] <=
1724 size_compressed_data = compressed_data_index_offsets.back() + 1;
1725 else
1726 size_compressed_data =
1727 std::max(size_compressed_data,
1728 compressed_data_index_offsets.back() + n_q_points_data);
1729
1730 ++face_index;
1731 }
1732
1733 if (update_flags_mapping & UpdateFlags::update_jacobians)
1734 {
1735 jacobians[0].resize(size_compressed_data);
1736 jacobians[0].shrink_to_fit();
1737 jacobians[1].resize(size_compressed_data);
1738 jacobians[1].shrink_to_fit();
1739 }
1740 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1741 {
1742 inverse_jacobians[0].resize(size_compressed_data);
1743 inverse_jacobians[0].shrink_to_fit();
1744 inverse_jacobians[1].resize(size_compressed_data);
1745 inverse_jacobians[1].shrink_to_fit();
1746 }
1747
1748 state = State::face_vector;
1749 is_reinitialized();
1750 }
1751
1752
1753
1754 template <int dim, int spacedim, typename Number>
1755 bool
1757 {
1758 return state == State::faces_on_cells_in_vector;
1759 }
1760
1761
1762
1763 template <int dim, int spacedim, typename Number>
1764 unsigned int
1766 const unsigned int geometry_index) const
1767 {
1768 return n_q_points_unvectorized[geometry_index];
1769 }
1770
1771
1772 template <int dim, int spacedim, typename Number>
1775 const unsigned int geometry_index) const
1776 {
1777 return cell_type[geometry_index];
1778 }
1779
1780
1781
1782 template <int dim, int spacedim, typename Number>
1785 const unsigned int cell_index) const
1786 {
1787 Assert(
1788 additional_data.store_cells,
1789 ExcMessage(
1790 "Cells have been not stored. You can enable this by Additional::store_cells."));
1791 return {triangulation.get(),
1792 cell_level_and_indices[cell_index].first,
1793 cell_level_and_indices[cell_index].second};
1794 }
1795
1796
1797
1798 template <int dim, int spacedim, typename Number>
1799 template <typename NumberType>
1800 unsigned int
1802 const unsigned int n_q_points_unvectorized)
1803 {
1804 const unsigned int n_lanes =
1806 const unsigned int n_filled_lanes_last_batch =
1807 n_q_points_unvectorized % n_lanes;
1808 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1809 if (n_filled_lanes_last_batch > 0)
1810 ++n_q_points;
1811 return n_q_points;
1812 }
1813
1814
1815
1816 template <int dim, int spacedim, typename Number>
1817 template <bool is_face>
1818 unsigned int
1820 const unsigned int cell_index,
1821 const unsigned int face_number) const
1822 {
1824 return 0;
1825
1826 const unsigned int compressed_cell_index =
1827 compute_compressed_cell_index(cell_index);
1828 if (!is_face)
1829 {
1830 Assert(state == State::cell_vector,
1831 ExcMessage(
1832 "This mapping info is not reinitialized for a cell vector!"));
1833 return compressed_cell_index;
1834 }
1835 else
1836 {
1838 ExcMessage(
1839 "cell_index has to be set if face_number is specified!"));
1840 Assert(state == State::faces_on_cells_in_vector ||
1841 state == State::face_vector,
1842 ExcMessage("This mapping info is not reinitialized for faces"
1843 " on cells in a vector!"));
1844 if (state == State::faces_on_cells_in_vector)
1845 return cell_index_offset[compressed_cell_index] + face_number;
1846 else if (state == State::face_vector)
1847 return cell_index;
1848 }
1849 }
1850
1851
1852
1853 template <int dim, int spacedim, typename Number>
1854 unsigned int
1856 const unsigned int cell_index) const
1857 {
1858 if (do_cell_index_compression)
1859 {
1860 Assert(cell_index_to_compressed_cell_index[cell_index] !=
1862 ExcMessage("Mapping info object was not initialized for this"
1863 " active cell index!"));
1864 return cell_index_to_compressed_cell_index[cell_index];
1865 }
1866 else
1867 return cell_index;
1868 }
1869
1870
1871 template <int dim, int spacedim, typename Number>
1872 void
1874 const unsigned int unit_points_index_offset,
1875 const unsigned int n_q_points,
1876 const unsigned int n_q_points_unvectorized,
1877 const std::vector<Point<dim>> &points)
1878 {
1879 const unsigned int n_lanes =
1881
1882 for (unsigned int q = 0; q < n_q_points; ++q)
1883 {
1884 const unsigned int offset = unit_points_index_offset + q;
1885 for (unsigned int v = 0;
1886 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1887 ++v)
1888 for (unsigned int d = 0; d < dim; ++d)
1890 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1891 }
1892 }
1893
1894
1895
1896 template <int dim, int spacedim, typename Number>
1897 void
1899 const unsigned int unit_points_index_offset,
1900 const unsigned int n_q_points,
1901 const unsigned int n_q_points_unvectorized,
1902 const std::vector<Point<dim - 1>> &points)
1903 {
1904 const unsigned int n_lanes =
1906
1907 for (unsigned int q = 0; q < n_q_points; ++q)
1908 {
1909 const unsigned int offset = unit_points_index_offset + q;
1910 for (unsigned int v = 0;
1911 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1912 ++v)
1913 for (unsigned int d = 0; d < dim - 1; ++d)
1915 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1916 }
1917 }
1918
1919
1920
1921 template <int dim, int spacedim, typename Number>
1922 void
1924 const unsigned int unit_points_index_offset,
1925 const unsigned int n_q_points,
1926 const unsigned int n_q_points_unvectorized,
1927 const MappingInfo::MappingData &mapping_data,
1928 const std::vector<double> &weights,
1929 const unsigned int compressed_unit_point_index_offset,
1930 const bool affine_cell,
1931 const bool is_interior)
1932 {
1933 const unsigned int n_lanes =
1935
1936 for (unsigned int q = 0; q < n_q_points; ++q)
1937 {
1938 const unsigned int offset = unit_points_index_offset + q;
1939 const unsigned int compressed_offset =
1940 compressed_unit_point_index_offset + q;
1941 for (unsigned int v = 0;
1942 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1943 ++v)
1944 {
1945 if (q == 0 || !affine_cell)
1946 {
1947 if (update_flags_mapping & UpdateFlags::update_jacobians)
1948 for (unsigned int d = 0; d < dim; ++d)
1949 for (unsigned int s = 0; s < spacedim; ++s)
1951 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1952 v) = mapping_data.jacobians[q * n_lanes + v][d][s];
1953 if (update_flags_mapping &
1955 for (unsigned int d = 0; d < dim; ++d)
1956 for (unsigned int s = 0; s < spacedim; ++s)
1958 inverse_jacobians[is_interior ? 0 : 1]
1959 [compressed_offset][d][s],
1960 v) =
1961 mapping_data.inverse_jacobians[q * n_lanes + v][d][s];
1962 }
1963
1964 if (is_interior)
1965 {
1966 if (update_flags_mapping & UpdateFlags::update_JxW_values)
1967 {
1968 if (additional_data.use_global_weights)
1969 {
1971 JxW_values[offset], v) = weights[q * n_lanes + v];
1972 }
1973 else
1974 {
1976 JxW_values[offset], v) =
1977 mapping_data.JxW_values[q * n_lanes + v];
1978 }
1979 }
1980 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
1981 for (unsigned int s = 0; s < spacedim; ++s)
1983 normal_vectors[offset][s], v) =
1984 mapping_data.normal_vectors[q * n_lanes + v][s];
1985 if (update_flags_mapping &
1987 for (unsigned int s = 0; s < spacedim; ++s)
1989 real_points[offset][s], v) =
1990 mapping_data.quadrature_points[q * n_lanes + v][s];
1991 }
1992 }
1993 }
1994 }
1995
1996
1997
1998 template <int dim, int spacedim, typename Number>
1999 void
2001 const unsigned int n_unit_point_batches)
2002 {
2003 unit_points.resize(n_unit_point_batches);
2004 }
2005
2006
2007
2008 template <int dim, int spacedim, typename Number>
2009 void
2011 const unsigned int n_unit_point_batches)
2012 {
2013 unit_points_faces.resize(n_unit_point_batches);
2014 }
2015
2016
2017
2018 template <int dim, int spacedim, typename Number>
2019 void
2021 const unsigned int n_data_point_batches,
2022 const bool is_face_centric)
2023 {
2024 if (update_flags_mapping & UpdateFlags::update_jacobians)
2025 {
2026 jacobians[0].resize(n_data_point_batches);
2027 if (is_face_centric)
2028 jacobians[1].resize(n_data_point_batches);
2029 }
2030 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2031 {
2032 inverse_jacobians[0].resize(n_data_point_batches);
2033 if (is_face_centric)
2034 inverse_jacobians[1].resize(n_data_point_batches);
2035 }
2036 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2037 JxW_values.resize(n_data_point_batches);
2038 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2039 normal_vectors.resize(n_data_point_batches);
2040 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2041 real_points.resize(n_data_point_batches);
2042 }
2043
2044
2045
2046 template <int dim, int spacedim, typename Number>
2047 inline const Point<
2048 dim,
2051 const unsigned int offset) const
2052 {
2053 return unit_points.data() + offset;
2054 }
2055
2056
2057
2058 template <int dim, int spacedim, typename Number>
2059 inline const Point<
2060 dim - 1,
2063 const unsigned int offset) const
2064 {
2065 return unit_points_faces.data() + offset;
2066 }
2067
2068
2069
2070 template <int dim, int spacedim, typename Number>
2071 inline const Point<spacedim, Number> *
2073 const unsigned int offset) const
2074 {
2075 return real_points.data() + offset;
2076 }
2077
2078
2079
2080 template <int dim, int spacedim, typename Number>
2081 unsigned int
2083 const unsigned int geometry_index) const
2084 {
2085 return unit_points_index[geometry_index];
2086 }
2087
2088
2089
2090 template <int dim, int spacedim, typename Number>
2091 unsigned int
2093 const unsigned int geometry_index) const
2094 {
2095 return data_index_offsets[geometry_index];
2096 }
2097
2098
2099 template <int dim, int spacedim, typename Number>
2100 unsigned int
2102 const unsigned int geometry_index) const
2103 {
2104 return compressed_data_index_offsets[geometry_index];
2105 }
2106
2107
2108
2109 template <int dim, int spacedim, typename Number>
2112 const bool is_interior) const
2113 {
2114 return jacobians[is_interior ? 0 : 1].data() + offset;
2115 }
2116
2117
2118
2119 template <int dim, int spacedim, typename Number>
2122 const unsigned int offset,
2123 const bool is_interior) const
2124 {
2125 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2126 }
2127
2128
2129
2130 template <int dim, int spacedim, typename Number>
2131 inline const Tensor<1, spacedim, Number> *
2133 const unsigned int offset) const
2134 {
2135 return normal_vectors.data() + offset;
2136 }
2137
2138
2139
2140 template <int dim, int spacedim, typename Number>
2141 inline const Number *
2142 MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
2143 {
2144 return JxW_values.data() + offset;
2145 }
2146
2147
2148
2149 template <int dim, int spacedim, typename Number>
2152 {
2153 return *mapping;
2154 }
2155
2156
2157
2158 template <int dim, int spacedim, typename Number>
2161 {
2162 return update_flags;
2163 }
2164
2165
2166
2167 template <int dim, int spacedim, typename Number>
2170 {
2171 return update_flags_mapping;
2172 }
2173
2174
2175
2176 template <int dim, int spacedim, typename Number>
2177 boost::signals2::connection
2179 const std::function<void()> &set_is_reinitialized)
2180 {
2181 return is_reinitialized.connect(set_is_reinitialized);
2182 }
2183
2184
2185
2186 template <int dim, int spacedim, typename Number>
2187 std::size_t
2189 {
2190 std::size_t memory = MemoryConsumption::memory_consumption(unit_points);
2191 memory += MemoryConsumption::memory_consumption(unit_points_faces);
2192 memory += MemoryConsumption::memory_consumption(unit_points_index);
2193 memory += cell_type.capacity() *
2195 memory += MemoryConsumption::memory_consumption(data_index_offsets);
2196 memory +=
2197 MemoryConsumption::memory_consumption(compressed_data_index_offsets);
2198 memory += MemoryConsumption::memory_consumption(JxW_values);
2199 memory += MemoryConsumption::memory_consumption(normal_vectors);
2200 memory += MemoryConsumption::memory_consumption(jacobians);
2201 memory += MemoryConsumption::memory_consumption(inverse_jacobians);
2202 memory += MemoryConsumption::memory_consumption(real_points);
2203 memory += MemoryConsumption::memory_consumption(n_q_points_unvectorized);
2204 memory += MemoryConsumption::memory_consumption(cell_index_offset);
2206 cell_index_to_compressed_cell_index);
2207 memory += MemoryConsumption::memory_consumption(cell_level_and_indices);
2208 memory += sizeof(*this);
2209 return memory;
2210 }
2211} // namespace NonMatching
2212
2214
2215#endif
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
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:318
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
boost::signals2::connection connect_is_reinitialized(const std::function< void()> &set_is_reinitialized)
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
SmartPointer< const Triangulation< dim, spacedim > > triangulation
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
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)
const SmartPointer< const Mapping< dim, spacedim > > mapping
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
boost::signals2::signal< void()> is_reinitialized
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 UpdateFlags required_update_flags(const SmartPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags)
static void compute_mapping_data_for_immersed_surface_quadrature(const SmartPointer< 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_quadrature(const SmartPointer< 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 void compute_mapping_data_for_face_quadrature(const SmartPointer< 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)
Definition point.h:111
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
bool empty() const
unsigned int size() const
static constexpr unsigned char default_combined_face_orientation()
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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::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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::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)