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 ReferenceCells::get_hypercube<dim>(),
131 cell->face_orientation(face_no),
132 cell->face_flip(face_no),
133 cell->face_rotation(face_no)),
136 mapping_q->fill_mapping_data_for_face_quadrature(
137 cell, face_no, quadrature, *internal_mapping_data, mapping_data);
141 auto internal_mapping_data =
142 mapping->get_face_data(update_flags_mapping,
145 mapping->fill_fe_face_values(cell,
148 *internal_mapping_data,
154 template <
int dim,
int spacedim = dim>
157 const double diameter,
161 const auto jac_0 = inverse_jacobians[0];
162 const double zero_tolerance_double =
163 1e4 / diameter * std::numeric_limits<double>::epsilon() * 1024.;
164 bool jacobian_constant =
true;
165 for (
unsigned int q = 1; q < inverse_jacobians.size(); ++q)
168 for (
unsigned int d = 0; d < dim; ++d)
169 for (
unsigned int e = 0; e < spacedim; ++e)
170 if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
171 jacobian_constant =
false;
172 if (!jacobian_constant)
178 bool cell_cartesian = jacobian_constant && dim == spacedim;
179 for (
unsigned int d = 0; d < dim; ++d)
180 for (
unsigned int e = 0; e < dim; ++e)
182 if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
184 cell_cartesian =
false;
190 return ::internal::MatrixFreeFunctions::GeometryType::cartesian;
191 else if (jacobian_constant)
192 return ::internal::MatrixFreeFunctions::GeometryType::affine;
194 return ::internal::MatrixFreeFunctions::GeometryType::general;
220 template <
int dim,
int spacedim = dim,
typename Number =
double>
228 Number>::vectorized_value_type;
324 template <
typename ContainerType>
327 const ContainerType &cell_iterator_range,
328 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
337 template <
typename ContainerType>
340 const ContainerType &cell_iterator_range,
348 template <
typename ContainerType>
351 const ContainerType &cell_iterator_range,
359 template <
typename ContainerType>
362 const ContainerType &cell_iterator_range,
370 template <
typename CellIteratorType>
372 reinit_faces(
const std::vector<std::pair<CellIteratorType, unsigned int>>
373 &face_iterator_range_interior,
409 const bool is_interior =
true)
const;
417 const bool is_interior =
true)
const;
461 boost::signals2::connection
467 template <
bool is_face>
489 const unsigned int geometry_index)
const;
525 template <
typename NumberType>
552 const bool is_face_centric =
false);
559 const unsigned int n_q_points,
568 const unsigned int n_q_points,
577 const unsigned int n_q_points,
580 const std::vector<double> &weights,
581 const unsigned int compressed_unit_point_index_offset,
582 const bool affine_cell,
583 const bool is_interior =
true);
594 template <
typename ContainerType,
typename QuadratureType>
597 const ContainerType &cell_iterator_range,
598 const std::vector<QuadratureType> &quadrature_vector,
599 const unsigned int n_unfiltered_cells,
645 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
686 std::vector<::internal::MatrixFreeFunctions::GeometryType>
cell_type;
719 std::array<AlignedVector<DerivativeForm<1, dim, spacedim, Number>>, 2>
729 std::array<AlignedVector<DerivativeForm<1, spacedim, dim, Number>>, 2>
787 template <
int dim,
int spacedim,
typename Number>
790 const unsigned int offset,
791 const bool is_interior)
const
793 const auto &face_pair = face_number[offset];
794 return is_interior ? face_pair.first : face_pair.second;
800 template <
int dim,
int spacedim,
typename Number>
806 , update_flags(update_flags)
808 , additional_data(additional_data)
835 template <
int dim,
int spacedim,
typename Number>
839 n_q_points_unvectorized.
clear();
840 unit_points_index.clear();
841 data_index_offsets.clear();
842 compressed_data_index_offsets.clear();
848 template <
int dim,
int spacedim,
typename Number>
852 const std::vector<
Point<dim>> &unit_points_in)
859 template <
int dim,
int spacedim,
typename Number>
865 quadrature.initialize(unit_points_in);
866 reinit(cell, quadrature);
871 template <
int dim,
int spacedim,
typename Number>
877 n_q_points_unvectorized.resize(1);
878 n_q_points_unvectorized[0] = quadrature.
size();
880 const unsigned int n_q_points =
881 compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
883 const unsigned int n_q_points_data =
884 compute_n_q_points<Number>(n_q_points_unvectorized[0]);
887 resize_unit_points(n_q_points);
888 resize_data_fields(n_q_points_data);
893 n_q_points_unvectorized[0],
900 update_flags_mapping,
904 internal_mapping_data,
908 if (!quadrature.
empty() &&
913 mapping_data.inverse_jacobians));
923 n_q_points_unvectorized[0],
930 unit_points_index.push_back(0);
931 data_index_offsets.push_back(0);
932 compressed_data_index_offsets.push_back(0);
934 state = State::single_cell;
940 template <
int dim,
int spacedim,
typename Number>
941 template <
typename ContainerType>
944 const ContainerType &cell_iterator_range,
945 const std::vector<std::vector<
Point<dim>>> &unit_points_vector,
946 const unsigned int n_unfiltered_cells)
948 const unsigned int n_cells = unit_points_vector.size();
950 std::distance(cell_iterator_range.begin(),
951 cell_iterator_range.end()));
953 std::vector<Quadrature<dim>> quadrature_vector(n_cells);
958 reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
963 template <
int dim,
int spacedim,
typename Number>
964 template <
typename ContainerType,
typename QuadratureType>
967 const ContainerType &cell_iterator_range,
968 const std::vector<QuadratureType> &quadrature_vector,
969 const unsigned int n_unfiltered_cells,
972 const QuadratureType &quadrature,
977 do_cell_index_compression =
980 const unsigned int n_cells = quadrature_vector.size();
982 std::distance(cell_iterator_range.begin(),
983 cell_iterator_range.end()));
985 n_q_points_unvectorized.reserve(n_cells);
987 cell_type.reserve(n_cells);
989 if (additional_data.store_cells)
990 cell_level_and_indices.resize(n_cells);
993 unit_points_index.reserve(n_cells + 1);
994 unit_points_index.push_back(0);
995 data_index_offsets.reserve(n_cells + 1);
996 data_index_offsets.push_back(0);
997 for (
const auto &quadrature : quadrature_vector)
999 const unsigned int n_points = quadrature.size();
1000 n_q_points_unvectorized.push_back(n_points);
1002 const unsigned int n_q_points =
1003 compute_n_q_points<VectorizedArrayType>(n_points);
1004 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1006 const unsigned int n_q_points_data =
1007 compute_n_q_points<Number>(n_points);
1008 data_index_offsets.push_back(data_index_offsets.back() +
1012 const unsigned int n_unit_points = unit_points_index.back();
1013 const unsigned int n_data_points = data_index_offsets.back();
1016 resize_unit_points(n_unit_points);
1017 resize_data_fields(n_data_points);
1019 if (do_cell_index_compression)
1020 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1024 unsigned int size_compressed_data = 0;
1026 for (
const auto &cell : cell_iterator_range)
1028 if (additional_data.store_cells)
1031 cell_level_and_indices[
cell_index] = {cell->level(), cell->index()};
1034 const auto &quadrature = quadrature_vector[
cell_index];
1035 const bool empty = quadrature.empty();
1038 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1040 store_unit_points(unit_points_index[
cell_index],
1043 quadrature.get_points());
1046 compute_mapping_data(cell, quadrature, mapping_data);
1049 const unsigned int n_q_points_data =
1050 compute_n_q_points<Number>(n_q_points_unvectorized[
cell_index]);
1056 cell_type.push_back(
1061 cell_type.push_back(
1067 const bool affine_cells =
1076 1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
1081 const auto comparison_result =
1092 comparison_result ==
1095 compressed_data_index_offsets.push_back(
1096 compressed_data_index_offsets.back());
1100 const unsigned int n_compressed_data_last_cell =
1104 compute_n_q_points<Number>(
1107 compressed_data_index_offsets.push_back(
1108 compressed_data_index_offsets.back() +
1109 n_compressed_data_last_cell);
1113 compressed_data_index_offsets.push_back(0);
1116 mapping_data_previous_cell = mapping_data;
1118 store_mapping_data(data_index_offsets[
cell_index],
1122 quadrature.get_weights(),
1131 size_compressed_data = compressed_data_index_offsets.back() + 1;
1133 size_compressed_data =
1135 compressed_data_index_offsets.back() + n_q_points_data);
1137 if (do_cell_index_compression)
1138 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1146 jacobians[0].resize(size_compressed_data);
1147 jacobians[0].shrink_to_fit();
1151 inverse_jacobians[0].resize(size_compressed_data);
1152 inverse_jacobians[0].shrink_to_fit();
1155 state = State::cell_vector;
1161 template <
int dim,
int spacedim,
typename Number>
1162 template <
typename ContainerType>
1165 const ContainerType &cell_iterator_range,
1167 const unsigned int n_unfiltered_cells)
1169 auto compute_mapping_data_for_cells =
1177 update_flags_mapping,
1181 internal_mapping_data,
1185 do_reinit_cells<ContainerType, Quadrature<dim>>(
1186 cell_iterator_range,
1189 compute_mapping_data_for_cells);
1194 template <
int dim,
int spacedim,
typename Number>
1195 template <
typename ContainerType>
1198 const ContainerType &cell_iterator_range,
1200 const unsigned int n_unfiltered_cells)
1203 additional_data.use_global_weights ==
false,
1205 "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
1212 auto compute_mapping_data_for_surface =
1219 update_flags_mapping,
1222 internal_mapping_data,
1226 do_reinit_cells<ContainerType, ImmersedSurfaceQuadrature<dim>>(
1227 cell_iterator_range,
1230 compute_mapping_data_for_surface);
1235 template <
int dim,
int spacedim,
typename Number>
1236 template <
typename ContainerType>
1239 const ContainerType &cell_iterator_range,
1241 const unsigned int n_unfiltered_cells)
1247 do_cell_index_compression =
1250 const unsigned int n_cells = quadrature_vector.size();
1252 std::distance(cell_iterator_range.begin(),
1253 cell_iterator_range.end()));
1256 cell_index_offset.resize(n_cells);
1257 unsigned int n_faces = 0;
1259 for (
const auto &cell : cell_iterator_range)
1262 n_faces += cell->n_faces();
1266 n_q_points_unvectorized.reserve(n_faces);
1268 cell_type.reserve(n_faces);
1271 unit_points_index.resize(n_faces + 1);
1272 data_index_offsets.resize(n_faces + 1);
1274 unsigned int n_unit_points = 0;
1275 unsigned int n_data_points = 0;
1276 for (
const auto &cell : cell_iterator_range)
1278 for (
const auto &f : cell->face_indices())
1280 const unsigned int current_face_index =
1283 unit_points_index[current_face_index] = n_unit_points;
1284 data_index_offsets[current_face_index] = n_data_points;
1286 const unsigned int n_points =
1288 n_q_points_unvectorized.push_back(n_points);
1290 const unsigned int n_q_points =
1291 compute_n_q_points<VectorizedArrayType>(n_points);
1292 n_unit_points += n_q_points;
1294 const unsigned int n_q_points_data =
1295 compute_n_q_points<Number>(n_points);
1296 n_data_points += n_q_points_data;
1301 unit_points_index[n_faces] = n_unit_points;
1302 data_index_offsets[n_faces] = n_data_points;
1305 if (do_cell_index_compression)
1306 cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
1311 resize_unit_points_faces(n_unit_points);
1312 resize_data_fields(n_data_points);
1316 bool first_set =
false;
1317 unsigned int size_compressed_data = 0;
1319 for (
const auto &cell : cell_iterator_range)
1321 const auto &quadratures_on_faces = quadrature_vector[
cell_index];
1323 Assert(quadratures_on_faces.size() == cell->n_faces(),
1327 for (
const auto &f : cell->face_indices())
1329 const auto &quadrature_on_face = quadratures_on_faces[f];
1330 const bool empty = quadrature_on_face.empty();
1332 const unsigned int current_face_index =
1336 const unsigned int n_q_points =
1337 compute_n_q_points<VectorizedArrayType>(
1338 n_q_points_unvectorized[current_face_index]);
1339 store_unit_points_faces(unit_points_index[current_face_index],
1341 n_q_points_unvectorized[current_face_index],
1342 quadrature_on_face.get_points());
1346 update_flags_mapping,
1350 internal_mapping_data,
1358 cell->diameter(), mapping_data.inverse_jacobians));
1362 mapping_data_first = mapping_data;
1367 cell_type.push_back(
1370 if (current_face_index > 0)
1373 const bool affine_cells =
1374 cell_type[current_face_index] <=
1376 cell_type[current_face_index - 1] <=
1382 1e4 / cell->diameter() *
1383 std::numeric_limits<double>::epsilon() * 1024.);
1387 const auto comparison_result =
1388 (!affine_cells || mapping_data.inverse_jacobians.empty() ||
1392 mapping_data.inverse_jacobians[0],
1398 comparison_result ==
1401 compressed_data_index_offsets.push_back(
1402 compressed_data_index_offsets.back());
1404 else if (first_set &&
1405 (cell_type[current_face_index] <=
1408 mapping_data.inverse_jacobians[0],
1411 double>::ComparisonResult::equal))
1413 compressed_data_index_offsets.push_back(0);
1417 const unsigned int n_compressed_data_last_cell =
1418 cell_type[current_face_index - 1] <=
1421 compute_n_q_points<Number>(
1422 n_q_points_unvectorized[current_face_index - 1]);
1424 compressed_data_index_offsets.push_back(
1425 compressed_data_index_offsets.back() +
1426 n_compressed_data_last_cell);
1430 compressed_data_index_offsets.push_back(0);
1433 mapping_data_previous_cell = mapping_data;
1435 const unsigned int n_q_points_data = compute_n_q_points<Number>(
1436 n_q_points_unvectorized[current_face_index]);
1437 store_mapping_data(data_index_offsets[current_face_index],
1439 n_q_points_unvectorized[current_face_index],
1441 quadrature_on_face.get_weights(),
1442 data_index_offsets[current_face_index],
1443 cell_type[current_face_index] <=
1448 if (cell_type[current_face_index] <=
1450 size_compressed_data = compressed_data_index_offsets.back() + 1;
1452 size_compressed_data =
1454 compressed_data_index_offsets.back() +
1457 if (do_cell_index_compression)
1458 cell_index_to_compressed_cell_index[cell->active_cell_index()] =
1466 jacobians[0].resize(size_compressed_data);
1467 jacobians[0].shrink_to_fit();
1471 inverse_jacobians[0].resize(size_compressed_data);
1472 inverse_jacobians[0].shrink_to_fit();
1475 state = State::faces_on_cells_in_vector;
1481 template <
int dim,
int spacedim,
typename Number>
1482 template <
typename CellIteratorType>
1485 const std::vector<std::pair<CellIteratorType, unsigned int>>
1486 &face_iterator_range_interior,
1491 do_cell_index_compression =
false;
1496 const unsigned int n_faces = quadrature_vector.size();
1498 std::distance(face_iterator_range_interior.begin(),
1499 face_iterator_range_interior.end()));
1501 n_q_points_unvectorized.reserve(n_faces);
1503 cell_type.reserve(n_faces);
1504 face_number.reserve(n_faces);
1507 unit_points_index.reserve(n_faces + 1);
1508 unit_points_index.push_back(0);
1509 data_index_offsets.reserve(n_faces + 1);
1510 data_index_offsets.push_back(0);
1511 for (
const auto &quadrature : quadrature_vector)
1513 const unsigned int n_points = quadrature.size();
1514 n_q_points_unvectorized.push_back(n_points);
1516 const unsigned int n_q_points =
1517 compute_n_q_points<VectorizedArrayType>(n_points);
1518 unit_points_index.push_back(unit_points_index.back() + n_q_points);
1520 const unsigned int n_q_points_data =
1521 compute_n_q_points<Number>(n_points);
1522 data_index_offsets.push_back(data_index_offsets.back() +
1526 const unsigned int n_unit_points = unit_points_index.back();
1527 const unsigned int n_data_points = data_index_offsets.back();
1530 resize_unit_points_faces(n_unit_points);
1531 resize_data_fields(n_data_points,
true);
1533 std::array<MappingData, 2> mapping_data;
1534 std::array<MappingData, 2> mapping_data_previous_cell;
1535 std::array<MappingData, 2> mapping_data_first;
1536 bool first_set =
false;
1537 unsigned int size_compressed_data = 0;
1538 unsigned int face_index = 0;
1539 for (
const auto &cell_and_f : face_iterator_range_interior)
1541 const auto &quadrature_on_face = quadrature_vector[face_index];
1542 const bool empty = quadrature_on_face.empty();
1545 const auto &cell_m = cell_and_f.first;
1546 const auto f_m = cell_and_f.second;
1549 const auto &cell_p =
1550 cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m);
1552 cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m);
1554 face_number.emplace_back(f_m, f_p);
1557 cell_m->combined_face_orientation(f_m) ==
1559 cell_p->combined_face_orientation(f_p) ==
1562 "Non standard face orientation is currently not implemented."));
1565 const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
1566 n_q_points_unvectorized[face_index]);
1567 store_unit_points_faces(unit_points_index[face_index],
1569 n_q_points_unvectorized[face_index],
1570 quadrature_on_face.get_points());
1575 update_flags_mapping,
1579 internal_mapping_data,
1585 update_flags_mapping,
1589 internal_mapping_data,
1599 cell_m->diameter(), mapping_data[0].inverse_jacobians),
1601 cell_m->diameter(), mapping_data[1].inverse_jacobians)));
1607 mapping_data_first = mapping_data;
1612 cell_type.push_back(
1618 const bool affine_cells =
1619 cell_type[face_index] <=
1621 cell_type[face_index - 1] <=
1627 1e4 / cell_m->diameter() *
1628 std::numeric_limits<double>::epsilon() * 1024.);
1632 const auto comparison_result_m =
1633 (!affine_cells || mapping_data[0].inverse_jacobians.empty() ||
1634 mapping_data_previous_cell[0].inverse_jacobians.empty()) ?
1637 mapping_data[0].inverse_jacobians[0],
1638 mapping_data_previous_cell[0].inverse_jacobians[0]);
1640 const auto comparison_result_p =
1641 (!affine_cells || mapping_data[1].inverse_jacobians.empty() ||
1642 mapping_data_previous_cell[1].inverse_jacobians.empty()) ?
1645 mapping_data[1].inverse_jacobians[0],
1646 mapping_data_previous_cell[1].inverse_jacobians[0]);
1651 comparison_result_m ==
1653 comparison_result_p ==
1656 compressed_data_index_offsets.push_back(
1657 compressed_data_index_offsets.back());
1659 else if (first_set &&
1660 (cell_type[face_index] <=
1663 mapping_data[0].inverse_jacobians[0],
1664 mapping_data_first[0].inverse_jacobians[0]) ==
1666 double>::ComparisonResult::equal) &&
1668 mapping_data[1].inverse_jacobians[0],
1669 mapping_data_first[1].inverse_jacobians[0]) ==
1672 compressed_data_index_offsets.push_back(0);
1676 const unsigned int n_compressed_data_last_cell =
1677 cell_type[face_index - 1] <=
1680 compute_n_q_points<Number>(
1681 n_q_points_unvectorized[face_index - 1]);
1683 compressed_data_index_offsets.push_back(
1684 compressed_data_index_offsets.back() +
1685 n_compressed_data_last_cell);
1689 compressed_data_index_offsets.push_back(0);
1692 mapping_data_previous_cell = mapping_data;
1694 const unsigned int n_q_points_data =
1695 compute_n_q_points<Number>(n_q_points_unvectorized[face_index]);
1698 store_mapping_data(data_index_offsets[face_index],
1700 n_q_points_unvectorized[face_index],
1702 quadrature_on_face.get_weights(),
1703 data_index_offsets[face_index],
1704 cell_type[face_index] <=
1710 store_mapping_data(data_index_offsets[face_index],
1712 n_q_points_unvectorized[face_index],
1714 quadrature_on_face.get_weights(),
1715 data_index_offsets[face_index],
1716 cell_type[face_index] <=
1722 if (cell_type[face_index] <=
1724 size_compressed_data = compressed_data_index_offsets.back() + 1;
1726 size_compressed_data =
1728 compressed_data_index_offsets.back() + n_q_points_data);
1735 jacobians[0].resize(size_compressed_data);
1736 jacobians[0].shrink_to_fit();
1737 jacobians[1].resize(size_compressed_data);
1738 jacobians[1].shrink_to_fit();
1742 inverse_jacobians[0].resize(size_compressed_data);
1743 inverse_jacobians[0].shrink_to_fit();
1744 inverse_jacobians[1].resize(size_compressed_data);
1745 inverse_jacobians[1].shrink_to_fit();
1748 state = State::face_vector;
1754 template <
int dim,
int spacedim,
typename Number>
1758 return state == State::faces_on_cells_in_vector;
1763 template <
int dim,
int spacedim,
typename Number>
1766 const unsigned int geometry_index)
const
1768 return n_q_points_unvectorized[geometry_index];
1772 template <
int dim,
int spacedim,
typename Number>
1775 const unsigned int geometry_index)
const
1777 return cell_type[geometry_index];
1782 template <
int dim,
int spacedim,
typename Number>
1788 additional_data.store_cells,
1790 "Cells have been not stored. You can enable this by Additional::store_cells."));
1798 template <
int dim,
int spacedim,
typename Number>
1799 template <
typename NumberType>
1802 const unsigned int n_q_points_unvectorized)
1804 const unsigned int n_lanes =
1806 const unsigned int n_filled_lanes_last_batch =
1807 n_q_points_unvectorized % n_lanes;
1808 unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
1809 if (n_filled_lanes_last_batch > 0)
1816 template <
int dim,
int spacedim,
typename Number>
1817 template <
bool is_face>
1821 const unsigned int face_number)
const
1826 const unsigned int compressed_cell_index =
1830 Assert(state == State::cell_vector,
1832 "This mapping info is not reinitialized for a cell vector!"));
1833 return compressed_cell_index;
1839 "cell_index has to be set if face_number is specified!"));
1840 Assert(state == State::faces_on_cells_in_vector ||
1841 state == State::face_vector,
1842 ExcMessage(
"This mapping info is not reinitialized for faces"
1843 " on cells in a vector!"));
1844 if (state == State::faces_on_cells_in_vector)
1845 return cell_index_offset[compressed_cell_index] + face_number;
1846 else if (state == State::face_vector)
1853 template <
int dim,
int spacedim,
typename Number>
1858 if (do_cell_index_compression)
1862 ExcMessage(
"Mapping info object was not initialized for this"
1863 " active cell index!"));
1864 return cell_index_to_compressed_cell_index[
cell_index];
1871 template <
int dim,
int spacedim,
typename Number>
1874 const unsigned int unit_points_index_offset,
1875 const unsigned int n_q_points,
1876 const unsigned int n_q_points_unvectorized,
1879 const unsigned int n_lanes =
1882 for (
unsigned int q = 0; q < n_q_points; ++q)
1884 const unsigned int offset = unit_points_index_offset + q;
1885 for (
unsigned int v = 0;
1886 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1888 for (
unsigned int d = 0; d < dim; ++d)
1890 unit_points[offset][d], v) = points[q * n_lanes + v][d];
1896 template <
int dim,
int spacedim,
typename Number>
1899 const unsigned int unit_points_index_offset,
1900 const unsigned int n_q_points,
1901 const unsigned int n_q_points_unvectorized,
1904 const unsigned int n_lanes =
1907 for (
unsigned int q = 0; q < n_q_points; ++q)
1909 const unsigned int offset = unit_points_index_offset + q;
1910 for (
unsigned int v = 0;
1911 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1913 for (
unsigned int d = 0; d < dim - 1; ++d)
1915 unit_points_faces[offset][d], v) = points[q * n_lanes + v][d];
1921 template <
int dim,
int spacedim,
typename Number>
1924 const unsigned int unit_points_index_offset,
1925 const unsigned int n_q_points,
1926 const unsigned int n_q_points_unvectorized,
1928 const std::vector<double> &weights,
1929 const unsigned int compressed_unit_point_index_offset,
1930 const bool affine_cell,
1931 const bool is_interior)
1933 const unsigned int n_lanes =
1936 for (
unsigned int q = 0; q < n_q_points; ++q)
1938 const unsigned int offset = unit_points_index_offset + q;
1939 const unsigned int compressed_offset =
1940 compressed_unit_point_index_offset + q;
1941 for (
unsigned int v = 0;
1942 v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
1945 if (q == 0 || !affine_cell)
1948 for (
unsigned int d = 0; d < dim; ++d)
1949 for (
unsigned int s = 0; s < spacedim; ++s)
1951 jacobians[is_interior ? 0 : 1][compressed_offset][d][s],
1952 v) = mapping_data.
jacobians[q * n_lanes + v][d][s];
1953 if (update_flags_mapping &
1955 for (
unsigned int d = 0; d < dim; ++d)
1956 for (
unsigned int s = 0; s < spacedim; ++s)
1958 inverse_jacobians[is_interior ? 0 : 1]
1959 [compressed_offset][d][s],
1968 if (additional_data.use_global_weights)
1971 JxW_values[offset], v) = weights[q * n_lanes + v];
1976 JxW_values[offset], v) =
1981 for (
unsigned int s = 0; s < spacedim; ++s)
1983 normal_vectors[offset][s], v) =
1985 if (update_flags_mapping &
1987 for (
unsigned int s = 0; s < spacedim; ++s)
1989 real_points[offset][s], v) =
1998 template <
int dim,
int spacedim,
typename Number>
2001 const unsigned int n_unit_point_batches)
2003 unit_points.resize(n_unit_point_batches);
2008 template <
int dim,
int spacedim,
typename Number>
2011 const unsigned int n_unit_point_batches)
2013 unit_points_faces.resize(n_unit_point_batches);
2018 template <
int dim,
int spacedim,
typename Number>
2021 const unsigned int n_data_point_batches,
2022 const bool is_face_centric)
2026 jacobians[0].resize(n_data_point_batches);
2027 if (is_face_centric)
2028 jacobians[1].resize(n_data_point_batches);
2032 inverse_jacobians[0].resize(n_data_point_batches);
2033 if (is_face_centric)
2034 inverse_jacobians[1].resize(n_data_point_batches);
2037 JxW_values.resize(n_data_point_batches);
2039 normal_vectors.resize(n_data_point_batches);
2041 real_points.resize(n_data_point_batches);
2046 template <
int dim,
int spacedim,
typename Number>
2051 const unsigned int offset)
const
2053 return unit_points.data() + offset;
2058 template <
int dim,
int spacedim,
typename Number>
2063 const unsigned int offset)
const
2065 return unit_points_faces.data() + offset;
2070 template <
int dim,
int spacedim,
typename Number>
2073 const unsigned int offset)
const
2075 return real_points.data() + offset;
2080 template <
int dim,
int spacedim,
typename Number>
2083 const unsigned int geometry_index)
const
2085 return unit_points_index[geometry_index];
2090 template <
int dim,
int spacedim,
typename Number>
2093 const unsigned int geometry_index)
const
2095 return data_index_offsets[geometry_index];
2099 template <
int dim,
int spacedim,
typename Number>
2102 const unsigned int geometry_index)
const
2104 return compressed_data_index_offsets[geometry_index];
2109 template <
int dim,
int spacedim,
typename Number>
2112 const bool is_interior)
const
2114 return jacobians[is_interior ? 0 : 1].data() + offset;
2119 template <
int dim,
int spacedim,
typename Number>
2122 const unsigned int offset,
2123 const bool is_interior)
const
2125 return inverse_jacobians[is_interior ? 0 : 1].data() + offset;
2130 template <
int dim,
int spacedim,
typename Number>
2133 const unsigned int offset)
const
2135 return normal_vectors.data() + offset;
2140 template <
int dim,
int spacedim,
typename Number>
2141 inline const Number *
2144 return JxW_values.data() + offset;
2149 template <
int dim,
int spacedim,
typename Number>
2158 template <
int dim,
int spacedim,
typename Number>
2162 return update_flags;
2167 template <
int dim,
int spacedim,
typename Number>
2171 return update_flags_mapping;
2176 template <
int dim,
int spacedim,
typename Number>
2177 boost::signals2::connection
2179 const std::function<
void()> &set_is_reinitialized)
2181 return is_reinitialized.connect(set_is_reinitialized);
2186 template <
int dim,
int spacedim,
typename Number>
2193 memory += cell_type.capacity() *
2206 cell_index_to_compressed_cell_index);
2208 memory +=
sizeof(*this);
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, 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
boost::signals2::connection connect_is_reinitialized(const std::function< void()> &set_is_reinitialized)
const Mapping< dim, spacedim > & get_mapping() const
void reinit_faces(const ContainerType &cell_iterator_range, const std::vector< std::vector< Quadrature< dim - 1 > > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
std::array< AlignedVector< DerivativeForm< 1, spacedim, dim, Number > >, 2 > inverse_jacobians
SmartPointer< const Triangulation< dim, spacedim > > triangulation
const Point< dim, VectorizedArrayType > * get_unit_point(const unsigned int offset) const
unsigned int get_n_q_points_unvectorized(const unsigned int geometry_index) const
unsigned int compute_compressed_cell_index(const unsigned int cell_index) const
const Point< dim - 1, VectorizedArrayType > * get_unit_point_faces(const unsigned int offset) const
Quadrature< dim > quadrature
std::vector< unsigned int > cell_index_offset
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Quadrature< dim > &quadrature)
void store_mapping_data(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const MappingData &mapping_data, const std::vector< double > &weights, const unsigned int compressed_unit_point_index_offset, const bool affine_cell, const bool is_interior=true)
const SmartPointer< const Mapping< dim, spacedim > > mapping
AlignedVector< Tensor< 1, spacedim, Number > > normal_vectors
unsigned int compute_geometry_index_offset(const unsigned int cell_index, const unsigned int face_number) const
std::vector<::internal::MatrixFreeFunctions::GeometryType > cell_type
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim > > &unit_points)
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< Quadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
MappingInfo & operator=(const MappingInfo &)=delete
boost::signals2::signal< void()> is_reinitialized
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > internal_mapping_data
void reinit_surface(const ContainerType &cell_iterator_range, const std::vector< ImmersedSurfaceQuadrature< dim > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_unit_points_faces(const unsigned int n_unit_point_batches)
const Tensor< 1, spacedim, Number > * get_normal_vector(const unsigned int offset) const
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const std::vector< Point< dim > > &unit_points)
const DerivativeForm< 1, dim, spacedim, Number > * get_jacobian(const unsigned int offset, const bool is_interior=true) const
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
void reinit_cells(const ContainerType &cell_iterator_range, const std::vector< std::vector< Point< dim > > > &unit_points_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
void resize_data_fields(const unsigned int n_data_point_batches, const bool is_face_centric=false)
void resize_unit_points(const unsigned int n_unit_point_batches)
std::size_t memory_consumption() const
UpdateFlags get_update_flags() const
void store_unit_points_faces(const unsigned int unit_points_index_offset, const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const std::vector< Point< dim - 1 > > &points)
unsigned int compute_unit_point_index_offset(const unsigned int geometry_index) const
std::vector< unsigned int > cell_index_to_compressed_cell_index
AlignedVector< Point< spacedim, Number > > real_points
unsigned int compute_data_index_offset(const unsigned int geometry_index) const
@ 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 UpdateFlags required_update_flags(const SmartPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags)
static void compute_mapping_data_for_immersed_surface_quadrature(const SmartPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ImmersedSurfaceQuadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_quadrature(const SmartPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, CellSimilarity::Similarity &cell_similarity, const Quadrature< dim > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
static void compute_mapping_data_for_face_quadrature(const SmartPointer< const Mapping< dim, spacedim > > &mapping, const UpdateFlags &update_flags_mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > &internal_mapping_data, MappingData &mapping_data)
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)