Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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::shared_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
74 // reuse internal_mapping_data for MappingQ to avoid memory allocations
75 if (const MappingQ<dim, spacedim> *mapping_q =
76 dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
77 {
78 (void)mapping_q;
79 auto &data =
80 dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
81 *internal_mapping_data);
82 data.initialize(update_flags_mapping,
83 quadrature,
84 quadrature.size());
85 }
86 else
87 {
88 internal_mapping_data =
89 mapping->get_data(update_flags_mapping, quadrature);
90 }
91
92 cell_similarity = mapping->fill_fe_values(cell,
93 cell_similarity,
94 quadrature,
95 *internal_mapping_data,
96 mapping_data);
97 }
98
99
100
101 static void
103 const SmartPointer<const Mapping<dim, spacedim>> mapping,
104 const UpdateFlags &update_flags_mapping,
106 const ImmersedSurfaceQuadrature<dim> &quadrature,
107 std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
108 internal_mapping_data,
109 MappingData &mapping_data)
110 {
111 mapping_data.initialize(quadrature.size(), update_flags_mapping);
112
113 // reuse internal_mapping_data for MappingQ to avoid memory allocations
114 if (dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
115 {
116 auto &data =
117 dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
118 *internal_mapping_data);
119 data.initialize(update_flags_mapping,
120 quadrature,
121 quadrature.size());
122 }
123 else
124 {
125 internal_mapping_data =
126 mapping->get_data(update_flags_mapping, quadrature);
127 }
128
129 mapping->fill_fe_immersed_surface_values(cell,
130 quadrature,
131 *internal_mapping_data,
132 mapping_data);
133 }
134
135
136
137 static void
139 const SmartPointer<const Mapping<dim, spacedim>> mapping,
140 const UpdateFlags &update_flags_mapping,
142 const unsigned int face_no,
143 const Quadrature<dim - 1> &quadrature,
144 std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
145 internal_mapping_data,
146 MappingData &mapping_data)
147 {
148 mapping_data.initialize(quadrature.size(), update_flags_mapping);
149
150 // reuse internal_mapping_data for MappingQ to avoid memory allocations
151 if (const MappingQ<dim, spacedim> *mapping_q =
152 dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
153 {
154 auto &data =
155 dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
156 *internal_mapping_data);
157 data.initialize_face(update_flags_mapping,
159 ReferenceCells::get_hypercube<dim>(),
160 quadrature,
161 face_no,
162 cell->face_orientation(face_no),
163 cell->face_flip(face_no),
164 cell->face_rotation(face_no)),
165 quadrature.size());
166
167 mapping_q->fill_mapping_data_for_face_quadrature(
168 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
169 }
170 else
171 {
172 auto internal_mapping_data =
173 mapping->get_face_data(update_flags_mapping,
174 hp::QCollection<dim - 1>(quadrature));
175
176 mapping->fill_fe_face_values(cell,
177 face_no,
178 hp::QCollection<dim - 1>(quadrature),
179 *internal_mapping_data,
180 mapping_data);
181 }
182 }
183 };
184
185 template <int dim, int spacedim = dim>
188 const double diameter,
190 &inverse_jacobians)
191 {
192 const auto jac_0 = inverse_jacobians[0];
193 const double zero_tolerance_double =
194 1e4 / diameter * std::numeric_limits<double>::epsilon() * 1024.;
195 bool jacobian_constant = true;
196 for (unsigned int q = 1; q < inverse_jacobians.size(); ++q)
197 {
198 const DerivativeForm<1, spacedim, dim> &jac = inverse_jacobians[q];
199 for (unsigned int d = 0; d < dim; ++d)
200 for (unsigned int e = 0; e < spacedim; ++e)
201 if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
202 jacobian_constant = false;
203 if (!jacobian_constant)
204 break;
205 }
206
207 // check whether the Jacobian is diagonal to machine
208 // accuracy
209 bool cell_cartesian = jacobian_constant && dim == spacedim;
210 for (unsigned int d = 0; d < dim; ++d)
211 for (unsigned int e = 0; e < dim; ++e)
212 if (d != e)
213 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
214 {
215 cell_cartesian = false;
216 break;
217 }
218
219 // return cell type
220 if (cell_cartesian)
221 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
222 else if (jacobian_constant)
223 return ::internal::MatrixFreeFunctions::GeometryType::affine;
224 else
225 return ::internal::MatrixFreeFunctions::GeometryType::general;
226 }
227 } // namespace internal
228
245 template <int dim, int spacedim = dim, typename Number = double>
247 {
248 public:
252 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
253 Number>::vectorized_value_type;
254
260 {
265 const bool store_cells = false)
268 {}
269
276
285 };
286
301 const AdditionalData additional_data = AdditionalData());
302
307 void
309 const std::vector<Point<dim>> &unit_points);
310
315 void
317 const ArrayView<const Point<dim>> &unit_points);
318
324 void
326 const Quadrature<dim> &quadrature);
327
338 template <typename ContainerType>
339 void
341 const ContainerType &cell_iterator_range,
342 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
343 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
344
351 template <typename ContainerType>
352 void
354 const ContainerType &cell_iterator_range,
355 const std::vector<Quadrature<dim>> &quadrature_vector,
356 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
357
362 template <typename ContainerType>
363 void
365 const ContainerType &cell_iterator_range,
366 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
367 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
368
373 template <typename ContainerType>
374 void
376 const ContainerType &cell_iterator_range,
377 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
378 const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
379
384 template <typename CellIteratorType>
385 void
386 reinit_faces(const std::vector<std::pair<CellIteratorType, unsigned int>>
387 &face_iterator_range_interior,
388 const std::vector<Quadrature<dim - 1>> &quadrature_vector);
389
394 bool
396
400 unsigned int
401 get_face_number(const unsigned int offset, const bool is_interior) const;
402
408 get_unit_point(const unsigned int offset) const;
409
414 const Point<dim - 1, VectorizedArrayType> *
415 get_unit_point_faces(const unsigned int offset) const;
416
422 get_jacobian(const unsigned int offset,
423 const bool is_interior = true) const;
424
430 get_inverse_jacobian(const unsigned int offset,
431 const bool is_interior = true) const;
432
438 get_normal_vector(const unsigned int offset) const;
439
444 const Number *
445 get_JxW(const unsigned int offset) const;
446
452 get_real_point(const unsigned int offset) const;
453
458 get_mapping() const;
459
465
471
475 boost::signals2::connection
476 connect_is_reinitialized(const std::function<void()> &set_is_reinitialized);
477
481 template <bool is_face>
482 unsigned int
484 const unsigned int face_number) const;
485
489 unsigned int
490 compute_unit_point_index_offset(const unsigned int geometry_index) const;
491
495 unsigned int
496 compute_data_index_offset(const unsigned int geometry_index) const;
497
501 unsigned int
503 const unsigned int geometry_index) const;
504
508 unsigned int
509 get_n_q_points_unvectorized(const unsigned int geometry_index) const;
510
515 get_cell_type(const unsigned int geometry_index) const;
516
523 get_cell_iterator(const unsigned int cell_index) const;
524
528 std::size_t
530
531 private:
534 spacedim>;
535
539 template <typename NumberType>
540 unsigned int
542
546 void
548
552 void
553 resize_unit_points(const unsigned int n_unit_point_batches);
554
558 void
559 resize_unit_points_faces(const unsigned int n_unit_point_batches);
560
564 void
565 resize_data_fields(const unsigned int n_data_point_batches,
566 const bool is_face_centric = false);
567
571 void
572 store_unit_points(const unsigned int unit_points_index_offset,
573 const unsigned int n_q_points,
574 const unsigned int n_q_points_unvectorized,
575 const std::vector<Point<dim>> &points);
576
580 void
581 store_unit_points_faces(const unsigned int unit_points_index_offset,
582 const unsigned int n_q_points,
583 const unsigned int n_q_points_unvectorized,
584 const std::vector<Point<dim - 1>> &points);
585
589 void
590 store_mapping_data(const unsigned int unit_points_index_offset,
591 const unsigned int n_q_points,
592 const unsigned int n_q_points_unvectorized,
593 const MappingData &mapping_data,
594 const std::vector<double> &weights,
595 const unsigned int compressed_unit_point_index_offset,
596 const bool affine_cell,
597 const bool is_interior = true);
598
602 unsigned int
603 compute_compressed_cell_index(const unsigned int cell_index) const;
604
608 template <typename ContainerType, typename QuadratureType>
609 void
611 const ContainerType &cell_iterator_range,
612 const std::vector<QuadratureType> &quadrature_vector,
613 const unsigned int n_unfiltered_cells,
614 const std::function<
615 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
616 const QuadratureType &quadrature,
617 MappingData &mapping_data)> &compute_mapping_data);
618
622 enum class State
623 {
624 invalid,
629 };
630
636
643
650
654 std::vector<unsigned int> unit_points_index;
655
659 std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
661
666
671
676
681
688 std::vector<::internal::MatrixFreeFunctions::GeometryType> cell_type;
689
693 std::vector<unsigned int> data_index_offsets;
694
698 std::vector<unsigned int> compressed_data_index_offsets;
699
707
714
721 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
723
731 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
733
740
744 std::vector<unsigned int> n_q_points_unvectorized;
745
750 std::vector<unsigned int> cell_index_offset;
751
756 std::vector<unsigned int> cell_index_to_compressed_cell_index;
757
762
767 boost::signals2::signal<void()> is_reinitialized;
768
775
780 std::vector<std::pair<int, int>> cell_level_and_indices;
781
782 std::vector<std::pair<unsigned char, unsigned char>> face_number;
783 };
784
785 template <int dim, int spacedim, typename Number>
786 inline unsigned int
788 const unsigned int offset,
789 const bool is_interior) const
790 {
791 const auto &face_pair = face_number[offset];
792 return is_interior ? face_pair.first : face_pair.second;
793 }
794
795 // ----------------------- template functions ----------------------
796
797
798 template <int dim, int spacedim, typename Number>
800 const Mapping<dim, spacedim> &mapping,
801 const UpdateFlags update_flags,
802 const AdditionalData additional_data)
803 : mapping(&mapping)
804 , update_flags(update_flags)
805 , update_flags_mapping(update_default)
806 , additional_data(additional_data)
807 {
808 // translate update flags
818
819 // always save quadrature points for now
821
824 this->mapping, update_flags_mapping);
825
826 // construct internal_mapping_data for MappingQ to be able to reuse it in
827 // reinit() calls to avoid memory allocations
828 if (const MappingQ<dim, spacedim> *mapping_q =
829 dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
830 {
832 std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
833 mapping_q->get_degree());
834 }
835 }
836
837
838
839 template <int dim, int spacedim, typename Number>
840 void
842 {
843 n_q_points_unvectorized.clear();
844 unit_points_index.clear();
845 data_index_offsets.clear();
846 compressed_data_index_offsets.clear();
847 cell_type.clear();
848 }
849
850
851
852 template <int dim, int spacedim, typename Number>
853 void
856 const std::vector<Point<dim>> &unit_points_in)
857 {
858 reinit(cell, Quadrature<dim>(unit_points_in));
859 }
860
861
862
863 template <int dim, int spacedim, typename Number>
864 void
867 const ArrayView<const Point<dim>> &unit_points_in)
868 {
869 reinit(cell,
870 std::vector<Point<dim>>(unit_points_in.begin(),
871 unit_points_in.end()));
872 }
873
874
875
876 template <int dim, int spacedim, typename Number>
877 void
880 const Quadrature<dim> &quadrature)
881 {
882 n_q_points_unvectorized.resize(1);
883 n_q_points_unvectorized[0] = quadrature.size();
884
885 const unsigned int n_q_points =
886 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
887
888 const unsigned int n_q_points_data =
889 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
890
891 // resize data vectors
892 resize_unit_points(n_q_points);
893 resize_data_fields(n_q_points_data);
894
895 // store unit points
896 store_unit_points(0,
897 n_q_points,
898 n_q_points_unvectorized[0],
899 quadrature.get_points());
900
901 // compute mapping data
902 MappingData mapping_data;
906 update_flags_mapping,
907 cell,
908 cell_similarity,
909 quadrature,
910 internal_mapping_data,
911 mapping_data);
912
913 // check for cartesian/affine cell
914 if (!quadrature.empty() &&
915 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
916 {
917 cell_type.push_back(
918 internal::compute_geometry_type(cell->diameter(),
919 mapping_data.inverse_jacobians));
920 }
921 else
922 cell_type.push_back(
924
925 // store mapping data
926 store_mapping_data(
927 0,
928 n_q_points_data,
929 n_q_points_unvectorized[0],
930 mapping_data,
931 quadrature.get_weights(),
932 0,
933 cell_type.back() <=
935
936 unit_points_index.push_back(0);
937 data_index_offsets.push_back(0);
938 compressed_data_index_offsets.push_back(0);
939
940 state = State::single_cell;
941 is_reinitialized();
942 }
943
944
945
946 template <int dim, int spacedim, typename Number>
947 template <typename ContainerType>
948 void
950 const ContainerType &cell_iterator_range,
951 const std::vector<std::vector<Point<dim>>> &unit_points_vector,
952 const unsigned int n_unfiltered_cells)
953 {
954 const unsigned int n_cells = unit_points_vector.size();
955 AssertDimension(n_cells,
956 std::distance(cell_iterator_range.begin(),
957 cell_iterator_range.end()));
958
959 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
960 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
961 quadrature_vector[cell_index] =
962 Quadrature<dim>(unit_points_vector[cell_index]);
963
964 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
965 }
966
967
968
969 template <int dim, int spacedim, typename Number>
970 template <typename ContainerType, typename QuadratureType>
971 void
973 const ContainerType &cell_iterator_range,
974 const std::vector<QuadratureType> &quadrature_vector,
975 const unsigned int n_unfiltered_cells,
976 const std::function<
977 void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
978 const QuadratureType &quadrature,
979 MappingData &mapping_data)> &compute_mapping_data)
980 {
981 clear();
982
983 do_cell_index_compression =
984 n_unfiltered_cells != numbers::invalid_unsigned_int;
985
986 const unsigned int n_cells = quadrature_vector.size();
987 AssertDimension(n_cells,
988 std::distance(cell_iterator_range.begin(),
989 cell_iterator_range.end()));
990
991 n_q_points_unvectorized.reserve(n_cells);
992
993 cell_type.reserve(n_cells);
994
995 if (additional_data.store_cells)
996 cell_level_and_indices.resize(n_cells);
997
998 // fill unit points index offset vector
999 unit_points_index.reserve(n_cells + 1);
1000 unit_points_index.push_back(0);
1001 data_index_offsets.reserve(n_cells + 1);
1002 data_index_offsets.push_back(0);
1003 for (const auto &quadrature : quadrature_vector)
1004 {
1005 const unsigned int n_points = quadrature.size();
1006 n_q_points_unvectorized.push_back(n_points);
1007
1008 const unsigned int n_q_points =
1009 compute_n_q_points<VectorizedArrayType>(n_points);
1010 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1011
1012 const unsigned int n_q_points_data =
1013 compute_n_q_points<Number>(n_points);
1014 data_index_offsets.push_back(data_index_offsets.back() +
1015 n_q_points_data);
1016 }
1017
1018 const unsigned int n_unit_points = unit_points_index.back();
1019 const unsigned int n_data_points = data_index_offsets.back();
1020
1021 // resize data vectors
1022 resize_unit_points(n_unit_points);
1023 resize_data_fields(n_data_points);
1024
1025 if (do_cell_index_compression)
1026 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1028
1029 MappingData mapping_data;
1030 MappingData mapping_data_previous_cell;
1031 unsigned int size_compressed_data = 0;
1032 unsigned int cell_index = 0;
1033 for (const auto &cell : cell_iterator_range)
1034 {
1035 if (additional_data.store_cells)
1036 {
1037 this->triangulation = &cell->get_triangulation();
1038 cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
1039 }
1040
1041 const auto &quadrature = quadrature_vector[cell_index];
1042 const bool empty = quadrature.empty();
1043
1044 // store unit points
1045 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1046 n_q_points_unvectorized[cell_index]);
1047 store_unit_points(unit_points_index[cell_index],
1048 n_q_points,
1049 n_q_points_unvectorized[cell_index],
1050 quadrature.get_points());
1051
1052 // compute mapping data
1053 compute_mapping_data(cell, quadrature, mapping_data);
1054
1055 // store mapping data
1056 const unsigned int n_q_points_data =
1057 compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
1058
1059 // check for cartesian/affine cell
1060 if (!empty &&
1061 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1062 {
1063 cell_type.push_back(
1064 internal::compute_geometry_type(cell->diameter(),
1065 mapping_data.inverse_jacobians));
1066 }
1067 else
1068 cell_type.push_back(
1070
1071 if (cell_index > 0)
1072 {
1073 // check if current and previous cell are affine
1074 const bool affine_cells =
1075 cell_type[cell_index] <=
1077 cell_type[cell_index - 1] <=
1079
1080 // create a comparator to compare inverse Jacobian of current
1081 // and previous cell
1083 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1084 1024.);
1085
1086 // we can only compare if current and previous cell have at least
1087 // one quadrature point and both cells are at least affine
1088 const auto comparison_result =
1089 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1090 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1092 comparator.compare(
1093 mapping_data.inverse_jacobians[0],
1094 mapping_data_previous_cell.inverse_jacobians[0]);
1095
1096 // we can compress the Jacobians and inverse Jacobians if
1097 // inverse Jacobians are equal and cells are affine
1098 if (affine_cells &&
1099 comparison_result ==
1101 {
1102 compressed_data_index_offsets.push_back(
1103 compressed_data_index_offsets.back());
1104 }
1105 else
1106 {
1107 const unsigned int n_compressed_data_last_cell =
1108 cell_type[cell_index - 1] <=
1110 1 :
1111 compute_n_q_points<Number>(
1112 n_q_points_unvectorized[cell_index - 1]);
1113
1114 compressed_data_index_offsets.push_back(
1115 compressed_data_index_offsets.back() +
1116 n_compressed_data_last_cell);
1117 }
1118 }
1119 else
1120 compressed_data_index_offsets.push_back(0);
1121
1122 // cache mapping_data from previous cell
1123 mapping_data_previous_cell = mapping_data;
1124
1125 store_mapping_data(data_index_offsets[cell_index],
1126 n_q_points_data,
1127 n_q_points_unvectorized[cell_index],
1128 mapping_data,
1129 quadrature.get_weights(),
1130 compressed_data_index_offsets[cell_index],
1131 cell_type[cell_index] <=
1133
1134 // update size of compressed data depending on cell type and handle
1135 // empty quadratures
1136 if (cell_type[cell_index] <=
1138 size_compressed_data = compressed_data_index_offsets.back() + 1;
1139 else
1140 size_compressed_data =
1141 std::max(size_compressed_data,
1142 compressed_data_index_offsets.back() + n_q_points_data);
1143
1144 if (do_cell_index_compression)
1145 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1146 cell_index;
1147
1148 ++cell_index;
1149 }
1150
1151 if (update_flags_mapping & UpdateFlags::update_jacobians)
1152 {
1153 jacobians[0].resize(size_compressed_data);
1154 jacobians[0].shrink_to_fit();
1155 }
1156 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1157 {
1158 inverse_jacobians[0].resize(size_compressed_data);
1159 inverse_jacobians[0].shrink_to_fit();
1160 }
1161
1162 state = State::cell_vector;
1163 is_reinitialized();
1164 }
1165
1166
1167
1168 template <int dim, int spacedim, typename Number>
1169 template <typename ContainerType>
1170 void
1172 const ContainerType &cell_iterator_range,
1173 const std::vector<Quadrature<dim>> &quadrature_vector,
1174 const unsigned int n_unfiltered_cells)
1175 {
1176 auto compute_mapping_data_for_cells =
1177 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1178 const Quadrature<dim> &quadrature,
1179 MappingData &mapping_data) {
1180 CellSimilarity::Similarity cell_similarity =
1184 update_flags_mapping,
1185 cell,
1186 cell_similarity,
1187 quadrature,
1188 internal_mapping_data,
1189 mapping_data);
1190 };
1191
1192 do_reinit_cells<ContainerType, Quadrature<dim>>(
1193 cell_iterator_range,
1194 quadrature_vector,
1195 n_unfiltered_cells,
1196 compute_mapping_data_for_cells);
1197 }
1198
1199
1200
1201 template <int dim, int spacedim, typename Number>
1202 template <typename ContainerType>
1203 void
1205 const ContainerType &cell_iterator_range,
1206 const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
1207 const unsigned int n_unfiltered_cells)
1208 {
1209 Assert(
1210 additional_data.use_global_weights == false,
1211 ExcMessage(
1212 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1213
1214 Assert(additional_data.store_cells == false, ExcNotImplemented());
1215
1216 if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
1217 update_flags_mapping |= update_covariant_transformation;
1218
1219 auto compute_mapping_data_for_surface =
1220 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1221 const ImmersedSurfaceQuadrature<dim> &quadrature,
1222 MappingData &mapping_data) {
1225 mapping,
1226 update_flags_mapping,
1227 cell,
1228 quadrature,
1229 internal_mapping_data,
1230 mapping_data);
1231 };
1232
1233 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1234 cell_iterator_range,
1235 quadrature_vector,
1236 n_unfiltered_cells,
1237 compute_mapping_data_for_surface);
1238 }
1239
1240
1241
1242 template <int dim, int spacedim, typename Number>
1243 template <typename ContainerType>
1244 void
1246 const ContainerType &cell_iterator_range,
1247 const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
1248 const unsigned int n_unfiltered_cells)
1249 {
1250 clear();
1251
1252 Assert(additional_data.store_cells == false, ExcNotImplemented());
1253
1254 do_cell_index_compression =
1255 n_unfiltered_cells != numbers::invalid_unsigned_int;
1256
1257 const unsigned int n_cells = quadrature_vector.size();
1258 AssertDimension(n_cells,
1259 std::distance(cell_iterator_range.begin(),
1260 cell_iterator_range.end()));
1261
1262 // fill cell index offset vector
1263 cell_index_offset.resize(n_cells);
1264 unsigned int n_faces = 0;
1265 unsigned int cell_index = 0;
1266 for (const auto &cell : cell_iterator_range)
1267 {
1268 cell_index_offset[cell_index] = n_faces;
1269 n_faces += cell->n_faces();
1270 ++cell_index;
1271 }
1272
1273 n_q_points_unvectorized.reserve(n_faces);
1274
1275 cell_type.reserve(n_faces);
1276
1277 // fill unit points index offset vector
1278 unit_points_index.resize(n_faces + 1);
1279 data_index_offsets.resize(n_faces + 1);
1280 cell_index = 0;
1281 unsigned int n_unit_points = 0;
1282 unsigned int n_data_points = 0;
1283 for (const auto &cell : cell_iterator_range)
1284 {
1285 for (const auto &f : cell->face_indices())
1286 {
1287 const unsigned int current_face_index =
1288 cell_index_offset[cell_index] + f;
1289
1290 unit_points_index[current_face_index] = n_unit_points;
1291 data_index_offsets[current_face_index] = n_data_points;
1292
1293 const unsigned int n_points =
1294 quadrature_vector[cell_index][f].size();
1295 n_q_points_unvectorized.push_back(n_points);
1296
1297 const unsigned int n_q_points =
1298 compute_n_q_points<VectorizedArrayType>(n_points);
1299 n_unit_points += n_q_points;
1300
1301 const unsigned int n_q_points_data =
1302 compute_n_q_points<Number>(n_points);
1303 n_data_points += n_q_points_data;
1304 }
1305
1306 ++cell_index;
1307 }
1308 unit_points_index[n_faces] = n_unit_points;
1309 data_index_offsets[n_faces] = n_data_points;
1310
1311 // compress indices
1312 if (do_cell_index_compression)
1313 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1315
1316 // fill unit points and mapping data for every face of all cells
1317 // resize data vectors
1318 resize_unit_points_faces(n_unit_points);
1319 resize_data_fields(n_data_points);
1320
1321 MappingData mapping_data;
1322 MappingData mapping_data_previous_cell;
1323 MappingData mapping_data_first;
1324 bool first_set = false;
1325 unsigned int size_compressed_data = 0;
1326 cell_index = 0;
1327 for (const auto &cell : cell_iterator_range)
1328 {
1329 const auto &quadratures_on_faces = quadrature_vector[cell_index];
1330
1331 Assert(quadratures_on_faces.size() == cell->n_faces(),
1332 ExcDimensionMismatch(quadratures_on_faces.size(),
1333 cell->n_faces()));
1334
1335 for (const auto &f : cell->face_indices())
1336 {
1337 const auto &quadrature_on_face = quadratures_on_faces[f];
1338 const bool empty = quadrature_on_face.empty();
1339
1340 const unsigned int current_face_index =
1341 cell_index_offset[cell_index] + f;
1342
1343 // store unit points
1344 const unsigned int n_q_points =
1345 compute_n_q_points<VectorizedArrayType>(
1346 n_q_points_unvectorized[current_face_index]);
1347 store_unit_points_faces(unit_points_index[current_face_index],
1348 n_q_points,
1349 n_q_points_unvectorized[current_face_index],
1350 quadrature_on_face.get_points());
1351
1354 update_flags_mapping,
1355 cell,
1356 f,
1357 quadrature_on_face,
1358 internal_mapping_data,
1359 mapping_data);
1360
1361 // check for cartesian/affine cell
1362 if (!empty &&
1363 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1364 {
1365 cell_type.push_back(internal::compute_geometry_type(
1366 cell->diameter(), mapping_data.inverse_jacobians));
1367
1368 if (!first_set)
1369 {
1370 mapping_data_first = mapping_data;
1371 first_set = true;
1372 }
1373 }
1374 else
1375 cell_type.push_back(
1377
1378 if (current_face_index > 0)
1379 {
1380 // check if current and previous cell are affine
1381 const bool affine_cells =
1382 cell_type[current_face_index] <=
1384 cell_type[current_face_index - 1] <=
1386
1387 // create a comparator to compare inverse Jacobian of current
1388 // and previous cell
1390 1e4 / cell->diameter() *
1391 std::numeric_limits<double>::epsilon() * 1024.);
1392
1393 // we can only compare if current and previous cell have at
1394 // least one quadrature point and both cells are at least affine
1395 const auto comparison_result =
1396 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1397 mapping_data_previous_cell.inverse_jacobians.empty()) ?
1399 comparator.compare(
1400 mapping_data.inverse_jacobians[0],
1401 mapping_data_previous_cell.inverse_jacobians[0]);
1402
1403 // we can compress the Jacobians and inverse Jacobians if
1404 // inverse Jacobians are equal and cells are affine
1405 if (affine_cells &&
1406 comparison_result ==
1408 {
1409 compressed_data_index_offsets.push_back(
1410 compressed_data_index_offsets.back());
1411 }
1412 else if (first_set &&
1413 (cell_type[current_face_index] <=
1415 (comparator.compare(
1416 mapping_data.inverse_jacobians[0],
1417 mapping_data_first.inverse_jacobians[0]) ==
1419 double>::ComparisonResult::equal))
1420 {
1421 compressed_data_index_offsets.push_back(0);
1422 }
1423 else
1424 {
1425 const unsigned int n_compressed_data_last_cell =
1426 cell_type[current_face_index - 1] <=
1428 1 :
1429 compute_n_q_points<Number>(
1430 n_q_points_unvectorized[current_face_index - 1]);
1431
1432 compressed_data_index_offsets.push_back(
1433 compressed_data_index_offsets.back() +
1434 n_compressed_data_last_cell);
1435 }
1436 }
1437 else
1438 compressed_data_index_offsets.push_back(0);
1439
1440 // cache mapping_data from previous cell
1441 mapping_data_previous_cell = mapping_data;
1442
1443 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1444 n_q_points_unvectorized[current_face_index]);
1445 store_mapping_data(data_index_offsets[current_face_index],
1446 n_q_points_data,
1447 n_q_points_unvectorized[current_face_index],
1448 mapping_data,
1449 quadrature_on_face.get_weights(),
1450 data_index_offsets[current_face_index],
1451 cell_type[current_face_index] <=
1453
1454 // update size of compressed data depending on cell type and handle
1455 // empty quadratures
1456 if (cell_type[current_face_index] <=
1458 size_compressed_data = compressed_data_index_offsets.back() + 1;
1459 else
1460 size_compressed_data =
1461 std::max(size_compressed_data,
1462 compressed_data_index_offsets.back() +
1463 n_q_points_data);
1464 }
1465 if (do_cell_index_compression)
1466 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1467 cell_index;
1468
1469 ++cell_index;
1470 }
1471
1472 if (update_flags_mapping & UpdateFlags::update_jacobians)
1473 {
1474 jacobians[0].resize(size_compressed_data);
1475 jacobians[0].shrink_to_fit();
1476 }
1477 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1478 {
1479 inverse_jacobians[0].resize(size_compressed_data);
1480 inverse_jacobians[0].shrink_to_fit();
1481 }
1482
1483 state = State::faces_on_cells_in_vector;
1484 is_reinitialized();
1485 }
1486
1487
1488
1489 template <int dim, int spacedim, typename Number>
1490 template <typename CellIteratorType>
1491 void
1493 const std::vector<std::pair<CellIteratorType, unsigned int>>
1494 &face_iterator_range_interior,
1495 const std::vector<Quadrature<dim - 1>> &quadrature_vector)
1496 {
1497 clear();
1498
1499 do_cell_index_compression = false;
1500
1501 Assert(additional_data.store_cells == false, ExcNotImplemented());
1502
1503
1504 const unsigned int n_faces = quadrature_vector.size();
1505 AssertDimension(n_faces,
1506 std::distance(face_iterator_range_interior.begin(),
1507 face_iterator_range_interior.end()));
1508
1509 n_q_points_unvectorized.reserve(n_faces);
1510
1511 cell_type.reserve(n_faces);
1512 face_number.reserve(n_faces);
1513
1514 // fill unit points index offset vector
1515 unit_points_index.reserve(n_faces + 1);
1516 unit_points_index.push_back(0);
1517 data_index_offsets.reserve(n_faces + 1);
1518 data_index_offsets.push_back(0);
1519 for (const auto &quadrature : quadrature_vector)
1520 {
1521 const unsigned int n_points = quadrature.size();
1522 n_q_points_unvectorized.push_back(n_points);
1523
1524 const unsigned int n_q_points =
1525 compute_n_q_points<VectorizedArrayType>(n_points);
1526 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1527
1528 const unsigned int n_q_points_data =
1529 compute_n_q_points<Number>(n_points);
1530 data_index_offsets.push_back(data_index_offsets.back() +
1531 n_q_points_data);
1532 }
1533
1534 const unsigned int n_unit_points = unit_points_index.back();
1535 const unsigned int n_data_points = data_index_offsets.back();
1536
1537 // resize data vectors
1538 resize_unit_points_faces(n_unit_points);
1539 resize_data_fields(n_data_points, true);
1540
1541 std::array<MappingData, 2> mapping_data;
1542 std::array<MappingData, 2> mapping_data_previous_cell;
1543 std::array<MappingData, 2> mapping_data_first;
1544 bool first_set = false;
1545 unsigned int size_compressed_data = 0;
1546 unsigned int face_index = 0;
1547 for (const auto &cell_and_f : face_iterator_range_interior)
1548 {
1549 const auto &quadrature_on_face = quadrature_vector[face_index];
1550 const bool empty = quadrature_on_face.empty();
1551
1552 // get interior cell and face number
1553 const auto &cell_m = cell_and_f.first;
1554 const auto f_m = cell_and_f.second;
1555
1556 // get exterior cell and face number
1557 const auto &cell_p =
1558 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1559 const auto f_p =
1560 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1561
1562 face_number.emplace_back(f_m, f_p);
1563
1564 Assert(
1565 cell_m->combined_face_orientation(f_m) ==
1567 cell_p->combined_face_orientation(f_p) ==
1569 ExcMessage(
1570 "Non standard face orientation is currently not implemented."));
1571
1572 // store unit points
1573 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1574 n_q_points_unvectorized[face_index]);
1575 store_unit_points_faces(unit_points_index[face_index],
1576 n_q_points,
1577 n_q_points_unvectorized[face_index],
1578 quadrature_on_face.get_points());
1579
1580 // compute mapping for interior face
1583 update_flags_mapping,
1584 cell_m,
1585 f_m,
1586 quadrature_on_face,
1587 internal_mapping_data,
1588 mapping_data[0]);
1589
1590 // compute mapping for exterior face
1593 update_flags_mapping,
1594 cell_p,
1595 f_p,
1596 quadrature_on_face,
1597 internal_mapping_data,
1598 mapping_data[1]);
1599
1600 // check for cartesian/affine cell
1601 if (!empty &&
1602 update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1603 {
1604 // select more general type of interior and exterior cell
1605 cell_type.push_back(std::max(
1607 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1609 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1610
1611 // cache mapping data of first cell pair with non-empty quadrature
1612 // on the face
1613 if (!first_set)
1614 {
1615 mapping_data_first = mapping_data;
1616 first_set = true;
1617 }
1618 }
1619 else
1620 cell_type.push_back(
1622
1623 if (face_index > 0)
1624 {
1625 // check if current and previous cell pairs are affine
1626 const bool affine_cells =
1627 cell_type[face_index] <=
1629 cell_type[face_index - 1] <=
1631
1632 // create a comparator to compare inverse Jacobian of current
1633 // and previous cell pair
1635 1e4 / cell_m->diameter() *
1636 std::numeric_limits<double>::epsilon() * 1024.);
1637
1638 // we can only compare if current and previous cell have at
1639 // least one quadrature point and both cells are at least affine
1640 const auto comparison_result_m =
1641 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1642 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1644 comparator.compare(
1645 mapping_data[0].inverse_jacobians[0],
1646 mapping_data_previous_cell[0].inverse_jacobians[0]);
1647
1648 const auto comparison_result_p =
1649 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1650 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1652 comparator.compare(
1653 mapping_data[1].inverse_jacobians[0],
1654 mapping_data_previous_cell[1].inverse_jacobians[0]);
1655
1656 // we can compress the Jacobians and inverse Jacobians if
1657 // inverse Jacobians are equal and cells are affine
1658 if (affine_cells &&
1659 comparison_result_m ==
1661 comparison_result_p ==
1663 {
1664 compressed_data_index_offsets.push_back(
1665 compressed_data_index_offsets.back());
1666 }
1667 else if (first_set &&
1668 (cell_type[face_index] <=
1670 (comparator.compare(
1671 mapping_data[0].inverse_jacobians[0],
1672 mapping_data_first[0].inverse_jacobians[0]) ==
1674 double>::ComparisonResult::equal) &&
1675 (comparator.compare(
1676 mapping_data[1].inverse_jacobians[0],
1677 mapping_data_first[1].inverse_jacobians[0]) ==
1679 {
1680 compressed_data_index_offsets.push_back(0);
1681 }
1682 else
1683 {
1684 const unsigned int n_compressed_data_last_cell =
1685 cell_type[face_index - 1] <=
1687 1 :
1688 compute_n_q_points<Number>(
1689 n_q_points_unvectorized[face_index - 1]);
1690
1691 compressed_data_index_offsets.push_back(
1692 compressed_data_index_offsets.back() +
1693 n_compressed_data_last_cell);
1694 }
1695 }
1696 else
1697 compressed_data_index_offsets.push_back(0);
1698
1699 // cache mapping_data from previous cell pair
1700 mapping_data_previous_cell = mapping_data;
1701
1702 const unsigned int n_q_points_data =
1703 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1704
1705 // store mapping data of interior face
1706 store_mapping_data(data_index_offsets[face_index],
1707 n_q_points_data,
1708 n_q_points_unvectorized[face_index],
1709 mapping_data[0],
1710 quadrature_on_face.get_weights(),
1711 data_index_offsets[face_index],
1712 cell_type[face_index] <=
1714 true);
1715
1716 // store only necessary mapping data for exterior face (Jacobians and
1717 // inverse Jacobians)
1718 store_mapping_data(data_index_offsets[face_index],
1719 n_q_points_data,
1720 n_q_points_unvectorized[face_index],
1721 mapping_data[1],
1722 quadrature_on_face.get_weights(),
1723 data_index_offsets[face_index],
1724 cell_type[face_index] <=
1726 false);
1727
1728 // update size of compressed data depending on cell type and handle
1729 // empty quadratures
1730 if (cell_type[face_index] <=
1732 size_compressed_data = compressed_data_index_offsets.back() + 1;
1733 else
1734 size_compressed_data =
1735 std::max(size_compressed_data,
1736 compressed_data_index_offsets.back() + n_q_points_data);
1737
1738 ++face_index;
1739 }
1740
1741 if (update_flags_mapping & UpdateFlags::update_jacobians)
1742 {
1743 jacobians[0].resize(size_compressed_data);
1744 jacobians[0].shrink_to_fit();
1745 jacobians[1].resize(size_compressed_data);
1746 jacobians[1].shrink_to_fit();
1747 }
1748 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1749 {
1750 inverse_jacobians[0].resize(size_compressed_data);
1751 inverse_jacobians[0].shrink_to_fit();
1752 inverse_jacobians[1].resize(size_compressed_data);
1753 inverse_jacobians[1].shrink_to_fit();
1754 }
1755
1756 state = State::face_vector;
1757 is_reinitialized();
1758 }
1759
1760
1761
1762 template <int dim, int spacedim, typename Number>
1763 bool
1765 {
1766 return state == State::faces_on_cells_in_vector;
1767 }
1768
1769
1770
1771 template <int dim, int spacedim, typename Number>
1772 unsigned int
1774 const unsigned int geometry_index) const
1775 {
1776 return n_q_points_unvectorized[geometry_index];
1777 }
1778
1779
1780 template <int dim, int spacedim, typename Number>
1783 const unsigned int geometry_index) const
1784 {
1785 return cell_type[geometry_index];
1786 }
1787
1788
1789
1790 template <int dim, int spacedim, typename Number>
1793 const unsigned int cell_index) const
1794 {
1795 Assert(
1796 additional_data.store_cells,
1797 ExcMessage(
1798 "Cells have been not stored. You can enable this by Additional::store_cells."));
1799 return {triangulation.get(),
1800 cell_level_and_indices[cell_index].first,
1801 cell_level_and_indices[cell_index].second};
1802 }
1803
1804
1805
1806 template <int dim, int spacedim, typename Number>
1807 template <typename NumberType>
1808 unsigned int
1810 const unsigned int n_q_points_unvectorized)
1811 {
1812 const unsigned int n_lanes =
1814 const unsigned int n_filled_lanes_last_batch =
1815 n_q_points_unvectorized % n_lanes;
1816 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1817 if (n_filled_lanes_last_batch > 0)
1818 ++n_q_points;
1819 return n_q_points;
1820 }
1821
1822
1823
1824 template <int dim, int spacedim, typename Number>
1825 template <bool is_face>
1826 unsigned int
1828 const unsigned int cell_index,
1829 const unsigned int face_number) const
1830 {
1832 return 0;
1833
1834 const unsigned int compressed_cell_index =
1835 compute_compressed_cell_index(cell_index);
1836 if (!is_face)
1837 {
1838 Assert(state == State::cell_vector,
1839 ExcMessage(
1840 "This mapping info is not reinitialized for a cell vector!"));
1841 return compressed_cell_index;
1842 }
1843 else
1844 {
1846 ExcMessage(
1847 "cell_index has to be set if face_number is specified!"));
1848 Assert(state == State::faces_on_cells_in_vector ||
1849 state == State::face_vector,
1850 ExcMessage("This mapping info is not reinitialized for faces"
1851 " on cells in a vector!"));
1852 if (state == State::faces_on_cells_in_vector)
1853 return cell_index_offset[compressed_cell_index] + face_number;
1854 else if (state == State::face_vector)
1855 return cell_index;
1856 }
1857 }
1858
1859
1860
1861 template <int dim, int spacedim, typename Number>
1862 unsigned int
1864 const unsigned int cell_index) const
1865 {
1866 if (do_cell_index_compression)
1867 {
1868 Assert(cell_index_to_compressed_cell_index[cell_index] !=
1870 ExcMessage("Mapping info object was not initialized for this"
1871 " active cell index!"));
1872 return cell_index_to_compressed_cell_index[cell_index];
1873 }
1874 else
1875 return cell_index;
1876 }
1877
1878
1879 template <int dim, int spacedim, typename Number>
1880 void
1882 const unsigned int unit_points_index_offset,
1883 const unsigned int n_q_points,
1884 const unsigned int n_q_points_unvectorized,
1885 const std::vector<Point<dim>> &points)
1886 {
1887 const unsigned int n_lanes =
1889
1890 for (unsigned int q = 0; q < n_q_points; ++q)
1891 {
1892 const unsigned int offset = unit_points_index_offset + q;
1893 for (unsigned int v = 0;
1894 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1895 ++v)
1896 for (unsigned int d = 0; d < dim; ++d)
1898 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1899 }
1900 }
1901
1902
1903
1904 template <int dim, int spacedim, typename Number>
1905 void
1907 const unsigned int unit_points_index_offset,
1908 const unsigned int n_q_points,
1909 const unsigned int n_q_points_unvectorized,
1910 const std::vector<Point<dim - 1>> &points)
1911 {
1912 const unsigned int n_lanes =
1914
1915 for (unsigned int q = 0; q < n_q_points; ++q)
1916 {
1917 const unsigned int offset = unit_points_index_offset + q;
1918 for (unsigned int v = 0;
1919 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1920 ++v)
1921 for (unsigned int d = 0; d < dim - 1; ++d)
1923 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1924 }
1925 }
1926
1927
1928
1929 template <int dim, int spacedim, typename Number>
1930 void
1932 const unsigned int unit_points_index_offset,
1933 const unsigned int n_q_points,
1934 const unsigned int n_q_points_unvectorized,
1935 const MappingInfo::MappingData &mapping_data,
1936 const std::vector<double> &weights,
1937 const unsigned int compressed_unit_point_index_offset,
1938 const bool affine_cell,
1939 const bool is_interior)
1940 {
1941 const unsigned int n_lanes =
1943
1944 for (unsigned int q = 0; q < n_q_points; ++q)
1945 {
1946 const unsigned int offset = unit_points_index_offset + q;
1947 const unsigned int compressed_offset =
1948 compressed_unit_point_index_offset + q;
1949 for (unsigned int v = 0;
1950 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1951 ++v)
1952 {
1953 if (q == 0 || !affine_cell)
1954 {
1955 if (update_flags_mapping & UpdateFlags::update_jacobians)
1956 for (unsigned int d = 0; d < dim; ++d)
1957 for (unsigned int s = 0; s < spacedim; ++s)
1959 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1960 v) = mapping_data.jacobians[q * n_lanes + v][d][s];
1961 if (update_flags_mapping &
1963 for (unsigned int d = 0; d < dim; ++d)
1964 for (unsigned int s = 0; s < spacedim; ++s)
1966 inverse_jacobians[is_interior ? 0 : 1]
1967 [compressed_offset][d][s],
1968 v) =
1969 mapping_data.inverse_jacobians[q * n_lanes + v][d][s];
1970 }
1971
1972 if (is_interior)
1973 {
1974 if (update_flags_mapping & UpdateFlags::update_JxW_values)
1975 {
1976 if (additional_data.use_global_weights)
1977 {
1979 JxW_values[offset], v) = weights[q * n_lanes + v];
1980 }
1981 else
1982 {
1984 JxW_values[offset], v) =
1985 mapping_data.JxW_values[q * n_lanes + v];
1986 }
1987 }
1988 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
1989 for (unsigned int s = 0; s < spacedim; ++s)
1991 normal_vectors[offset][s], v) =
1992 mapping_data.normal_vectors[q * n_lanes + v][s];
1993 if (update_flags_mapping &
1995 for (unsigned int s = 0; s < spacedim; ++s)
1997 real_points[offset][s], v) =
1998 mapping_data.quadrature_points[q * n_lanes + v][s];
1999 }
2000 }
2001 }
2002 }
2003
2004
2005
2006 template <int dim, int spacedim, typename Number>
2007 void
2009 const unsigned int n_unit_point_batches)
2010 {
2011 unit_points.resize(n_unit_point_batches);
2012 }
2013
2014
2015
2016 template <int dim, int spacedim, typename Number>
2017 void
2019 const unsigned int n_unit_point_batches)
2020 {
2021 unit_points_faces.resize(n_unit_point_batches);
2022 }
2023
2024
2025
2026 template <int dim, int spacedim, typename Number>
2027 void
2029 const unsigned int n_data_point_batches,
2030 const bool is_face_centric)
2031 {
2032 if (update_flags_mapping & UpdateFlags::update_jacobians)
2033 {
2034 jacobians[0].resize(n_data_point_batches);
2035 if (is_face_centric)
2036 jacobians[1].resize(n_data_point_batches);
2037 }
2038 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2039 {
2040 inverse_jacobians[0].resize(n_data_point_batches);
2041 if (is_face_centric)
2042 inverse_jacobians[1].resize(n_data_point_batches);
2043 }
2044 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2045 JxW_values.resize(n_data_point_batches);
2046 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2047 normal_vectors.resize(n_data_point_batches);
2048 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2049 real_points.resize(n_data_point_batches);
2050 }
2051
2052
2053
2054 template <int dim, int spacedim, typename Number>
2055 inline const Point<
2056 dim,
2059 const unsigned int offset) const
2060 {
2061 return unit_points.data() + offset;
2062 }
2063
2064
2065
2066 template <int dim, int spacedim, typename Number>
2067 inline const Point<
2068 dim - 1,
2071 const unsigned int offset) const
2072 {
2073 return unit_points_faces.data() + offset;
2074 }
2075
2076
2077
2078 template <int dim, int spacedim, typename Number>
2079 inline const Point<spacedim, Number> *
2081 const unsigned int offset) const
2082 {
2083 return real_points.data() + offset;
2084 }
2085
2086
2087
2088 template <int dim, int spacedim, typename Number>
2089 unsigned int
2091 const unsigned int geometry_index) const
2092 {
2093 return unit_points_index[geometry_index];
2094 }
2095
2096
2097
2098 template <int dim, int spacedim, typename Number>
2099 unsigned int
2101 const unsigned int geometry_index) const
2102 {
2103 return data_index_offsets[geometry_index];
2104 }
2105
2106
2107 template <int dim, int spacedim, typename Number>
2108 unsigned int
2110 const unsigned int geometry_index) const
2111 {
2112 return compressed_data_index_offsets[geometry_index];
2113 }
2114
2115
2116
2117 template <int dim, int spacedim, typename Number>
2120 const bool is_interior) const
2121 {
2122 return jacobians[is_interior ? 0 : 1].data() + offset;
2123 }
2124
2125
2126
2127 template <int dim, int spacedim, typename Number>
2130 const unsigned int offset,
2131 const bool is_interior) const
2132 {
2133 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2134 }
2135
2136
2137
2138 template <int dim, int spacedim, typename Number>
2139 inline const Tensor<1, spacedim, Number> *
2141 const unsigned int offset) const
2142 {
2143 return normal_vectors.data() + offset;
2144 }
2145
2146
2147
2148 template <int dim, int spacedim, typename Number>
2149 inline const Number *
2150 MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
2151 {
2152 return JxW_values.data() + offset;
2153 }
2154
2155
2156
2157 template <int dim, int spacedim, typename Number>
2160 {
2161 return *mapping;
2162 }
2163
2164
2165
2166 template <int dim, int spacedim, typename Number>
2169 {
2170 return update_flags;
2171 }
2172
2173
2174
2175 template <int dim, int spacedim, typename Number>
2178 {
2179 return update_flags_mapping;
2180 }
2181
2182
2183
2184 template <int dim, int spacedim, typename Number>
2185 boost::signals2::connection
2187 const std::function<void()> &set_is_reinitialized)
2188 {
2189 return is_reinitialized.connect(set_is_reinitialized);
2190 }
2191
2192
2193
2194 template <int dim, int spacedim, typename Number>
2195 std::size_t
2197 {
2198 std::size_t memory = MemoryConsumption::memory_consumption(unit_points);
2199 memory += MemoryConsumption::memory_consumption(unit_points_faces);
2200 memory += MemoryConsumption::memory_consumption(unit_points_index);
2201 memory += cell_type.capacity() *
2203 memory += MemoryConsumption::memory_consumption(data_index_offsets);
2204 memory +=
2205 MemoryConsumption::memory_consumption(compressed_data_index_offsets);
2206 memory += MemoryConsumption::memory_consumption(JxW_values);
2207 memory += MemoryConsumption::memory_consumption(normal_vectors);
2208 memory += MemoryConsumption::memory_consumption(jacobians);
2209 memory += MemoryConsumption::memory_consumption(inverse_jacobians);
2210 memory += MemoryConsumption::memory_consumption(real_points);
2211 memory += MemoryConsumption::memory_consumption(n_q_points_unvectorized);
2212 memory += MemoryConsumption::memory_consumption(cell_index_offset);
2214 cell_index_to_compressed_cell_index);
2215 memory += MemoryConsumption::memory_consumption(cell_level_and_indices);
2216 memory += sizeof(*this);
2217 return memory;
2218 }
2219} // namespace NonMatching
2220
2222
2223#endif
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:161
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:81
Abstract base class for mapping classes.
Definition mapping.h:316
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
std::shared_ptr< typename Mapping< dim, spacedim >::InternalDataBase > internal_mapping_data
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
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)
boost::signals2::signal< void()> is_reinitialized
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)
void reinit_faces(const std::vector< std::pair< CellIteratorType, unsigned int > > &face_iterator_range_interior, const std::vector< Quadrature< dim - 1 > > &quadrature_vector)
UpdateFlags update_flags_mapping
void do_reinit_cells(const ContainerType &cell_iterator_range, const std::vector< QuadratureType > &quadrature_vector, const unsigned int n_unfiltered_cells, const std::function< void(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const QuadratureType &quadrature, MappingData &mapping_data)> &compute_mapping_data)
UpdateFlags get_update_flags_mapping() const
Triangulation< dim, spacedim >::cell_iterator get_cell_iterator(const unsigned int cell_index) const
static void compute_mapping_data_for_immersed_surface_quadrature(const SmartPointer< const Mapping< dim, spacedim > > mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ImmersedSurfaceQuadrature< dim > &quadrature, std::shared_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::shared_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::shared_ptr< typename Mapping< dim, spacedim >::InternalDataBase > internal_mapping_data, MappingData &mapping_data)
static UpdateFlags required_update_flags(const SmartPointer< const Mapping< dim, spacedim > > mapping, const UpdateFlags &update_flags)
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)