16#ifndef dealii_non_matching_mapping_info_h
17#define dealii_non_matching_mapping_info_h
45 template <
int dim,
int spacedim = dim>
58 return mapping->requires_update_flags(update_flags);
69 internal_mapping_data,
81 *internal_mapping_data);
88 internal_mapping_data =
89 mapping->get_data(update_flags_mapping, quadrature);
92 cell_similarity = mapping->fill_fe_values(cell,
95 *internal_mapping_data,
108 internal_mapping_data,
118 *internal_mapping_data);
125 internal_mapping_data =
126 mapping->get_data(update_flags_mapping, quadrature);
129 mapping->fill_fe_immersed_surface_values(cell,
131 *internal_mapping_data,
142 const unsigned int face_no,
145 internal_mapping_data,
156 *internal_mapping_data);
159 ReferenceCells::get_hypercube<dim>(),
162 cell->face_orientation(face_no),
163 cell->face_flip(face_no),
164 cell->face_rotation(face_no)),
167 mapping_q->fill_mapping_data_for_face_quadrature(
168 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
172 auto internal_mapping_data =
173 mapping->get_face_data(update_flags_mapping,
176 mapping->fill_fe_face_values(cell,
179 *internal_mapping_data,
185 template <
int dim,
int spacedim = dim>
188 const double diameter,
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)
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)
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)
213 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
215 cell_cartesian =
false;
221 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
222 else if (jacobian_constant)
223 return ::internal::MatrixFreeFunctions::GeometryType::affine;
225 return ::internal::MatrixFreeFunctions::GeometryType::general;
245 template <
int dim,
int spacedim = dim,
typename Number =
double>
253 Number>::vectorized_value_type;
338 template <
typename ContainerType>
341 const ContainerType &cell_iterator_range,
342 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
351 template <
typename ContainerType>
354 const ContainerType &cell_iterator_range,
362 template <
typename ContainerType>
365 const ContainerType &cell_iterator_range,
373 template <
typename ContainerType>
376 const ContainerType &cell_iterator_range,
384 template <
typename CellIteratorType>
386 reinit_faces(
const std::vector<std::pair<CellIteratorType, unsigned int>>
387 &face_iterator_range_interior,
423 const bool is_interior =
true)
const;
431 const bool is_interior =
true)
const;
475 boost::signals2::connection
481 template <
bool is_face>
503 const unsigned int geometry_index)
const;
539 template <
typename NumberType>
566 const bool is_face_centric =
false);
573 const unsigned int n_q_points,
582 const unsigned int n_q_points,
591 const unsigned int n_q_points,
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);
608 template <
typename ContainerType,
typename QuadratureType>
611 const ContainerType &cell_iterator_range,
612 const std::vector<QuadratureType> &quadrature_vector,
613 const unsigned int n_unfiltered_cells,
616 const QuadratureType &quadrature,
617 MappingData &mapping_data)> &compute_mapping_data);
659 std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
688 std::vector<::internal::MatrixFreeFunctions::GeometryType>
cell_type;
721 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
731 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
785 template <
int dim,
int spacedim,
typename Number>
788 const unsigned int offset,
789 const bool is_interior)
const
791 const auto &face_pair = face_number[offset];
792 return is_interior ? face_pair.first : face_pair.second;
798 template <
int dim,
int spacedim,
typename Number>
804 , update_flags(update_flags)
806 , additional_data(additional_data)
832 std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
833 mapping_q->get_degree());
839 template <
int dim,
int spacedim,
typename Number>
843 n_q_points_unvectorized.clear();
844 unit_points_index.clear();
845 data_index_offsets.clear();
846 compressed_data_index_offsets.clear();
852 template <
int dim,
int spacedim,
typename Number>
856 const std::vector<
Point<dim>> &unit_points_in)
863 template <
int dim,
int spacedim,
typename Number>
870 std::vector<
Point<dim>>(unit_points_in.begin(),
871 unit_points_in.end()));
876 template <
int dim,
int spacedim,
typename Number>
882 n_q_points_unvectorized.resize(1);
883 n_q_points_unvectorized[0] = quadrature.
size();
885 const unsigned int n_q_points =
886 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
888 const unsigned int n_q_points_data =
889 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
892 resize_unit_points(n_q_points);
893 resize_data_fields(n_q_points_data);
898 n_q_points_unvectorized[0],
906 update_flags_mapping,
910 internal_mapping_data,
914 if (!quadrature.
empty() &&
929 n_q_points_unvectorized[0],
936 unit_points_index.push_back(0);
937 data_index_offsets.push_back(0);
938 compressed_data_index_offsets.push_back(0);
940 state = State::single_cell;
946 template <
int dim,
int spacedim,
typename Number>
947 template <
typename ContainerType>
950 const ContainerType &cell_iterator_range,
951 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
952 const unsigned int n_unfiltered_cells)
954 const unsigned int n_cells = unit_points_vector.size();
956 std::distance(cell_iterator_range.begin(),
957 cell_iterator_range.end()));
959 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
964 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
969 template <
int dim,
int spacedim,
typename Number>
970 template <
typename ContainerType,
typename QuadratureType>
973 const ContainerType &cell_iterator_range,
974 const std::vector<QuadratureType> &quadrature_vector,
975 const unsigned int n_unfiltered_cells,
978 const QuadratureType &quadrature,
983 do_cell_index_compression =
986 const unsigned int n_cells = quadrature_vector.size();
988 std::distance(cell_iterator_range.begin(),
989 cell_iterator_range.end()));
991 n_q_points_unvectorized.reserve(n_cells);
993 cell_type.reserve(n_cells);
995 if (additional_data.store_cells)
996 cell_level_and_indices.resize(n_cells);
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)
1005 const unsigned int n_points = quadrature.size();
1006 n_q_points_unvectorized.push_back(n_points);
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);
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() +
1018 const unsigned int n_unit_points = unit_points_index.back();
1019 const unsigned int n_data_points = data_index_offsets.back();
1022 resize_unit_points(n_unit_points);
1023 resize_data_fields(n_data_points);
1025 if (do_cell_index_compression)
1026 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1031 unsigned int size_compressed_data = 0;
1033 for (
const auto &cell : cell_iterator_range)
1035 if (additional_data.store_cells)
1038 cell_level_and_indices[
cell_index] = {cell->level(), cell->index()};
1041 const auto &quadrature = quadrature_vector[
cell_index];
1042 const bool empty = quadrature.empty();
1045 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1047 store_unit_points(unit_points_index[
cell_index],
1050 quadrature.get_points());
1053 compute_mapping_data(cell, quadrature, mapping_data);
1056 const unsigned int n_q_points_data =
1057 compute_n_q_points<Number>(n_q_points_unvectorized[
cell_index]);
1063 cell_type.push_back(
1068 cell_type.push_back(
1074 const bool affine_cells =
1083 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1088 const auto comparison_result =
1099 comparison_result ==
1102 compressed_data_index_offsets.push_back(
1103 compressed_data_index_offsets.back());
1107 const unsigned int n_compressed_data_last_cell =
1111 compute_n_q_points<Number>(
1114 compressed_data_index_offsets.push_back(
1115 compressed_data_index_offsets.back() +
1116 n_compressed_data_last_cell);
1120 compressed_data_index_offsets.push_back(0);
1123 mapping_data_previous_cell = mapping_data;
1125 store_mapping_data(data_index_offsets[
cell_index],
1129 quadrature.get_weights(),
1138 size_compressed_data = compressed_data_index_offsets.back() + 1;
1140 size_compressed_data =
1142 compressed_data_index_offsets.back() + n_q_points_data);
1144 if (do_cell_index_compression)
1145 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1153 jacobians[0].resize(size_compressed_data);
1154 jacobians[0].shrink_to_fit();
1158 inverse_jacobians[0].resize(size_compressed_data);
1159 inverse_jacobians[0].shrink_to_fit();
1162 state = State::cell_vector;
1168 template <
int dim,
int spacedim,
typename Number>
1169 template <
typename ContainerType>
1172 const ContainerType &cell_iterator_range,
1174 const unsigned int n_unfiltered_cells)
1176 auto compute_mapping_data_for_cells =
1184 update_flags_mapping,
1188 internal_mapping_data,
1192 do_reinit_cells<ContainerType, Quadrature<dim>>(
1193 cell_iterator_range,
1196 compute_mapping_data_for_cells);
1201 template <
int dim,
int spacedim,
typename Number>
1202 template <
typename ContainerType>
1205 const ContainerType &cell_iterator_range,
1207 const unsigned int n_unfiltered_cells)
1210 additional_data.use_global_weights ==
false,
1212 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1219 auto compute_mapping_data_for_surface =
1226 update_flags_mapping,
1229 internal_mapping_data,
1233 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1234 cell_iterator_range,
1237 compute_mapping_data_for_surface);
1242 template <
int dim,
int spacedim,
typename Number>
1243 template <
typename ContainerType>
1246 const ContainerType &cell_iterator_range,
1248 const unsigned int n_unfiltered_cells)
1254 do_cell_index_compression =
1257 const unsigned int n_cells = quadrature_vector.size();
1259 std::distance(cell_iterator_range.begin(),
1260 cell_iterator_range.end()));
1263 cell_index_offset.resize(n_cells);
1264 unsigned int n_faces = 0;
1266 for (
const auto &cell : cell_iterator_range)
1269 n_faces += cell->n_faces();
1273 n_q_points_unvectorized.reserve(n_faces);
1275 cell_type.reserve(n_faces);
1278 unit_points_index.resize(n_faces + 1);
1279 data_index_offsets.resize(n_faces + 1);
1281 unsigned int n_unit_points = 0;
1282 unsigned int n_data_points = 0;
1283 for (
const auto &cell : cell_iterator_range)
1285 for (
const auto &f : cell->face_indices())
1287 const unsigned int current_face_index =
1290 unit_points_index[current_face_index] = n_unit_points;
1291 data_index_offsets[current_face_index] = n_data_points;
1293 const unsigned int n_points =
1295 n_q_points_unvectorized.push_back(n_points);
1297 const unsigned int n_q_points =
1298 compute_n_q_points<VectorizedArrayType>(n_points);
1299 n_unit_points += n_q_points;
1301 const unsigned int n_q_points_data =
1302 compute_n_q_points<Number>(n_points);
1303 n_data_points += n_q_points_data;
1308 unit_points_index[n_faces] = n_unit_points;
1309 data_index_offsets[n_faces] = n_data_points;
1312 if (do_cell_index_compression)
1313 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1318 resize_unit_points_faces(n_unit_points);
1319 resize_data_fields(n_data_points);
1324 bool first_set =
false;
1325 unsigned int size_compressed_data = 0;
1327 for (
const auto &cell : cell_iterator_range)
1329 const auto &quadratures_on_faces = quadrature_vector[
cell_index];
1331 Assert(quadratures_on_faces.size() == cell->n_faces(),
1335 for (
const auto &f : cell->face_indices())
1337 const auto &quadrature_on_face = quadratures_on_faces[f];
1338 const bool empty = quadrature_on_face.empty();
1340 const unsigned int current_face_index =
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],
1349 n_q_points_unvectorized[current_face_index],
1350 quadrature_on_face.get_points());
1354 update_flags_mapping,
1358 internal_mapping_data,
1370 mapping_data_first = mapping_data;
1375 cell_type.push_back(
1378 if (current_face_index > 0)
1381 const bool affine_cells =
1382 cell_type[current_face_index] <=
1384 cell_type[current_face_index - 1] <=
1390 1e4 / cell->diameter() *
1391 std::numeric_limits<double>::epsilon() * 1024.);
1395 const auto comparison_result =
1406 comparison_result ==
1409 compressed_data_index_offsets.push_back(
1410 compressed_data_index_offsets.back());
1412 else if (first_set &&
1413 (cell_type[current_face_index] <=
1419 double>::ComparisonResult::equal))
1421 compressed_data_index_offsets.push_back(0);
1425 const unsigned int n_compressed_data_last_cell =
1426 cell_type[current_face_index - 1] <=
1429 compute_n_q_points<Number>(
1430 n_q_points_unvectorized[current_face_index - 1]);
1432 compressed_data_index_offsets.push_back(
1433 compressed_data_index_offsets.back() +
1434 n_compressed_data_last_cell);
1438 compressed_data_index_offsets.push_back(0);
1441 mapping_data_previous_cell = mapping_data;
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],
1447 n_q_points_unvectorized[current_face_index],
1449 quadrature_on_face.get_weights(),
1450 data_index_offsets[current_face_index],
1451 cell_type[current_face_index] <=
1456 if (cell_type[current_face_index] <=
1458 size_compressed_data = compressed_data_index_offsets.back() + 1;
1460 size_compressed_data =
1462 compressed_data_index_offsets.back() +
1465 if (do_cell_index_compression)
1466 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1474 jacobians[0].resize(size_compressed_data);
1475 jacobians[0].shrink_to_fit();
1479 inverse_jacobians[0].resize(size_compressed_data);
1480 inverse_jacobians[0].shrink_to_fit();
1483 state = State::faces_on_cells_in_vector;
1489 template <
int dim,
int spacedim,
typename Number>
1490 template <
typename CellIteratorType>
1493 const std::vector<std::pair<CellIteratorType, unsigned int>>
1494 &face_iterator_range_interior,
1499 do_cell_index_compression =
false;
1504 const unsigned int n_faces = quadrature_vector.size();
1506 std::distance(face_iterator_range_interior.begin(),
1507 face_iterator_range_interior.end()));
1509 n_q_points_unvectorized.reserve(n_faces);
1511 cell_type.reserve(n_faces);
1512 face_number.reserve(n_faces);
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)
1521 const unsigned int n_points = quadrature.size();
1522 n_q_points_unvectorized.push_back(n_points);
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);
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() +
1534 const unsigned int n_unit_points = unit_points_index.back();
1535 const unsigned int n_data_points = data_index_offsets.back();
1538 resize_unit_points_faces(n_unit_points);
1539 resize_data_fields(n_data_points,
true);
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)
1549 const auto &quadrature_on_face = quadrature_vector[face_index];
1550 const bool empty = quadrature_on_face.empty();
1553 const auto &cell_m = cell_and_f.first;
1554 const auto f_m = cell_and_f.second;
1557 const auto &cell_p =
1558 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1560 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1562 face_number.emplace_back(f_m, f_p);
1565 cell_m->combined_face_orientation(f_m) ==
1567 cell_p->combined_face_orientation(f_p) ==
1570 "Non standard face orientation is currently not implemented."));
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],
1577 n_q_points_unvectorized[face_index],
1578 quadrature_on_face.get_points());
1583 update_flags_mapping,
1587 internal_mapping_data,
1593 update_flags_mapping,
1597 internal_mapping_data,
1607 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1609 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1615 mapping_data_first = mapping_data;
1620 cell_type.push_back(
1626 const bool affine_cells =
1627 cell_type[face_index] <=
1629 cell_type[face_index - 1] <=
1635 1e4 / cell_m->diameter() *
1636 std::numeric_limits<double>::epsilon() * 1024.);
1640 const auto comparison_result_m =
1641 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1642 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1645 mapping_data[0].inverse_jacobians[0],
1646 mapping_data_previous_cell[0].inverse_jacobians[0]);
1648 const auto comparison_result_p =
1649 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1650 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1653 mapping_data[1].inverse_jacobians[0],
1654 mapping_data_previous_cell[1].inverse_jacobians[0]);
1659 comparison_result_m ==
1661 comparison_result_p ==
1664 compressed_data_index_offsets.push_back(
1665 compressed_data_index_offsets.back());
1667 else if (first_set &&
1668 (cell_type[face_index] <=
1671 mapping_data[0].inverse_jacobians[0],
1672 mapping_data_first[0].inverse_jacobians[0]) ==
1674 double>::ComparisonResult::equal) &&
1676 mapping_data[1].inverse_jacobians[0],
1677 mapping_data_first[1].inverse_jacobians[0]) ==
1680 compressed_data_index_offsets.push_back(0);
1684 const unsigned int n_compressed_data_last_cell =
1685 cell_type[face_index - 1] <=
1688 compute_n_q_points<Number>(
1689 n_q_points_unvectorized[face_index - 1]);
1691 compressed_data_index_offsets.push_back(
1692 compressed_data_index_offsets.back() +
1693 n_compressed_data_last_cell);
1697 compressed_data_index_offsets.push_back(0);
1700 mapping_data_previous_cell = mapping_data;
1702 const unsigned int n_q_points_data =
1703 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1706 store_mapping_data(data_index_offsets[face_index],
1708 n_q_points_unvectorized[face_index],
1710 quadrature_on_face.get_weights(),
1711 data_index_offsets[face_index],
1712 cell_type[face_index] <=
1718 store_mapping_data(data_index_offsets[face_index],
1720 n_q_points_unvectorized[face_index],
1722 quadrature_on_face.get_weights(),
1723 data_index_offsets[face_index],
1724 cell_type[face_index] <=
1730 if (cell_type[face_index] <=
1732 size_compressed_data = compressed_data_index_offsets.back() + 1;
1734 size_compressed_data =
1736 compressed_data_index_offsets.back() + n_q_points_data);
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();
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();
1756 state = State::face_vector;
1762 template <
int dim,
int spacedim,
typename Number>
1766 return state == State::faces_on_cells_in_vector;
1771 template <
int dim,
int spacedim,
typename Number>
1774 const unsigned int geometry_index)
const
1776 return n_q_points_unvectorized[geometry_index];
1780 template <
int dim,
int spacedim,
typename Number>
1783 const unsigned int geometry_index)
const
1785 return cell_type[geometry_index];
1790 template <
int dim,
int spacedim,
typename Number>
1796 additional_data.store_cells,
1798 "Cells have been not stored. You can enable this by Additional::store_cells."));
1806 template <
int dim,
int spacedim,
typename Number>
1807 template <
typename NumberType>
1810 const unsigned int n_q_points_unvectorized)
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)
1824 template <
int dim,
int spacedim,
typename Number>
1825 template <
bool is_face>
1829 const unsigned int face_number)
const
1834 const unsigned int compressed_cell_index =
1838 Assert(state == State::cell_vector,
1840 "This mapping info is not reinitialized for a cell vector!"));
1841 return compressed_cell_index;
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)
1861 template <
int dim,
int spacedim,
typename Number>
1866 if (do_cell_index_compression)
1870 ExcMessage(
"Mapping info object was not initialized for this"
1871 " active cell index!"));
1872 return cell_index_to_compressed_cell_index[
cell_index];
1879 template <
int dim,
int spacedim,
typename Number>
1882 const unsigned int unit_points_index_offset,
1883 const unsigned int n_q_points,
1884 const unsigned int n_q_points_unvectorized,
1887 const unsigned int n_lanes =
1890 for (
unsigned int q = 0; q < n_q_points; ++q)
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;
1896 for (
unsigned int d = 0; d < dim; ++d)
1898 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1904 template <
int dim,
int spacedim,
typename Number>
1907 const unsigned int unit_points_index_offset,
1908 const unsigned int n_q_points,
1909 const unsigned int n_q_points_unvectorized,
1912 const unsigned int n_lanes =
1915 for (
unsigned int q = 0; q < n_q_points; ++q)
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;
1921 for (
unsigned int d = 0; d < dim - 1; ++d)
1923 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1929 template <
int dim,
int spacedim,
typename Number>
1932 const unsigned int unit_points_index_offset,
1933 const unsigned int n_q_points,
1934 const unsigned int n_q_points_unvectorized,
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)
1941 const unsigned int n_lanes =
1944 for (
unsigned int q = 0; q < n_q_points; ++q)
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;
1953 if (q == 0 || !affine_cell)
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],
1976 if (additional_data.use_global_weights)
1979 JxW_values[offset], v) = weights[q * n_lanes + v];
1984 JxW_values[offset], v) =
1989 for (
unsigned int s = 0; s < spacedim; ++s)
1991 normal_vectors[offset][s], v) =
1993 if (update_flags_mapping &
1995 for (
unsigned int s = 0; s < spacedim; ++s)
1997 real_points[offset][s], v) =
2006 template <
int dim,
int spacedim,
typename Number>
2009 const unsigned int n_unit_point_batches)
2011 unit_points.resize(n_unit_point_batches);
2016 template <
int dim,
int spacedim,
typename Number>
2019 const unsigned int n_unit_point_batches)
2021 unit_points_faces.resize(n_unit_point_batches);
2026 template <
int dim,
int spacedim,
typename Number>
2029 const unsigned int n_data_point_batches,
2030 const bool is_face_centric)
2034 jacobians[0].resize(n_data_point_batches);
2035 if (is_face_centric)
2036 jacobians[1].resize(n_data_point_batches);
2040 inverse_jacobians[0].resize(n_data_point_batches);
2041 if (is_face_centric)
2042 inverse_jacobians[1].resize(n_data_point_batches);
2045 JxW_values.resize(n_data_point_batches);
2047 normal_vectors.resize(n_data_point_batches);
2049 real_points.resize(n_data_point_batches);
2054 template <
int dim,
int spacedim,
typename Number>
2059 const unsigned int offset)
const
2061 return unit_points.data() + offset;
2066 template <
int dim,
int spacedim,
typename Number>
2071 const unsigned int offset)
const
2073 return unit_points_faces.data() + offset;
2078 template <
int dim,
int spacedim,
typename Number>
2081 const unsigned int offset)
const
2083 return real_points.data() + offset;
2088 template <
int dim,
int spacedim,
typename Number>
2091 const unsigned int geometry_index)
const
2093 return unit_points_index[geometry_index];
2098 template <
int dim,
int spacedim,
typename Number>
2101 const unsigned int geometry_index)
const
2103 return data_index_offsets[geometry_index];
2107 template <
int dim,
int spacedim,
typename Number>
2110 const unsigned int geometry_index)
const
2112 return compressed_data_index_offsets[geometry_index];
2117 template <
int dim,
int spacedim,
typename Number>
2120 const bool is_interior)
const
2122 return jacobians[is_interior ? 0 : 1].data() + offset;
2127 template <
int dim,
int spacedim,
typename Number>
2130 const unsigned int offset,
2131 const bool is_interior)
const
2133 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2138 template <
int dim,
int spacedim,
typename Number>
2141 const unsigned int offset)
const
2143 return normal_vectors.data() + offset;
2148 template <
int dim,
int spacedim,
typename Number>
2149 inline const Number *
2152 return JxW_values.data() + offset;
2157 template <
int dim,
int spacedim,
typename Number>
2166 template <
int dim,
int spacedim,
typename Number>
2170 return update_flags;
2175 template <
int dim,
int spacedim,
typename Number>
2179 return update_flags_mapping;
2184 template <
int dim,
int spacedim,
typename Number>
2185 boost::signals2::connection
2187 const std::function<
void()> &set_is_reinitialized)
2189 return is_reinitialized.connect(set_is_reinitialized);
2194 template <
int dim,
int spacedim,
typename Number>
2201 memory += cell_type.capacity() *
2214 cell_index_to_compressed_cell_index);
2216 memory +=
sizeof(*this);
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Abstract base class for mapping classes.
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
@ faces_on_cells_in_vector
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
bool do_cell_index_compression
std::vector< unsigned int > data_index_offsets
bool is_face_state() const
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)
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
static constexpr unsigned char default_combined_face_orientation()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
@ 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
::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)