Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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 template <bool is_face>
462 unsigned int
464 const unsigned int face_number) const;
465
469 unsigned int
470 compute_unit_point_index_offset(const unsigned int geometry_index) const;
471
475 unsigned int
476 compute_data_index_offset(const unsigned int geometry_index) const;
477
481 unsigned int
483 const unsigned int geometry_index) const;
484
488 unsigned int
489 get_n_q_points_unvectorized(const unsigned int geometry_index) const;
490
495 get_cell_type(const unsigned int geometry_index) const;
496
503 get_cell_iterator(const unsigned int cell_index) const;
504
508 std::size_t
510
511 private:
514 spacedim>;
515
519 template <typename NumberType>
520 unsigned int
522
526 void
528
532 void
533 resize_unit_points(const unsigned int n_unit_point_batches);
534
538 void
539 resize_unit_points_faces(const unsigned int n_unit_point_batches);
540
544 void
545 resize_data_fields(const unsigned int n_data_point_batches,
546 const bool is_face_centric = false);
547
551 void
552 store_unit_points(const unsigned int unit_points_index_offset,
553 const unsigned int n_q_points,
554 const unsigned int n_q_points_unvectorized,
555 const std::vector<Point<dim>> &points);
556
560 void
561 store_unit_points_faces(const unsigned int unit_points_index_offset,
562 const unsigned int n_q_points,
563 const unsigned int n_q_points_unvectorized,
564 const std::vector<Point<dim - 1>> &points);
565
569 void
570 store_mapping_data(const unsigned int unit_points_index_offset,
571 const unsigned int n_q_points,
572 const unsigned int n_q_points_unvectorized,
574 const std::vector<double> &weights,
575 const unsigned int compressed_unit_point_index_offset,
576 const bool affine_cell,
577 const bool is_interior = true);
578
582 unsigned int
583 compute_compressed_cell_index(const unsigned int cell_index) const;
584
588 template <typename ContainerType, typename QuadratureType>
589 void
591 const ContainerType &cell_iterator_range,
592 const std::vector<QuadratureType> &quadrature_vector,
593 const unsigned int n_unfiltered_cells,
594 const std::function<
595 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
596 const QuadratureType &quadrature,
597 MappingData &mapping_data)> &compute_mapping_data);
598
602 enum class State
603 {
604 invalid,
609 };
610
616
623
630
634 std::vector<unsigned int> unit_points_index;
635
639 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
641
646
653
658
663
668
673
680 std::vector<::internal::MatrixFreeFunctions::GeometryType> cell_type;
681
685 std::vector<unsigned int> data_index_offsets;
686
690 std::vector<unsigned int> compressed_data_index_offsets;
691
699
706
713 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
715
723 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
725
732
736 std::vector<unsigned int> n_q_points_unvectorized;
737
742 std::vector<unsigned int> cell_index_offset;
743
748 std::vector<unsigned int> cell_index_to_compressed_cell_index;
749
754
761
766 std::vector<std::pair<int, int>> cell_level_and_indices;
767
772 std::vector<std::pair<unsigned char, unsigned char>> face_number;
773 };
774
775 template <int dim, int spacedim, typename Number>
776 inline unsigned int
778 const unsigned int offset,
779 const bool is_interior) const
780 {
781 const auto &face_pair = face_number[offset];
782 return is_interior ? face_pair.first : face_pair.second;
783 }
784
785 // ----------------------- template functions ----------------------
786
787
788 template <int dim, int spacedim, typename Number>
790 const Mapping<dim, spacedim> &mapping,
791 const UpdateFlags update_flags,
792 const AdditionalData additional_data)
793 : mapping(&mapping)
794 , update_flags(update_flags)
795 , update_flags_mapping(update_default)
796 , additional_data(additional_data)
797 {
798 // translate update flags
808
809 // always save quadrature points for now
811
814 this->mapping, update_flags_mapping);
815
816 // construct internal_mapping_data for mappings for reuse in reinit()
817 // calls to avoid frequent memory allocations
819 }
820
821
822
823 template <int dim, int spacedim, typename Number>
824 void
826 {
827 n_q_points_unvectorized.clear();
828 unit_points_index.clear();
829 data_index_offsets.clear();
830 compressed_data_index_offsets.clear();
831 cell_type.clear();
832 }
833
834
835
836 template <int dim, int spacedim, typename Number>
837 void
840 const std::vector<Point<dim>> &unit_points_in)
841 {
842 reinit(cell, make_array_view(unit_points_in));
843 }
844
845
846
847 template <int dim, int spacedim, typename Number>
848 void
851 const ArrayView<const Point<dim>> &unit_points_in)
852 {
853 quadrature.initialize(unit_points_in);
854 reinit(cell, quadrature);
855 }
856
857
858
859 template <int dim, int spacedim, typename Number>
860 void
863 const Quadrature<dim> &quadrature)
864 {
865 n_q_points_unvectorized.resize(1);
866 n_q_points_unvectorized[0] = quadrature.size();
867
868 const unsigned int n_q_points =
869 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
870
871 const unsigned int n_q_points_data =
872 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
873
874 // resize data vectors
875 resize_unit_points(n_q_points);
876 resize_data_fields(n_q_points_data);
877
878 // store unit points
879 store_unit_points(0,
880 n_q_points,
881 n_q_points_unvectorized[0],
882 quadrature.get_points());
883
884 // compute mapping data
888 update_flags_mapping,
889 cell,
890 cell_similarity,
891 quadrature,
892 internal_mapping_data,
893 mapping_data);
894
895 // check for cartesian/affine cell
896 if (!quadrature.empty() &&
897 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
898 {
899 cell_type.push_back(
900 internal::compute_geometry_type(cell->diameter(),
901 mapping_data.inverse_jacobians));
902 }
903 else
904 cell_type.push_back(
906
907 // store mapping data
908 store_mapping_data(
909 0,
910 n_q_points_data,
911 n_q_points_unvectorized[0],
912 mapping_data,
913 quadrature.get_weights(),
914 0,
915 cell_type.back() <=
917
918 unit_points_index.push_back(0);
919 data_index_offsets.push_back(0);
920 compressed_data_index_offsets.push_back(0);
921
922 state = State::single_cell;
923 }
924
925
926
927 template <int dim, int spacedim, typename Number>
928 template <typename ContainerType>
929 void
931 const ContainerType &cell_iterator_range,
932 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
933 const unsigned int n_unfiltered_cells)
934 {
935 const unsigned int n_cells = unit_points_vector.size();
936 AssertDimension(n_cells,
937 std::distance(cell_iterator_range.begin(),
938 cell_iterator_range.end()));
939
940 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
941 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
942 quadrature_vector[cell_index] =
943 Quadrature<dim>(unit_points_vector[cell_index]);
944
945 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
946 }
947
948
949
950 template <int dim, int spacedim, typename Number>
951 template <typename ContainerType, typename QuadratureType>
952 void
954 const ContainerType &cell_iterator_range,
955 const std::vector<QuadratureType> &quadrature_vector,
956 const unsigned int n_unfiltered_cells,
957 const std::function<
958 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
959 const QuadratureType &quadrature,
960 MappingData &mapping_data)> &compute_mapping_data)
961 {
962 clear();
963
964 do_cell_index_compression =
965 n_unfiltered_cells != numbers::invalid_unsigned_int;
966
967 const unsigned int n_cells = quadrature_vector.size();
968 AssertDimension(n_cells,
969 std::distance(cell_iterator_range.begin(),
970 cell_iterator_range.end()));
971
972 n_q_points_unvectorized.reserve(n_cells);
973
974 cell_type.reserve(n_cells);
975
976 if (additional_data.store_cells)
977 cell_level_and_indices.resize(n_cells);
978
979 // fill unit points index offset vector
980 unit_points_index.reserve(n_cells + 1);
981 unit_points_index.push_back(0);
982 data_index_offsets.reserve(n_cells + 1);
983 data_index_offsets.push_back(0);
984 for (const auto &quadrature : quadrature_vector)
985 {
986 const unsigned int n_points = quadrature.size();
987 n_q_points_unvectorized.push_back(n_points);
988
989 const unsigned int n_q_points =
990 compute_n_q_points<VectorizedArrayType>(n_points);
991 unit_points_index.push_back(unit_points_index.back() + n_q_points);
992
993 const unsigned int n_q_points_data =
994 compute_n_q_points<Number>(n_points);
995 data_index_offsets.push_back(data_index_offsets.back() +
996 n_q_points_data);
997 }
998
999 const unsigned int n_unit_points = unit_points_index.back();
1000 const unsigned int n_data_points = data_index_offsets.back();
1001
1002 // resize data vectors
1003 resize_unit_points(n_unit_points);
1004 resize_data_fields(n_data_points);
1005
1006 if (do_cell_index_compression)
1007 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1009
1010 MappingData mapping_data_previous_cell;
1011 unsigned int size_compressed_data = 0;
1012 unsigned int cell_index = 0;
1013 for (const auto &cell : cell_iterator_range)
1014 {
1015 if (additional_data.store_cells)
1016 {
1017 this->triangulation = &cell->get_triangulation();
1018 cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
1019 }
1020
1021 const auto &quadrature = quadrature_vector[cell_index];
1022 const bool empty = quadrature.empty();
1023
1024 // store unit points
1025 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1026 n_q_points_unvectorized[cell_index]);
1027 store_unit_points(unit_points_index[cell_index],
1028 n_q_points,
1029 n_q_points_unvectorized[cell_index],
1030 quadrature.get_points());
1031
1032 // compute mapping data
1033 compute_mapping_data(cell, quadrature, mapping_data);
1034
1035 // store mapping data
1036 const unsigned int n_q_points_data =
1037 compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
1038
1039 // check for cartesian/affine cell
1040 if (!empty &&
1041 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1042 {
1043 cell_type.push_back(
1044 internal::compute_geometry_type(cell->diameter(),
1045 mapping_data.inverse_jacobians));
1046 }
1047 else
1048 cell_type.push_back(
1050
1051 if (cell_index > 0)
1052 {
1053 // check if current and previous cell are affine
1054 const bool affine_cells =
1055 cell_type[cell_index] <=
1057 cell_type[cell_index - 1] <=
1059
1060 // create a comparator to compare inverse Jacobian of current
1061 // and previous cell
1063 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1064 1024.);
1065
1066 // we can only compare if current and previous cell have at least
1067 // one quadrature point and both cells are at least affine
1068 const auto comparison_result =
1069 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1070 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1072 comparator.compare(
1073 mapping_data.inverse_jacobians[0],
1074 mapping_data_previous_cell.inverse_jacobians[0]);
1075
1076 // we can compress the Jacobians and inverse Jacobians if
1077 // inverse Jacobians are equal and cells are affine
1078 if (affine_cells &&
1079 comparison_result ==
1081 {
1082 compressed_data_index_offsets.push_back(
1083 compressed_data_index_offsets.back());
1084 }
1085 else
1086 {
1087 const unsigned int n_compressed_data_last_cell =
1088 cell_type[cell_index - 1] <=
1090 1 :
1091 compute_n_q_points<Number>(
1092 n_q_points_unvectorized[cell_index - 1]);
1093
1094 compressed_data_index_offsets.push_back(
1095 compressed_data_index_offsets.back() +
1096 n_compressed_data_last_cell);
1097 }
1098 }
1099 else
1100 compressed_data_index_offsets.push_back(0);
1101
1102 // cache mapping_data from previous cell
1103 mapping_data_previous_cell = mapping_data;
1104
1105 store_mapping_data(data_index_offsets[cell_index],
1106 n_q_points_data,
1107 n_q_points_unvectorized[cell_index],
1108 mapping_data,
1109 quadrature.get_weights(),
1110 compressed_data_index_offsets[cell_index],
1111 cell_type[cell_index] <=
1113
1114 // update size of compressed data depending on cell type and handle
1115 // empty quadratures
1116 if (cell_type[cell_index] <=
1118 size_compressed_data = compressed_data_index_offsets.back() + 1;
1119 else
1120 size_compressed_data =
1121 std::max(size_compressed_data,
1122 compressed_data_index_offsets.back() + n_q_points_data);
1123
1124 if (do_cell_index_compression)
1125 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1126 cell_index;
1127
1128 ++cell_index;
1129 }
1130
1131 if (update_flags_mapping & UpdateFlags::update_jacobians)
1132 {
1133 jacobians[0].resize(size_compressed_data);
1134 jacobians[0].shrink_to_fit();
1135 }
1136 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1137 {
1138 inverse_jacobians[0].resize(size_compressed_data);
1139 inverse_jacobians[0].shrink_to_fit();
1140 }
1141
1142 state = State::cell_vector;
1143 }
1144
1145
1146
1147 template <int dim, int spacedim, typename Number>
1148 template <typename ContainerType>
1149 void
1151 const ContainerType &cell_iterator_range,
1152 const std::vector<Quadrature<dim>> &quadrature_vector,
1153 const unsigned int n_unfiltered_cells)
1154 {
1155 auto compute_mapping_data_for_cells =
1156 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1157 const Quadrature<dim> &quadrature,
1158 MappingData &mapping_data) {
1159 CellSimilarity::Similarity cell_similarity =
1163 update_flags_mapping,
1164 cell,
1165 cell_similarity,
1166 quadrature,
1167 internal_mapping_data,
1168 mapping_data);
1169 };
1170
1171 do_reinit_cells<ContainerType, Quadrature<dim>>(
1172 cell_iterator_range,
1173 quadrature_vector,
1174 n_unfiltered_cells,
1175 compute_mapping_data_for_cells);
1176 }
1177
1178
1179
1180 template <int dim, int spacedim, typename Number>
1181 template <typename ContainerType>
1182 void
1184 const ContainerType &cell_iterator_range,
1185 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
1186 const unsigned int n_unfiltered_cells)
1187 {
1188 Assert(
1189 additional_data.use_global_weights == false,
1190 ExcMessage(
1191 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1192
1193 Assert(additional_data.store_cells == false, ExcNotImplemented());
1194
1195 if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
1196 update_flags_mapping |= update_covariant_transformation;
1197
1198 auto compute_mapping_data_for_surface =
1199 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1200 const ImmersedSurfaceQuadrature<dim> &quadrature,
1201 MappingData &mapping_data) {
1204 mapping,
1205 update_flags_mapping,
1206 cell,
1207 quadrature,
1208 internal_mapping_data,
1209 mapping_data);
1210 };
1211
1212 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1213 cell_iterator_range,
1214 quadrature_vector,
1215 n_unfiltered_cells,
1216 compute_mapping_data_for_surface);
1217 }
1218
1219
1220
1221 template <int dim, int spacedim, typename Number>
1222 template <typename ContainerType>
1223 void
1225 const ContainerType &cell_iterator_range,
1226 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
1227 const unsigned int n_unfiltered_cells)
1228 {
1229 clear();
1230
1231 Assert(additional_data.store_cells == false, ExcNotImplemented());
1232
1233 do_cell_index_compression =
1234 n_unfiltered_cells != numbers::invalid_unsigned_int;
1235
1236 const unsigned int n_cells = quadrature_vector.size();
1237 AssertDimension(n_cells,
1238 std::distance(cell_iterator_range.begin(),
1239 cell_iterator_range.end()));
1240
1241 // fill cell index offset vector
1242 cell_index_offset.resize(n_cells);
1243 unsigned int n_faces = 0;
1244 unsigned int cell_index = 0;
1245 for (const auto &cell : cell_iterator_range)
1246 {
1247 cell_index_offset[cell_index] = n_faces;
1248 n_faces += cell->n_faces();
1249 ++cell_index;
1250 }
1251
1252 n_q_points_unvectorized.reserve(n_faces);
1253
1254 cell_type.reserve(n_faces);
1255
1256 // fill unit points index offset vector
1257 unit_points_index.resize(n_faces + 1);
1258 data_index_offsets.resize(n_faces + 1);
1259 cell_index = 0;
1260 unsigned int n_unit_points = 0;
1261 unsigned int n_data_points = 0;
1262 for (const auto &cell : cell_iterator_range)
1263 {
1264 for (const auto &f : cell->face_indices())
1265 {
1266 const unsigned int current_face_index =
1267 cell_index_offset[cell_index] + f;
1268
1269 unit_points_index[current_face_index] = n_unit_points;
1270 data_index_offsets[current_face_index] = n_data_points;
1271
1272 const unsigned int n_points =
1273 quadrature_vector[cell_index][f].size();
1274 n_q_points_unvectorized.push_back(n_points);
1275
1276 const unsigned int n_q_points =
1277 compute_n_q_points<VectorizedArrayType>(n_points);
1278 n_unit_points += n_q_points;
1279
1280 const unsigned int n_q_points_data =
1281 compute_n_q_points<Number>(n_points);
1282 n_data_points += n_q_points_data;
1283 }
1284
1285 ++cell_index;
1286 }
1287 unit_points_index[n_faces] = n_unit_points;
1288 data_index_offsets[n_faces] = n_data_points;
1289
1290 // compress indices
1291 if (do_cell_index_compression)
1292 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1294
1295 // fill unit points and mapping data for every face of all cells
1296 // resize data vectors
1297 resize_unit_points_faces(n_unit_points);
1298 resize_data_fields(n_data_points);
1299
1300 MappingData mapping_data_previous_cell;
1301 MappingData mapping_data_first;
1302 bool first_set = false;
1303 unsigned int size_compressed_data = 0;
1304 cell_index = 0;
1305 for (const auto &cell : cell_iterator_range)
1306 {
1307 const auto &quadratures_on_faces = quadrature_vector[cell_index];
1308
1309 Assert(quadratures_on_faces.size() == cell->n_faces(),
1310 ExcDimensionMismatch(quadratures_on_faces.size(),
1311 cell->n_faces()));
1312
1313 for (const auto &f : cell->face_indices())
1314 {
1315 const auto &quadrature_on_face = quadratures_on_faces[f];
1316 const bool empty = quadrature_on_face.empty();
1317
1318 const unsigned int current_face_index =
1319 cell_index_offset[cell_index] + f;
1320
1321 // store unit points
1322 const unsigned int n_q_points =
1323 compute_n_q_points<VectorizedArrayType>(
1324 n_q_points_unvectorized[current_face_index]);
1325 store_unit_points_faces(unit_points_index[current_face_index],
1326 n_q_points,
1327 n_q_points_unvectorized[current_face_index],
1328 quadrature_on_face.get_points());
1329
1332 update_flags_mapping,
1333 cell,
1334 f,
1335 quadrature_on_face,
1336 internal_mapping_data,
1337 mapping_data);
1338
1339 // check for cartesian/affine cell
1340 if (!empty &&
1341 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1342 {
1343 cell_type.push_back(internal::compute_geometry_type(
1344 cell->diameter(), mapping_data.inverse_jacobians));
1345
1346 if (!first_set)
1347 {
1348 mapping_data_first = mapping_data;
1349 first_set = true;
1350 }
1351 }
1352 else
1353 cell_type.push_back(
1355
1356 if (current_face_index > 0)
1357 {
1358 // check if current and previous cell are affine
1359 const bool affine_cells =
1360 cell_type[current_face_index] <=
1362 cell_type[current_face_index - 1] <=
1364
1365 // create a comparator to compare inverse Jacobian of current
1366 // and previous cell
1368 1e4 / cell->diameter() *
1369 std::numeric_limits<double>::epsilon() * 1024.);
1370
1371 // we can only compare if current and previous cell have at
1372 // least one quadrature point and both cells are at least affine
1373 const auto comparison_result =
1374 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1375 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1377 comparator.compare(
1378 mapping_data.inverse_jacobians[0],
1379 mapping_data_previous_cell.inverse_jacobians[0]);
1380
1381 // we can compress the Jacobians and inverse Jacobians if
1382 // inverse Jacobians are equal and cells are affine
1383 if (affine_cells &&
1384 comparison_result ==
1386 {
1387 compressed_data_index_offsets.push_back(
1388 compressed_data_index_offsets.back());
1389 }
1390 else if (first_set &&
1391 (cell_type[current_face_index] <=
1393 (comparator.compare(
1394 mapping_data.inverse_jacobians[0],
1395 mapping_data_first.inverse_jacobians[0]) ==
1397 double>::ComparisonResult::equal))
1398 {
1399 compressed_data_index_offsets.push_back(0);
1400 }
1401 else
1402 {
1403 const unsigned int n_compressed_data_last_cell =
1404 cell_type[current_face_index - 1] <=
1406 1 :
1407 compute_n_q_points<Number>(
1408 n_q_points_unvectorized[current_face_index - 1]);
1409
1410 compressed_data_index_offsets.push_back(
1411 compressed_data_index_offsets.back() +
1412 n_compressed_data_last_cell);
1413 }
1414 }
1415 else
1416 compressed_data_index_offsets.push_back(0);
1417
1418 // cache mapping_data from previous cell
1419 mapping_data_previous_cell = mapping_data;
1420
1421 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1422 n_q_points_unvectorized[current_face_index]);
1423 store_mapping_data(data_index_offsets[current_face_index],
1424 n_q_points_data,
1425 n_q_points_unvectorized[current_face_index],
1426 mapping_data,
1427 quadrature_on_face.get_weights(),
1428 data_index_offsets[current_face_index],
1429 cell_type[current_face_index] <=
1431
1432 // update size of compressed data depending on cell type and handle
1433 // empty quadratures
1434 if (cell_type[current_face_index] <=
1436 size_compressed_data = compressed_data_index_offsets.back() + 1;
1437 else
1438 size_compressed_data =
1439 std::max(size_compressed_data,
1440 compressed_data_index_offsets.back() +
1441 n_q_points_data);
1442 }
1443 if (do_cell_index_compression)
1444 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1445 cell_index;
1446
1447 ++cell_index;
1448 }
1449
1450 if (update_flags_mapping & UpdateFlags::update_jacobians)
1451 {
1452 jacobians[0].resize(size_compressed_data);
1453 jacobians[0].shrink_to_fit();
1454 }
1455 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1456 {
1457 inverse_jacobians[0].resize(size_compressed_data);
1458 inverse_jacobians[0].shrink_to_fit();
1459 }
1460
1461 state = State::faces_on_cells_in_vector;
1462 }
1463
1464
1465
1466 template <int dim, int spacedim, typename Number>
1467 template <typename CellIteratorType>
1468 void
1470 const std::vector<std::pair<CellIteratorType, unsigned int>>
1471 &face_iterator_range_interior,
1472 const std::vector<Quadrature<dim - 1>> &quadrature_vector)
1473 {
1474 clear();
1475
1476 do_cell_index_compression = false;
1477
1478 Assert(additional_data.store_cells == false, ExcNotImplemented());
1479
1480
1481 const unsigned int n_faces = quadrature_vector.size();
1482 AssertDimension(n_faces,
1483 std::distance(face_iterator_range_interior.begin(),
1484 face_iterator_range_interior.end()));
1485
1486 n_q_points_unvectorized.reserve(n_faces);
1487
1488 cell_type.reserve(n_faces);
1489 face_number.reserve(n_faces);
1490
1491 // fill unit points index offset vector
1492 unit_points_index.reserve(n_faces + 1);
1493 unit_points_index.push_back(0);
1494 data_index_offsets.reserve(n_faces + 1);
1495 data_index_offsets.push_back(0);
1496 for (const auto &quadrature : quadrature_vector)
1497 {
1498 const unsigned int n_points = quadrature.size();
1499 n_q_points_unvectorized.push_back(n_points);
1500
1501 const unsigned int n_q_points =
1502 compute_n_q_points<VectorizedArrayType>(n_points);
1503 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1504
1505 const unsigned int n_q_points_data =
1506 compute_n_q_points<Number>(n_points);
1507 data_index_offsets.push_back(data_index_offsets.back() +
1508 n_q_points_data);
1509 }
1510
1511 const unsigned int n_unit_points = unit_points_index.back();
1512 const unsigned int n_data_points = data_index_offsets.back();
1513
1514 // resize data vectors
1515 resize_unit_points_faces(n_unit_points);
1516 resize_data_fields(n_data_points, true);
1517
1518 std::array<MappingData, 2> mapping_data;
1519 std::array<MappingData, 2> mapping_data_previous_cell;
1520 std::array<MappingData, 2> mapping_data_first;
1521 bool first_set = false;
1522 unsigned int size_compressed_data = 0;
1523 unsigned int face_index = 0;
1524 for (const auto &cell_and_f : face_iterator_range_interior)
1525 {
1526 const auto &quadrature_on_face = quadrature_vector[face_index];
1527 const bool empty = quadrature_on_face.empty();
1528
1529 // get interior cell and face number
1530 const auto &cell_m = cell_and_f.first;
1531 const auto f_m = cell_and_f.second;
1532
1533 // get exterior cell and face number
1534 const auto &cell_p =
1535 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1536 const auto f_p =
1537 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1538
1539 face_number.emplace_back(f_m, f_p);
1540
1541 Assert(
1542 cell_m->combined_face_orientation(f_m) ==
1544 cell_p->combined_face_orientation(f_p) ==
1546 ExcMessage(
1547 "Non standard face orientation is currently not implemented."));
1548
1549 // store unit points
1550 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1551 n_q_points_unvectorized[face_index]);
1552 store_unit_points_faces(unit_points_index[face_index],
1553 n_q_points,
1554 n_q_points_unvectorized[face_index],
1555 quadrature_on_face.get_points());
1556
1557 // compute mapping for interior face
1560 update_flags_mapping,
1561 cell_m,
1562 f_m,
1563 quadrature_on_face,
1564 internal_mapping_data,
1565 mapping_data[0]);
1566
1567 // compute mapping for exterior face
1570 update_flags_mapping,
1571 cell_p,
1572 f_p,
1573 quadrature_on_face,
1574 internal_mapping_data,
1575 mapping_data[1]);
1576
1577 // check for cartesian/affine cell
1578 if (!empty &&
1579 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1580 {
1581 // select more general type of interior and exterior cell
1582 cell_type.push_back(std::max(
1584 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1586 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1587
1588 // cache mapping data of first cell pair with non-empty quadrature
1589 // on the face
1590 if (!first_set)
1591 {
1592 mapping_data_first = mapping_data;
1593 first_set = true;
1594 }
1595 }
1596 else
1597 cell_type.push_back(
1599
1600 if (face_index > 0)
1601 {
1602 // check if current and previous cell pairs are affine
1603 const bool affine_cells =
1604 cell_type[face_index] <=
1606 cell_type[face_index - 1] <=
1608
1609 // create a comparator to compare inverse Jacobian of current
1610 // and previous cell pair
1612 1e4 / cell_m->diameter() *
1613 std::numeric_limits<double>::epsilon() * 1024.);
1614
1615 // we can only compare if current and previous cell have at
1616 // least one quadrature point and both cells are at least affine
1617 const auto comparison_result_m =
1618 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1619 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1621 comparator.compare(
1622 mapping_data[0].inverse_jacobians[0],
1623 mapping_data_previous_cell[0].inverse_jacobians[0]);
1624
1625 const auto comparison_result_p =
1626 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1627 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1629 comparator.compare(
1630 mapping_data[1].inverse_jacobians[0],
1631 mapping_data_previous_cell[1].inverse_jacobians[0]);
1632
1633 // we can compress the Jacobians and inverse Jacobians if
1634 // inverse Jacobians are equal and cells are affine
1635 if (affine_cells &&
1636 comparison_result_m ==
1638 comparison_result_p ==
1640 {
1641 compressed_data_index_offsets.push_back(
1642 compressed_data_index_offsets.back());
1643 }
1644 else if (first_set &&
1645 (cell_type[face_index] <=
1647 (comparator.compare(
1648 mapping_data[0].inverse_jacobians[0],
1649 mapping_data_first[0].inverse_jacobians[0]) ==
1651 double>::ComparisonResult::equal) &&
1652 (comparator.compare(
1653 mapping_data[1].inverse_jacobians[0],
1654 mapping_data_first[1].inverse_jacobians[0]) ==
1656 {
1657 compressed_data_index_offsets.push_back(0);
1658 }
1659 else
1660 {
1661 const unsigned int n_compressed_data_last_cell =
1662 cell_type[face_index - 1] <=
1664 1 :
1665 compute_n_q_points<Number>(
1666 n_q_points_unvectorized[face_index - 1]);
1667
1668 compressed_data_index_offsets.push_back(
1669 compressed_data_index_offsets.back() +
1670 n_compressed_data_last_cell);
1671 }
1672 }
1673 else
1674 compressed_data_index_offsets.push_back(0);
1675
1676 // cache mapping_data from previous cell pair
1677 mapping_data_previous_cell = mapping_data;
1678
1679 const unsigned int n_q_points_data =
1680 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1681
1682 // store mapping data of interior face
1683 store_mapping_data(data_index_offsets[face_index],
1684 n_q_points_data,
1685 n_q_points_unvectorized[face_index],
1686 mapping_data[0],
1687 quadrature_on_face.get_weights(),
1688 data_index_offsets[face_index],
1689 cell_type[face_index] <=
1691 true);
1692
1693 // store only necessary mapping data for exterior face (Jacobians and
1694 // inverse Jacobians)
1695 store_mapping_data(data_index_offsets[face_index],
1696 n_q_points_data,
1697 n_q_points_unvectorized[face_index],
1698 mapping_data[1],
1699 quadrature_on_face.get_weights(),
1700 data_index_offsets[face_index],
1701 cell_type[face_index] <=
1703 false);
1704
1705 // update size of compressed data depending on cell type and handle
1706 // empty quadratures
1707 if (cell_type[face_index] <=
1709 size_compressed_data = compressed_data_index_offsets.back() + 1;
1710 else
1711 size_compressed_data =
1712 std::max(size_compressed_data,
1713 compressed_data_index_offsets.back() + n_q_points_data);
1714
1715 ++face_index;
1716 }
1717
1718 if (update_flags_mapping & UpdateFlags::update_jacobians)
1719 {
1720 jacobians[0].resize(size_compressed_data);
1721 jacobians[0].shrink_to_fit();
1722 jacobians[1].resize(size_compressed_data);
1723 jacobians[1].shrink_to_fit();
1724 }
1725 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1726 {
1727 inverse_jacobians[0].resize(size_compressed_data);
1728 inverse_jacobians[0].shrink_to_fit();
1729 inverse_jacobians[1].resize(size_compressed_data);
1730 inverse_jacobians[1].shrink_to_fit();
1731 }
1732
1733 state = State::face_vector;
1734 }
1735
1736
1737
1738 template <int dim, int spacedim, typename Number>
1739 bool
1741 {
1742 return state == State::faces_on_cells_in_vector;
1743 }
1744
1745
1746
1747 template <int dim, int spacedim, typename Number>
1748 unsigned int
1750 const unsigned int geometry_index) const
1751 {
1752 return n_q_points_unvectorized[geometry_index];
1753 }
1754
1755
1756 template <int dim, int spacedim, typename Number>
1759 const unsigned int geometry_index) const
1760 {
1761 return cell_type[geometry_index];
1762 }
1763
1764
1765
1766 template <int dim, int spacedim, typename Number>
1769 const unsigned int cell_index) const
1770 {
1771 Assert(
1772 additional_data.store_cells,
1773 ExcMessage(
1774 "Cells have been not stored. You can enable this by Additional::store_cells."));
1775 return {triangulation.get(),
1776 cell_level_and_indices[cell_index].first,
1777 cell_level_and_indices[cell_index].second};
1778 }
1779
1780
1781
1782 template <int dim, int spacedim, typename Number>
1783 template <typename NumberType>
1784 unsigned int
1786 const unsigned int n_q_points_unvectorized)
1787 {
1788 const unsigned int n_lanes =
1790 const unsigned int n_filled_lanes_last_batch =
1791 n_q_points_unvectorized % n_lanes;
1792 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1793 if (n_filled_lanes_last_batch > 0)
1794 ++n_q_points;
1795 return n_q_points;
1796 }
1797
1798
1799
1800 template <int dim, int spacedim, typename Number>
1801 template <bool is_face>
1802 unsigned int
1804 const unsigned int cell_index,
1805 const unsigned int face_number) const
1806 {
1808 return 0;
1809
1810 const unsigned int compressed_cell_index =
1811 compute_compressed_cell_index(cell_index);
1812 if (!is_face)
1813 {
1814 Assert(state == State::cell_vector,
1815 ExcMessage(
1816 "This mapping info is not reinitialized for a cell vector!"));
1817 return compressed_cell_index;
1818 }
1819 else
1820 {
1822 ExcMessage(
1823 "cell_index has to be set if face_number is specified!"));
1824 Assert(state == State::faces_on_cells_in_vector ||
1825 state == State::face_vector,
1826 ExcMessage("This mapping info is not reinitialized for faces"
1827 " on cells in a vector!"));
1828 if (state == State::faces_on_cells_in_vector)
1829 return cell_index_offset[compressed_cell_index] + face_number;
1830 else if (state == State::face_vector)
1831 return cell_index;
1832 }
1833 }
1834
1835
1836
1837 template <int dim, int spacedim, typename Number>
1838 unsigned int
1840 const unsigned int cell_index) const
1841 {
1842 if (do_cell_index_compression)
1843 {
1844 Assert(cell_index_to_compressed_cell_index[cell_index] !=
1846 ExcMessage("Mapping info object was not initialized for this"
1847 " active cell index!"));
1848 return cell_index_to_compressed_cell_index[cell_index];
1849 }
1850 else
1851 return cell_index;
1852 }
1853
1854
1855 template <int dim, int spacedim, typename Number>
1856 void
1858 const unsigned int unit_points_index_offset,
1859 const unsigned int n_q_points,
1860 const unsigned int n_q_points_unvectorized,
1861 const std::vector<Point<dim>> &points)
1862 {
1863 const unsigned int n_lanes =
1865
1866 for (unsigned int q = 0; q < n_q_points; ++q)
1867 {
1868 const unsigned int offset = unit_points_index_offset + q;
1869 for (unsigned int v = 0;
1870 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1871 ++v)
1872 for (unsigned int d = 0; d < dim; ++d)
1874 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1875 }
1876 }
1877
1878
1879
1880 template <int dim, int spacedim, typename Number>
1881 void
1883 const unsigned int unit_points_index_offset,
1884 const unsigned int n_q_points,
1885 const unsigned int n_q_points_unvectorized,
1886 const std::vector<Point<dim - 1>> &points)
1887 {
1888 const unsigned int n_lanes =
1890
1891 for (unsigned int q = 0; q < n_q_points; ++q)
1892 {
1893 const unsigned int offset = unit_points_index_offset + q;
1894 for (unsigned int v = 0;
1895 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1896 ++v)
1897 for (unsigned int d = 0; d < dim - 1; ++d)
1899 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1900 }
1901 }
1902
1903
1904
1905 template <int dim, int spacedim, typename Number>
1906 void
1908 const unsigned int unit_points_index_offset,
1909 const unsigned int n_q_points,
1910 const unsigned int n_q_points_unvectorized,
1911 const MappingInfo::MappingData &mapping_data,
1912 const std::vector<double> &weights,
1913 const unsigned int compressed_unit_point_index_offset,
1914 const bool affine_cell,
1915 const bool is_interior)
1916 {
1917 const unsigned int n_lanes =
1919
1920 for (unsigned int q = 0; q < n_q_points; ++q)
1921 {
1922 const unsigned int offset = unit_points_index_offset + q;
1923 const unsigned int compressed_offset =
1924 compressed_unit_point_index_offset + q;
1925 for (unsigned int v = 0;
1926 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1927 ++v)
1928 {
1929 if (q == 0 || !affine_cell)
1930 {
1931 if (update_flags_mapping & UpdateFlags::update_jacobians)
1932 for (unsigned int d = 0; d < dim; ++d)
1933 for (unsigned int s = 0; s < spacedim; ++s)
1935 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1936 v) = mapping_data.jacobians[q * n_lanes + v][d][s];
1937 if (update_flags_mapping &
1939 for (unsigned int d = 0; d < dim; ++d)
1940 for (unsigned int s = 0; s < spacedim; ++s)
1942 inverse_jacobians[is_interior ? 0 : 1]
1943 [compressed_offset][d][s],
1944 v) =
1945 mapping_data.inverse_jacobians[q * n_lanes + v][d][s];
1946 }
1947
1948 if (is_interior)
1949 {
1950 if (update_flags_mapping & UpdateFlags::update_JxW_values)
1951 {
1952 if (additional_data.use_global_weights)
1953 {
1955 JxW_values[offset], v) = weights[q * n_lanes + v];
1956 }
1957 else
1958 {
1960 JxW_values[offset], v) =
1961 mapping_data.JxW_values[q * n_lanes + v];
1962 }
1963 }
1964 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
1965 for (unsigned int s = 0; s < spacedim; ++s)
1967 normal_vectors[offset][s], v) =
1968 mapping_data.normal_vectors[q * n_lanes + v][s];
1969 if (update_flags_mapping &
1971 for (unsigned int s = 0; s < spacedim; ++s)
1973 real_points[offset][s], v) =
1974 mapping_data.quadrature_points[q * n_lanes + v][s];
1975 }
1976 }
1977 }
1978 }
1979
1980
1981
1982 template <int dim, int spacedim, typename Number>
1983 void
1985 const unsigned int n_unit_point_batches)
1986 {
1987 unit_points.resize(n_unit_point_batches);
1988 }
1989
1990
1991
1992 template <int dim, int spacedim, typename Number>
1993 void
1995 const unsigned int n_unit_point_batches)
1996 {
1997 unit_points_faces.resize(n_unit_point_batches);
1998 }
1999
2000
2001
2002 template <int dim, int spacedim, typename Number>
2003 void
2005 const unsigned int n_data_point_batches,
2006 const bool is_face_centric)
2007 {
2008 if (update_flags_mapping & UpdateFlags::update_jacobians)
2009 {
2010 jacobians[0].resize(n_data_point_batches);
2011 if (is_face_centric)
2012 jacobians[1].resize(n_data_point_batches);
2013 }
2014 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2015 {
2016 inverse_jacobians[0].resize(n_data_point_batches);
2017 if (is_face_centric)
2018 inverse_jacobians[1].resize(n_data_point_batches);
2019 }
2020 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2021 JxW_values.resize(n_data_point_batches);
2022 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2023 normal_vectors.resize(n_data_point_batches);
2024 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2025 real_points.resize(n_data_point_batches);
2026 }
2027
2028
2029
2030 template <int dim, int spacedim, typename Number>
2031 inline const Point<
2032 dim,
2035 const unsigned int offset) const
2036 {
2037 return unit_points.data() + offset;
2038 }
2039
2040
2041
2042 template <int dim, int spacedim, typename Number>
2043 inline const Point<
2044 dim - 1,
2047 const unsigned int offset) const
2048 {
2049 return unit_points_faces.data() + offset;
2050 }
2051
2052
2053
2054 template <int dim, int spacedim, typename Number>
2055 inline const Point<spacedim, Number> *
2057 const unsigned int offset) const
2058 {
2059 return real_points.data() + offset;
2060 }
2061
2062
2063
2064 template <int dim, int spacedim, typename Number>
2065 unsigned int
2067 const unsigned int geometry_index) const
2068 {
2069 return unit_points_index[geometry_index];
2070 }
2071
2072
2073
2074 template <int dim, int spacedim, typename Number>
2075 unsigned int
2077 const unsigned int geometry_index) const
2078 {
2079 return data_index_offsets[geometry_index];
2080 }
2081
2082
2083 template <int dim, int spacedim, typename Number>
2084 unsigned int
2086 const unsigned int geometry_index) const
2087 {
2088 return compressed_data_index_offsets[geometry_index];
2089 }
2090
2091
2092
2093 template <int dim, int spacedim, typename Number>
2096 const bool is_interior) const
2097 {
2098 return jacobians[is_interior ? 0 : 1].data() + offset;
2099 }
2100
2101
2102
2103 template <int dim, int spacedim, typename Number>
2106 const unsigned int offset,
2107 const bool is_interior) const
2108 {
2109 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2110 }
2111
2112
2113
2114 template <int dim, int spacedim, typename Number>
2115 inline const Tensor<1, spacedim, Number> *
2117 const unsigned int offset) const
2118 {
2119 return normal_vectors.data() + offset;
2120 }
2121
2122
2123
2124 template <int dim, int spacedim, typename Number>
2125 inline const Number *
2126 MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
2127 {
2128 return JxW_values.data() + offset;
2129 }
2130
2131
2132
2133 template <int dim, int spacedim, typename Number>
2136 {
2137 return *mapping;
2138 }
2139
2140
2141
2142 template <int dim, int spacedim, typename Number>
2145 {
2146 return update_flags;
2147 }
2148
2149
2150
2151 template <int dim, int spacedim, typename Number>
2154 {
2155 return update_flags_mapping;
2156 }
2157
2158
2159
2160 template <int dim, int spacedim, typename Number>
2161 std::size_t
2163 {
2164 std::size_t memory = MemoryConsumption::memory_consumption(unit_points);
2165 memory += MemoryConsumption::memory_consumption(unit_points_faces);
2166 memory += MemoryConsumption::memory_consumption(unit_points_index);
2167 memory += cell_type.capacity() *
2169 memory += MemoryConsumption::memory_consumption(data_index_offsets);
2170 memory +=
2171 MemoryConsumption::memory_consumption(compressed_data_index_offsets);
2172 memory += MemoryConsumption::memory_consumption(JxW_values);
2173 memory += MemoryConsumption::memory_consumption(normal_vectors);
2174 memory += MemoryConsumption::memory_consumption(jacobians);
2175 memory += MemoryConsumption::memory_consumption(inverse_jacobians);
2176 memory += MemoryConsumption::memory_consumption(real_points);
2177 memory += MemoryConsumption::memory_consumption(n_q_points_unvectorized);
2178 memory += MemoryConsumption::memory_consumption(cell_index_offset);
2180 cell_index_to_compressed_cell_index);
2181 memory += MemoryConsumption::memory_consumption(cell_level_and_indices);
2182 memory += sizeof(*this);
2183 return memory;
2184 }
2185} // namespace NonMatching
2186
2188
2189#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: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
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
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)