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,
73 internal_mapping_data->reinit(update_flags_mapping, quadrature);
75 cell_similarity = mapping->fill_fe_values(cell,
78 *internal_mapping_data,
91 &internal_mapping_data,
94 mapping_data.
initialize(quadrature.size(), update_flags_mapping);
96 internal_mapping_data->reinit(update_flags_mapping, quadrature);
98 mapping->fill_fe_immersed_surface_values(cell,
100 *internal_mapping_data,
111 const unsigned int face_no,
114 &internal_mapping_data,
125 *internal_mapping_data);
128 cell->reference_cell(),
131 cell->combined_face_orientation(face_no)),
134 mapping_q->fill_mapping_data_for_face_quadrature(
135 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
139 auto internal_mapping_data =
140 mapping->get_face_data(update_flags_mapping,
143 mapping->fill_fe_face_values(cell,
146 *internal_mapping_data,
152 template <
int dim,
int spacedim = dim>
155 const double diameter,
159 const auto jac_0 = inverse_jacobians[0];
160 const double zero_tolerance_double =
161 1e4 / diameter * std::numeric_limits<double>::epsilon() * 1024.;
162 bool jacobian_constant =
true;
163 for (
unsigned int q = 1; q < inverse_jacobians.size(); ++q)
166 for (
unsigned int d = 0; d < dim; ++d)
167 for (
unsigned int e = 0; e < spacedim; ++e)
168 if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
169 jacobian_constant =
false;
170 if (!jacobian_constant)
176 bool cell_cartesian = jacobian_constant && dim == spacedim;
177 for (
unsigned int d = 0; d < dim; ++d)
178 for (
unsigned int e = 0; e < dim; ++e)
180 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
182 cell_cartesian =
false;
188 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
189 else if (jacobian_constant)
190 return ::internal::MatrixFreeFunctions::GeometryType::affine;
192 return ::internal::MatrixFreeFunctions::GeometryType::general;
218 template <
int dim,
int spacedim = dim,
typename Number =
double>
226 Number>::vectorized_value_type;
322 template <
typename ContainerType>
325 const ContainerType &cell_iterator_range,
326 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
335 template <
typename ContainerType>
338 const ContainerType &cell_iterator_range,
346 template <
typename ContainerType>
349 const ContainerType &cell_iterator_range,
357 template <
typename ContainerType>
360 const ContainerType &cell_iterator_range,
368 template <
typename CellIteratorType>
370 reinit_faces(
const std::vector<std::pair<CellIteratorType, unsigned int>>
371 &face_iterator_range_interior,
407 const bool is_interior =
true)
const;
415 const bool is_interior =
true)
const;
459 template <
bool is_face>
481 const unsigned int geometry_index)
const;
517 template <
typename NumberType>
544 const bool is_face_centric =
false);
551 const unsigned int n_q_points,
560 const unsigned int n_q_points,
569 const unsigned int n_q_points,
572 const std::vector<double> &weights,
573 const unsigned int compressed_unit_point_index_offset,
574 const bool affine_cell,
575 const bool is_interior =
true);
586 template <
typename ContainerType,
typename QuadratureType>
589 const ContainerType &cell_iterator_range,
590 const std::vector<QuadratureType> &quadrature_vector,
591 const unsigned int n_unfiltered_cells,
637 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
678 std::vector<::internal::MatrixFreeFunctions::GeometryType>
cell_type;
711 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
721 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
773 template <
int dim,
int spacedim,
typename Number>
776 const unsigned int offset,
777 const bool is_interior)
const
779 const auto &face_pair = face_number[offset];
780 return is_interior ? face_pair.first : face_pair.second;
786 template <
int dim,
int spacedim,
typename Number>
792 , update_flags(update_flags)
794 , additional_data(additional_data)
821 template <
int dim,
int spacedim,
typename Number>
825 n_q_points_unvectorized.
clear();
826 unit_points_index.clear();
827 data_index_offsets.clear();
828 compressed_data_index_offsets.clear();
834 template <
int dim,
int spacedim,
typename Number>
838 const std::vector<
Point<dim>> &unit_points_in)
845 template <
int dim,
int spacedim,
typename Number>
851 quadrature.initialize(unit_points_in);
852 reinit(cell, quadrature);
857 template <
int dim,
int spacedim,
typename Number>
863 n_q_points_unvectorized.resize(1);
864 n_q_points_unvectorized[0] = quadrature.
size();
866 const unsigned int n_q_points =
867 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
869 const unsigned int n_q_points_data =
870 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
873 resize_unit_points(n_q_points);
874 resize_data_fields(n_q_points_data);
879 n_q_points_unvectorized[0],
886 update_flags_mapping,
890 internal_mapping_data,
894 if (!quadrature.
empty() &&
899 mapping_data.inverse_jacobians));
909 n_q_points_unvectorized[0],
916 unit_points_index.push_back(0);
917 data_index_offsets.push_back(0);
918 compressed_data_index_offsets.push_back(0);
920 state = State::single_cell;
925 template <
int dim,
int spacedim,
typename Number>
926 template <
typename ContainerType>
929 const ContainerType &cell_iterator_range,
930 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
931 const unsigned int n_unfiltered_cells)
933 const unsigned int n_cells = unit_points_vector.size();
935 std::distance(cell_iterator_range.begin(),
936 cell_iterator_range.end()));
938 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
943 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
948 template <
int dim,
int spacedim,
typename Number>
949 template <
typename ContainerType,
typename QuadratureType>
952 const ContainerType &cell_iterator_range,
953 const std::vector<QuadratureType> &quadrature_vector,
954 const unsigned int n_unfiltered_cells,
957 const QuadratureType &quadrature,
962 do_cell_index_compression =
965 const unsigned int n_cells = quadrature_vector.size();
967 std::distance(cell_iterator_range.begin(),
968 cell_iterator_range.end()));
970 n_q_points_unvectorized.reserve(n_cells);
972 cell_type.reserve(n_cells);
974 if (additional_data.store_cells)
975 cell_level_and_indices.resize(n_cells);
978 unit_points_index.reserve(n_cells + 1);
979 unit_points_index.push_back(0);
980 data_index_offsets.reserve(n_cells + 1);
981 data_index_offsets.push_back(0);
982 for (
const auto &quadrature : quadrature_vector)
984 const unsigned int n_points = quadrature.size();
985 n_q_points_unvectorized.push_back(n_points);
987 const unsigned int n_q_points =
988 compute_n_q_points<VectorizedArrayType>(n_points);
989 unit_points_index.push_back(unit_points_index.back() + n_q_points);
991 const unsigned int n_q_points_data =
992 compute_n_q_points<Number>(n_points);
993 data_index_offsets.push_back(data_index_offsets.back() +
997 const unsigned int n_unit_points = unit_points_index.back();
998 const unsigned int n_data_points = data_index_offsets.back();
1001 resize_unit_points(n_unit_points);
1002 resize_data_fields(n_data_points);
1004 if (do_cell_index_compression)
1005 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1009 unsigned int size_compressed_data = 0;
1011 for (
const auto &cell : cell_iterator_range)
1013 if (additional_data.store_cells)
1016 cell_level_and_indices[
cell_index] = {cell->level(), cell->index()};
1019 const auto &quadrature = quadrature_vector[
cell_index];
1020 const bool empty = quadrature.empty();
1023 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1025 store_unit_points(unit_points_index[
cell_index],
1028 quadrature.get_points());
1031 compute_mapping_data(cell, quadrature, mapping_data);
1034 const unsigned int n_q_points_data =
1035 compute_n_q_points<Number>(n_q_points_unvectorized[
cell_index]);
1041 cell_type.push_back(
1046 cell_type.push_back(
1052 const bool affine_cells =
1061 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1066 const auto comparison_result =
1077 comparison_result ==
1080 compressed_data_index_offsets.push_back(
1081 compressed_data_index_offsets.back());
1085 const unsigned int n_compressed_data_last_cell =
1089 compute_n_q_points<Number>(
1092 compressed_data_index_offsets.push_back(
1093 compressed_data_index_offsets.back() +
1094 n_compressed_data_last_cell);
1098 compressed_data_index_offsets.push_back(0);
1101 mapping_data_previous_cell = mapping_data;
1103 store_mapping_data(data_index_offsets[
cell_index],
1107 quadrature.get_weights(),
1116 size_compressed_data = compressed_data_index_offsets.back() + 1;
1118 size_compressed_data =
1120 compressed_data_index_offsets.back() + n_q_points_data);
1122 if (do_cell_index_compression)
1123 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1131 jacobians[0].resize(size_compressed_data);
1132 jacobians[0].shrink_to_fit();
1136 inverse_jacobians[0].resize(size_compressed_data);
1137 inverse_jacobians[0].shrink_to_fit();
1140 state = State::cell_vector;
1145 template <
int dim,
int spacedim,
typename Number>
1146 template <
typename ContainerType>
1149 const ContainerType &cell_iterator_range,
1151 const unsigned int n_unfiltered_cells)
1153 auto compute_mapping_data_for_cells =
1161 update_flags_mapping,
1165 internal_mapping_data,
1169 do_reinit_cells<ContainerType, Quadrature<dim>>(
1170 cell_iterator_range,
1173 compute_mapping_data_for_cells);
1178 template <
int dim,
int spacedim,
typename Number>
1179 template <
typename ContainerType>
1182 const ContainerType &cell_iterator_range,
1184 const unsigned int n_unfiltered_cells)
1187 additional_data.use_global_weights ==
false,
1189 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1196 auto compute_mapping_data_for_surface =
1203 update_flags_mapping,
1206 internal_mapping_data,
1210 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1211 cell_iterator_range,
1214 compute_mapping_data_for_surface);
1219 template <
int dim,
int spacedim,
typename Number>
1220 template <
typename ContainerType>
1223 const ContainerType &cell_iterator_range,
1225 const unsigned int n_unfiltered_cells)
1231 do_cell_index_compression =
1234 const unsigned int n_cells = quadrature_vector.size();
1236 std::distance(cell_iterator_range.begin(),
1237 cell_iterator_range.end()));
1240 cell_index_offset.resize(n_cells);
1241 unsigned int n_faces = 0;
1243 for (
const auto &cell : cell_iterator_range)
1246 n_faces += cell->n_faces();
1250 n_q_points_unvectorized.reserve(n_faces);
1252 cell_type.reserve(n_faces);
1255 unit_points_index.resize(n_faces + 1);
1256 data_index_offsets.resize(n_faces + 1);
1258 unsigned int n_unit_points = 0;
1259 unsigned int n_data_points = 0;
1260 for (
const auto &cell : cell_iterator_range)
1262 for (
const auto &f : cell->face_indices())
1264 const unsigned int current_face_index =
1267 unit_points_index[current_face_index] = n_unit_points;
1268 data_index_offsets[current_face_index] = n_data_points;
1270 const unsigned int n_points =
1272 n_q_points_unvectorized.push_back(n_points);
1274 const unsigned int n_q_points =
1275 compute_n_q_points<VectorizedArrayType>(n_points);
1276 n_unit_points += n_q_points;
1278 const unsigned int n_q_points_data =
1279 compute_n_q_points<Number>(n_points);
1280 n_data_points += n_q_points_data;
1285 unit_points_index[n_faces] = n_unit_points;
1286 data_index_offsets[n_faces] = n_data_points;
1289 if (do_cell_index_compression)
1290 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1295 resize_unit_points_faces(n_unit_points);
1296 resize_data_fields(n_data_points);
1300 bool first_set =
false;
1301 unsigned int size_compressed_data = 0;
1303 for (
const auto &cell : cell_iterator_range)
1305 const auto &quadratures_on_faces = quadrature_vector[
cell_index];
1307 Assert(quadratures_on_faces.size() == cell->n_faces(),
1311 for (
const auto &f : cell->face_indices())
1313 const auto &quadrature_on_face = quadratures_on_faces[f];
1314 const bool empty = quadrature_on_face.empty();
1316 const unsigned int current_face_index =
1320 const unsigned int n_q_points =
1321 compute_n_q_points<VectorizedArrayType>(
1322 n_q_points_unvectorized[current_face_index]);
1323 store_unit_points_faces(unit_points_index[current_face_index],
1325 n_q_points_unvectorized[current_face_index],
1326 quadrature_on_face.get_points());
1330 update_flags_mapping,
1334 internal_mapping_data,
1342 cell->diameter(), mapping_data.inverse_jacobians));
1346 mapping_data_first = mapping_data;
1351 cell_type.push_back(
1354 if (current_face_index > 0)
1357 const bool affine_cells =
1358 cell_type[current_face_index] <=
1360 cell_type[current_face_index - 1] <=
1366 1e4 / cell->diameter() *
1367 std::numeric_limits<double>::epsilon() * 1024.);
1371 const auto comparison_result =
1372 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1376 mapping_data.inverse_jacobians[0],
1382 comparison_result ==
1385 compressed_data_index_offsets.push_back(
1386 compressed_data_index_offsets.back());
1388 else if (first_set &&
1389 (cell_type[current_face_index] <=
1392 mapping_data.inverse_jacobians[0],
1395 double>::ComparisonResult::equal))
1397 compressed_data_index_offsets.push_back(0);
1401 const unsigned int n_compressed_data_last_cell =
1402 cell_type[current_face_index - 1] <=
1405 compute_n_q_points<Number>(
1406 n_q_points_unvectorized[current_face_index - 1]);
1408 compressed_data_index_offsets.push_back(
1409 compressed_data_index_offsets.back() +
1410 n_compressed_data_last_cell);
1414 compressed_data_index_offsets.push_back(0);
1417 mapping_data_previous_cell = mapping_data;
1419 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1420 n_q_points_unvectorized[current_face_index]);
1421 store_mapping_data(data_index_offsets[current_face_index],
1423 n_q_points_unvectorized[current_face_index],
1425 quadrature_on_face.get_weights(),
1426 data_index_offsets[current_face_index],
1427 cell_type[current_face_index] <=
1432 if (cell_type[current_face_index] <=
1434 size_compressed_data = compressed_data_index_offsets.back() + 1;
1436 size_compressed_data =
1438 compressed_data_index_offsets.back() +
1441 if (do_cell_index_compression)
1442 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1450 jacobians[0].resize(size_compressed_data);
1451 jacobians[0].shrink_to_fit();
1455 inverse_jacobians[0].resize(size_compressed_data);
1456 inverse_jacobians[0].shrink_to_fit();
1459 state = State::faces_on_cells_in_vector;
1464 template <
int dim,
int spacedim,
typename Number>
1465 template <
typename CellIteratorType>
1468 const std::vector<std::pair<CellIteratorType, unsigned int>>
1469 &face_iterator_range_interior,
1474 do_cell_index_compression =
false;
1479 const unsigned int n_faces = quadrature_vector.size();
1481 std::distance(face_iterator_range_interior.begin(),
1482 face_iterator_range_interior.end()));
1484 n_q_points_unvectorized.reserve(n_faces);
1486 cell_type.reserve(n_faces);
1487 face_number.reserve(n_faces);
1490 unit_points_index.reserve(n_faces + 1);
1491 unit_points_index.push_back(0);
1492 data_index_offsets.reserve(n_faces + 1);
1493 data_index_offsets.push_back(0);
1494 for (
const auto &quadrature : quadrature_vector)
1496 const unsigned int n_points = quadrature.size();
1497 n_q_points_unvectorized.push_back(n_points);
1499 const unsigned int n_q_points =
1500 compute_n_q_points<VectorizedArrayType>(n_points);
1501 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1503 const unsigned int n_q_points_data =
1504 compute_n_q_points<Number>(n_points);
1505 data_index_offsets.push_back(data_index_offsets.back() +
1509 const unsigned int n_unit_points = unit_points_index.back();
1510 const unsigned int n_data_points = data_index_offsets.back();
1513 resize_unit_points_faces(n_unit_points);
1514 resize_data_fields(n_data_points,
true);
1516 std::array<MappingData, 2> mapping_data;
1517 std::array<MappingData, 2> mapping_data_previous_cell;
1518 std::array<MappingData, 2> mapping_data_first;
1519 bool first_set =
false;
1520 unsigned int size_compressed_data = 0;
1521 unsigned int face_index = 0;
1522 for (
const auto &cell_and_f : face_iterator_range_interior)
1524 const auto &quadrature_on_face = quadrature_vector[face_index];
1525 const bool empty = quadrature_on_face.empty();
1528 const auto &cell_m = cell_and_f.first;
1529 const auto f_m = cell_and_f.second;
1532 const auto &cell_p =
1533 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1535 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1537 face_number.emplace_back(f_m, f_p);
1540 cell_m->combined_face_orientation(f_m) ==
1542 cell_p->combined_face_orientation(f_p) ==
1545 "Non standard face orientation is currently not implemented."));
1548 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1549 n_q_points_unvectorized[face_index]);
1550 store_unit_points_faces(unit_points_index[face_index],
1552 n_q_points_unvectorized[face_index],
1553 quadrature_on_face.get_points());
1558 update_flags_mapping,
1562 internal_mapping_data,
1568 update_flags_mapping,
1572 internal_mapping_data,
1582 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1584 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1590 mapping_data_first = mapping_data;
1595 cell_type.push_back(
1601 const bool affine_cells =
1602 cell_type[face_index] <=
1604 cell_type[face_index - 1] <=
1610 1e4 / cell_m->diameter() *
1611 std::numeric_limits<double>::epsilon() * 1024.);
1615 const auto comparison_result_m =
1616 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1617 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1620 mapping_data[0].inverse_jacobians[0],
1621 mapping_data_previous_cell[0].inverse_jacobians[0]);
1623 const auto comparison_result_p =
1624 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1625 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1628 mapping_data[1].inverse_jacobians[0],
1629 mapping_data_previous_cell[1].inverse_jacobians[0]);
1634 comparison_result_m ==
1636 comparison_result_p ==
1639 compressed_data_index_offsets.push_back(
1640 compressed_data_index_offsets.back());
1642 else if (first_set &&
1643 (cell_type[face_index] <=
1646 mapping_data[0].inverse_jacobians[0],
1647 mapping_data_first[0].inverse_jacobians[0]) ==
1649 double>::ComparisonResult::equal) &&
1651 mapping_data[1].inverse_jacobians[0],
1652 mapping_data_first[1].inverse_jacobians[0]) ==
1655 compressed_data_index_offsets.push_back(0);
1659 const unsigned int n_compressed_data_last_cell =
1660 cell_type[face_index - 1] <=
1663 compute_n_q_points<Number>(
1664 n_q_points_unvectorized[face_index - 1]);
1666 compressed_data_index_offsets.push_back(
1667 compressed_data_index_offsets.back() +
1668 n_compressed_data_last_cell);
1672 compressed_data_index_offsets.push_back(0);
1675 mapping_data_previous_cell = mapping_data;
1677 const unsigned int n_q_points_data =
1678 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1681 store_mapping_data(data_index_offsets[face_index],
1683 n_q_points_unvectorized[face_index],
1685 quadrature_on_face.get_weights(),
1686 data_index_offsets[face_index],
1687 cell_type[face_index] <=
1693 store_mapping_data(data_index_offsets[face_index],
1695 n_q_points_unvectorized[face_index],
1697 quadrature_on_face.get_weights(),
1698 data_index_offsets[face_index],
1699 cell_type[face_index] <=
1705 if (cell_type[face_index] <=
1707 size_compressed_data = compressed_data_index_offsets.back() + 1;
1709 size_compressed_data =
1711 compressed_data_index_offsets.back() + n_q_points_data);
1718 jacobians[0].resize(size_compressed_data);
1719 jacobians[0].shrink_to_fit();
1720 jacobians[1].resize(size_compressed_data);
1721 jacobians[1].shrink_to_fit();
1725 inverse_jacobians[0].resize(size_compressed_data);
1726 inverse_jacobians[0].shrink_to_fit();
1727 inverse_jacobians[1].resize(size_compressed_data);
1728 inverse_jacobians[1].shrink_to_fit();
1731 state = State::face_vector;
1736 template <
int dim,
int spacedim,
typename Number>
1740 return state == State::faces_on_cells_in_vector;
1745 template <
int dim,
int spacedim,
typename Number>
1748 const unsigned int geometry_index)
const
1750 return n_q_points_unvectorized[geometry_index];
1754 template <
int dim,
int spacedim,
typename Number>
1757 const unsigned int geometry_index)
const
1759 return cell_type[geometry_index];
1764 template <
int dim,
int spacedim,
typename Number>
1770 additional_data.store_cells,
1772 "Cells have been not stored. You can enable this by Additional::store_cells."));
1780 template <
int dim,
int spacedim,
typename Number>
1781 template <
typename NumberType>
1784 const unsigned int n_q_points_unvectorized)
1786 const unsigned int n_lanes =
1788 const unsigned int n_filled_lanes_last_batch =
1789 n_q_points_unvectorized % n_lanes;
1790 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1791 if (n_filled_lanes_last_batch > 0)
1798 template <
int dim,
int spacedim,
typename Number>
1799 template <
bool is_face>
1803 const unsigned int face_number)
const
1808 const unsigned int compressed_cell_index =
1812 Assert(state == State::cell_vector,
1814 "This mapping info is not reinitialized for a cell vector!"));
1815 return compressed_cell_index;
1821 "cell_index has to be set if face_number is specified!"));
1822 Assert(state == State::faces_on_cells_in_vector ||
1823 state == State::face_vector,
1824 ExcMessage(
"This mapping info is not reinitialized for faces"
1825 " on cells in a vector!"));
1826 if (state == State::faces_on_cells_in_vector)
1827 return cell_index_offset[compressed_cell_index] + face_number;
1828 else if (state == State::face_vector)
1835 template <
int dim,
int spacedim,
typename Number>
1840 if (do_cell_index_compression)
1844 ExcMessage(
"Mapping info object was not initialized for this"
1845 " active cell index!"));
1846 return cell_index_to_compressed_cell_index[
cell_index];
1853 template <
int dim,
int spacedim,
typename Number>
1856 const unsigned int unit_points_index_offset,
1857 const unsigned int n_q_points,
1858 const unsigned int n_q_points_unvectorized,
1861 const unsigned int n_lanes =
1864 for (
unsigned int q = 0; q < n_q_points; ++q)
1866 const unsigned int offset = unit_points_index_offset + q;
1867 for (
unsigned int v = 0;
1868 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1870 for (
unsigned int d = 0; d < dim; ++d)
1872 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1878 template <
int dim,
int spacedim,
typename Number>
1881 const unsigned int unit_points_index_offset,
1882 const unsigned int n_q_points,
1883 const unsigned int n_q_points_unvectorized,
1886 const unsigned int n_lanes =
1889 for (
unsigned int q = 0; q < n_q_points; ++q)
1891 const unsigned int offset = unit_points_index_offset + q;
1892 for (
unsigned int v = 0;
1893 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1895 for (
unsigned int d = 0; d < dim - 1; ++d)
1897 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1903 template <
int dim,
int spacedim,
typename Number>
1906 const unsigned int unit_points_index_offset,
1907 const unsigned int n_q_points,
1908 const unsigned int n_q_points_unvectorized,
1910 const std::vector<double> &weights,
1911 const unsigned int compressed_unit_point_index_offset,
1912 const bool affine_cell,
1913 const bool is_interior)
1915 const unsigned int n_lanes =
1918 for (
unsigned int q = 0; q < n_q_points; ++q)
1920 const unsigned int offset = unit_points_index_offset + q;
1921 const unsigned int compressed_offset =
1922 compressed_unit_point_index_offset + q;
1923 for (
unsigned int v = 0;
1924 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1927 if (q == 0 || !affine_cell)
1930 for (
unsigned int d = 0; d < dim; ++d)
1931 for (
unsigned int s = 0; s < spacedim; ++s)
1933 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1934 v) = mapping_data.
jacobians[q * n_lanes + v][d][s];
1935 if (update_flags_mapping &
1937 for (
unsigned int d = 0; d < dim; ++d)
1938 for (
unsigned int s = 0; s < spacedim; ++s)
1940 inverse_jacobians[is_interior ? 0 : 1]
1941 [compressed_offset][d][s],
1950 if (additional_data.use_global_weights)
1953 JxW_values[offset], v) = weights[q * n_lanes + v];
1958 JxW_values[offset], v) =
1963 for (
unsigned int s = 0; s < spacedim; ++s)
1965 normal_vectors[offset][s], v) =
1967 if (update_flags_mapping &
1969 for (
unsigned int s = 0; s < spacedim; ++s)
1971 real_points[offset][s], v) =
1980 template <
int dim,
int spacedim,
typename Number>
1983 const unsigned int n_unit_point_batches)
1985 unit_points.resize(n_unit_point_batches);
1990 template <
int dim,
int spacedim,
typename Number>
1993 const unsigned int n_unit_point_batches)
1995 unit_points_faces.resize(n_unit_point_batches);
2000 template <
int dim,
int spacedim,
typename Number>
2003 const unsigned int n_data_point_batches,
2004 const bool is_face_centric)
2008 jacobians[0].resize(n_data_point_batches);
2009 if (is_face_centric)
2010 jacobians[1].resize(n_data_point_batches);
2014 inverse_jacobians[0].resize(n_data_point_batches);
2015 if (is_face_centric)
2016 inverse_jacobians[1].resize(n_data_point_batches);
2019 JxW_values.resize(n_data_point_batches);
2021 normal_vectors.resize(n_data_point_batches);
2023 real_points.resize(n_data_point_batches);
2028 template <
int dim,
int spacedim,
typename Number>
2033 const unsigned int offset)
const
2035 return unit_points.data() + offset;
2040 template <
int dim,
int spacedim,
typename Number>
2045 const unsigned int offset)
const
2047 return unit_points_faces.data() + offset;
2052 template <
int dim,
int spacedim,
typename Number>
2055 const unsigned int offset)
const
2057 return real_points.data() + offset;
2062 template <
int dim,
int spacedim,
typename Number>
2065 const unsigned int geometry_index)
const
2067 return unit_points_index[geometry_index];
2072 template <
int dim,
int spacedim,
typename Number>
2075 const unsigned int geometry_index)
const
2077 return data_index_offsets[geometry_index];
2081 template <
int dim,
int spacedim,
typename Number>
2084 const unsigned int geometry_index)
const
2086 return compressed_data_index_offsets[geometry_index];
2091 template <
int dim,
int spacedim,
typename Number>
2094 const bool is_interior)
const
2096 return jacobians[is_interior ? 0 : 1].data() + offset;
2101 template <
int dim,
int spacedim,
typename Number>
2104 const unsigned int offset,
2105 const bool is_interior)
const
2107 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2112 template <
int dim,
int spacedim,
typename Number>
2115 const unsigned int offset)
const
2117 return normal_vectors.data() + offset;
2122 template <
int dim,
int spacedim,
typename Number>
2123 inline const Number *
2126 return JxW_values.data() + offset;
2131 template <
int dim,
int spacedim,
typename Number>
2140 template <
int dim,
int spacedim,
typename Number>
2144 return update_flags;
2149 template <
int dim,
int spacedim,
typename Number>
2153 return update_flags_mapping;
2158 template <
int dim,
int spacedim,
typename Number>
2165 memory += cell_type.capacity() *
2178 cell_index_to_compressed_cell_index);
2180 memory +=
sizeof(*this);
370 reinit_faces(
const std::vector<std::pair<CellIteratorType, unsigned int>> {
…}
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
void initialize_face(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
unsigned int get_face_number(const unsigned int offset, const bool is_interior) const
const Mapping< dim, spacedim > & get_mapping() const
void reinit_faces(const ContainerType &cell_iterator_range, const std::vector< std::vector< Quadrature< dim - 1 > > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
std::array< AlignedVector< DerivativeForm< 1, spacedim, dim, Number > >, 2 > inverse_jacobians
const Point< dim, VectorizedArrayType > * get_unit_point(const unsigned int offset) const
unsigned int get_n_q_points_unvectorized(const unsigned int geometry_index) const
const ObserverPointer< const Mapping< dim, spacedim > > mapping
unsigned int compute_compressed_cell_index(const unsigned int cell_index) const
const Point< dim - 1, VectorizedArrayType > * get_unit_point_faces(const unsigned int offset) const
Quadrature< dim > quadrature
std::vector< unsigned int > cell_index_offset
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Quadrature< dim > &quadrature)
void store_mapping_data(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const MappingData &mapping_data, const std::vector< double > &weights, const unsigned int compressed_unit_point_index_offset, const bool affine_cell, const bool is_interior=true)
AlignedVector< Tensor< 1, spacedim, Number > > normal_vectors
unsigned int compute_geometry_index_offset(const unsigned int cell_index, const unsigned int face_number) const
std::vector<::internal::MatrixFreeFunctions::GeometryType > cell_type
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim > > &unit_points)
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< Quadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
MappingInfo & operator=(const MappingInfo &)=delete
ObserverPointer< const Triangulation< dim, spacedim > > triangulation
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > internal_mapping_data
void reinit_surface(const ContainerType &cell_iterator_range, const std::vector< ImmersedSurfaceQuadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_unit_points_faces(const unsigned int n_unit_point_batches)
const Tensor< 1, spacedim, Number > * get_normal_vector(const unsigned int offset) const
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const std::vector< Point< dim > > &unit_points)
const DerivativeForm< 1, dim, spacedim, Number > * get_jacobian(const unsigned int offset, const bool is_interior=true) const
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< std::vector< Point< dim > > > &unit_points_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_data_fields(const unsigned int n_data_point_batches, const bool is_face_centric=false)
void resize_unit_points(const unsigned int n_unit_point_batches)
std::size_t memory_consumption() const
UpdateFlags get_update_flags() const
void store_unit_points_faces(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const std::vector< Point< dim - 1 > > &points)
unsigned int compute_unit_point_index_offset(const unsigned int geometry_index) const
std::vector< unsigned int > cell_index_to_compressed_cell_index
AlignedVector< Point< spacedim, Number > > real_points
unsigned int compute_data_index_offset(const unsigned int geometry_index) const
@ 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)
MappingInfo(const MappingInfo &)=delete
void reinit_faces(const std::vector< std::pair< CellIteratorType, unsigned int > > &face_iterator_range_interior, const std::vector< Quadrature< dim - 1 > > &quadrature_vector)
UpdateFlags update_flags_mapping
void do_reinit_cells(const ContainerType &cell_iterator_range, const std::vector< QuadratureType > &quadrature_vector, const unsigned int n_unfiltered_cells, const std::function< void(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const QuadratureType &quadrature, MappingData &mapping_data)> &compute_mapping_data)
UpdateFlags get_update_flags_mapping() const
Triangulation< dim, spacedim >::cell_iterator get_cell_iterator(const unsigned int cell_index) const
static void compute_mapping_data_for_immersed_surface_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ImmersedSurfaceQuadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_face_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_quadrature(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, CellSimilarity::Similarity &cell_similarity, const Quadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static UpdateFlags required_update_flags(const ObserverPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#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::vector< index_type > data
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
::internal::MatrixFreeFunctions::GeometryType compute_geometry_type(const double diameter, const std::vector< DerivativeForm< 1, spacedim, dim, double > > &inverse_jacobians)
constexpr unsigned int invalid_unsigned_int
constexpr types::geometric_orientation default_geometric_orientation
::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)