22#include <boost/container/small_vector.hpp>
36 reflect(
const std::vector<
Point<2>> &points)
38 std::vector<Point<2>> q_points;
39 q_points.reserve(points.size());
41 q_points.emplace_back(p[1], p[0]);
50 rotate(
const std::vector<
Point<2>> &points,
const unsigned int n_times)
52 std::vector<Point<2>> q_points;
53 q_points.reserve(points.size());
59 q_points.push_back(p);
64 q_points.emplace_back(1.0 - p[1], p[0]);
69 q_points.emplace_back(1.0 - p[0], 1.0 - p[1]);
74 q_points.emplace_back(p[1], 1.0 - p[0]);
88 project_to_hex_face_and_append(
90 const unsigned int face_no,
93 const unsigned int subface_no = 0)
97 const double const_value = face_no % 2;
102 const unsigned int xi_index = (1 + face_no / 2) % 3,
103 eta_index = (2 + face_no / 2) % 3,
104 const_index = face_no / 2;
111 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
112 eta_translation = 0.0;
121 xi_translation = subface_no % 2 * 0.5;
125 eta_translation = subface_no % 2 * 0.5;
130 xi_translation =
int(subface_no % 2) * 0.5;
131 eta_translation =
int(subface_no / 2) * 0.5;
143 cell_point[xi_index] = xi_scale * p[0] + xi_translation;
144 cell_point[eta_index] = eta_scale * p[1] + eta_translation;
145 cell_point[const_index] = const_value;
146 q_points.push_back(cell_point);
150 std::vector<Point<2>>
151 mutate_points_with_offset(
152 const std::vector<
Point<2>> &points,
165 switch (combined_orientation)
168 return reflect(points);
170 return rotate(reflect(points), 3);
172 return rotate(reflect(points), 2);
174 return rotate(reflect(points), 1);
178 return rotate(points, 1);
180 return rotate(points, 2);
182 return rotate(points, 3);
194 combined_orientation),
198 std::pair<unsigned int, RefinementCase<2>>
199 select_subface_no_and_refinement_case(
200 const unsigned int subface_no,
201 const bool face_orientation,
202 const bool face_flip,
203 const bool face_rotation,
206 constexpr int dim = 3;
276 static const unsigned int
312 equivalent_refine_case[ref_case][subface_no];
313 const unsigned int equ_subface_no =
314 equivalent_subface_number[ref_case][subface_no];
321 (face_orientation == face_rotation ?
322 ref_case_permutation[equ_ref_case] :
325 const unsigned int final_subface_no =
335 return std::make_pair(final_subface_no, final_ref_case);
347 const unsigned int face_no,
351 const auto face_quadrature =
356 q_points = face_quadrature.get_points();
365 const unsigned int face_no,
369 const auto face_quadrature =
374 q_points = face_quadrature.get_points();
383 const unsigned int face_no,
387 (void)reference_cell;
393 internal::QProjector::project_to_hex_face_and_append(quadrature.
get_points(),
404 const unsigned int face_no,
405 const bool face_orientation,
406 const bool face_flip,
407 const bool face_rotation)
424 const unsigned int face_no,
431 std::vector<Point<dim>> q_points;
432 std::vector<double> q_weights = quadrature.
get_weights();
433 q_points.reserve(quadrature.
size());
434 if constexpr (dim == 1)
437 q_points.emplace_back(
static_cast<double>(face_no));
439 else if constexpr (dim == 2)
444 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
447 q_points.emplace_back(quadrature.
point(p)[0], 0);
448 else if (face_no == 1)
449 q_points.emplace_back(1.0 - quadrature.
point(p)[0],
450 quadrature.
point(p)[0]);
451 else if (face_no == 2)
452 q_points.emplace_back(0, 1.0 - quadrature.
point(p)[0]);
457 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
460 q_points.emplace_back(0, quadrature.
point(p)[0]);
461 else if (face_no == 1)
462 q_points.emplace_back(1, quadrature.
point(p)[0]);
463 else if (face_no == 2)
464 q_points.emplace_back(quadrature.
point(p)[0], 0);
465 else if (face_no == 3)
466 q_points.emplace_back(quadrature.
point(p)[0], 1);
475 std::reverse(q_points.begin(), q_points.end());
476 std::reverse(q_weights.begin(), q_weights.end());
479 else if constexpr (dim == 3)
483 internal::QProjector::mutate_quadrature(quadrature, orientation);
500 const unsigned int face_no,
506 (void)reference_cell;
508 const unsigned int dim = 1;
512 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
521 const unsigned int face_no,
522 const unsigned int subface_no,
527 const auto face_quadrature =
528 project_to_subface(reference_cell,
534 q_points = face_quadrature.get_points();
543 const unsigned int face_no,
544 const unsigned int subface_no,
549 (void)reference_cell;
551 const auto face_quadrature =
552 project_to_subface(reference_cell,
558 q_points = face_quadrature.get_points();
568 const unsigned int face_no,
569 const unsigned int subface_no,
590 const unsigned int face_no,
591 const unsigned int subface_no,
599 reference_cell.face_reference_cell(face_no)
600 .template n_children<dim - 1>(ref_case));
602 std::vector<Point<dim>> q_points;
603 std::vector<double> q_weights = quadrature.
get_weights();
604 q_points.reserve(quadrature.
size());
606 if constexpr (dim == 1)
609 q_points.emplace_back(
static_cast<double>(face_no));
611 else if constexpr (dim == 2)
616 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
621 q_points.emplace_back(quadrature.
point(p)[0] / 2, 0);
623 q_points.emplace_back(0.5 + quadrature.
point(p)[0] / 2, 0);
625 else if (face_no == 1)
628 q_points.emplace_back(1 - quadrature.
point(p)[0] / 2,
629 quadrature.
point(p)[0] / 2);
631 q_points.emplace_back(0.5 - quadrature.
point(p)[0] / 2,
632 0.5 + quadrature.
point(p)[0] / 2);
634 else if (face_no == 2)
637 q_points.emplace_back(0, 1 - quadrature.
point(p)[0] / 2);
639 q_points.emplace_back(0, 0.5 - quadrature.
point(p)[0] / 2);
645 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
650 q_points.emplace_back(0, quadrature.
point(p)[0] / 2);
652 q_points.emplace_back(0, quadrature.
point(p)[0] / 2 + 0.5);
654 else if (face_no == 1)
657 q_points.emplace_back(1, quadrature.
point(p)[0] / 2);
659 q_points.emplace_back(1, quadrature.
point(p)[0] / 2 + 0.5);
661 else if (face_no == 2)
664 q_points.emplace_back(quadrature.
point(p)[0] / 2, 0);
666 q_points.emplace_back(quadrature.
point(p)[0] / 2 + 0.5, 0);
668 else if (face_no == 3)
671 q_points.emplace_back(quadrature.
point(p)[0] / 2, 1);
673 q_points.emplace_back(quadrature.
point(p)[0] / 2 + 0.5, 1);
683 std::reverse(q_points.begin(), q_points.end());
684 std::reverse(q_weights.begin(), q_weights.end());
687 else if constexpr (dim == 3)
690 internal::QProjector::project_to_hex_face_and_append(
691 quadrature.
get_points(), face_no, q_points, ref_case, subface_no);
709 const auto process = [&](
const std::vector<std::vector<Point<dim>>> &faces) {
711 std::vector<Point<dim>> points;
712 std::vector<double> weights;
715 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
718 reference_cell.face_reference_cell(face_no);
722 orientation < reference_cell.n_face_orientations(face_no);
725 const auto &face = faces[face_no];
746 const boost::container::small_vector<Point<dim>, 8> support_points =
753 const auto &sub_quadrature_points =
754 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
755 const auto &sub_quadrature_weights =
756 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
759 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
764 for (
const unsigned int i :
766 mapped_point += support_points[i] *
768 sub_quadrature_points[j], i);
770 points.push_back(mapped_point);
774 const double scaling = reference_cell.face_measure(face_no) /
775 face_reference_cell.
volume();
776 weights.push_back(sub_quadrature_weights[j] * scaling);
785 std::vector<std::vector<Point<dim>>> face_vertex_locations(
786 reference_cell.n_faces());
787 for (
const unsigned int f : reference_cell.face_indices())
789 face_vertex_locations[f].resize(
790 reference_cell.face_reference_cell(f).n_vertices());
791 for (
const unsigned int v :
792 reference_cell.face_reference_cell(f).vertex_indices())
793 face_vertex_locations[f][v] =
794 reference_cell.face_vertex_location<dim>(f, v);
797 return process(face_vertex_locations);
808 (void)reference_cell;
810 const unsigned int dim = 1;
817 std::vector<Point<dim>> q_points;
818 q_points.reserve(n_points * n_faces * subfaces_per_face);
819 std::vector<Point<dim>> help(n_points);
823 for (
unsigned int face = 0; face < n_faces; ++face)
824 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
826 project_to_subface(reference_cell, quadrature, face, subface, help);
827 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
831 std::vector<double> weights;
832 weights.reserve(n_points * n_faces * subfaces_per_face);
833 for (
unsigned int face = 0; face < n_faces; ++face)
834 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
837 std::back_inserter(weights));
839 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
841 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
860 const unsigned int dim = 2;
862 const unsigned int n_points = quadrature.
size(),
868 std::vector<Point<dim>> q_points;
869 q_points.reserve(n_points * n_faces * subfaces_per_face);
870 std::vector<Point<dim>> help(n_points);
874 for (
unsigned int face = 0; face < n_faces; ++face)
875 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
877 project_to_subface(reference_cell, quadrature, face, subface, help);
878 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
882 std::vector<double> weights;
883 weights.reserve(n_points * n_faces * subfaces_per_face);
884 for (
unsigned int face = 0; face < n_faces; ++face)
885 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
888 std::back_inserter(weights));
890 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
892 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
911 const unsigned int dim = 3;
913 const unsigned int n_points = quadrature.
size(),
915 total_subfaces_per_face = 2 + 2 + 4;
918 std::vector<Point<dim>> q_points;
919 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
921 std::vector<double> weights;
922 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
926 for (
unsigned char offset = 0; offset < 8; ++offset)
928 const auto mutation =
929 internal::QProjector::mutate_points_with_offset(quadrature.
get_points(),
933 for (
unsigned int face = 0; face < n_faces; ++face)
937 for (
unsigned int subface = 0;
942 internal::QProjector::project_to_hex_face_and_append(
952 std::back_inserter(weights));
956 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
958 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
970 const unsigned int child_no)
972 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
974 (void)reference_cell;
978 const unsigned int n_q_points = quadrature.
size();
980 std::vector<Point<dim>> q_points(n_q_points);
981 for (
unsigned int i = 0; i < n_q_points; ++i)
989 std::vector<double> weights = quadrature.
get_weights();
990 for (
unsigned int i = 0; i < n_q_points; ++i)
1003 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1005 (void)reference_cell;
1007 const unsigned int n_points = quadrature.
size(),
1010 std::vector<Point<dim>> q_points(n_points * n_children);
1011 std::vector<double> weights(n_points * n_children);
1015 for (
unsigned int child = 0; child < n_children; ++child)
1018 project_to_child(reference_cell, quadrature, child);
1019 for (
unsigned int i = 0; i < n_points; ++i)
1021 q_points[child * n_points + i] = help.
point(i);
1022 weights[child * n_points + i] = help.
weight(i);
1037 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1039 (void)reference_cell;
1041 const unsigned int n = quadrature.
size();
1042 std::vector<Point<dim>> points(n);
1043 std::vector<double> weights(n);
1044 const double length = p1.
distance(p2);
1046 for (
unsigned int k = 0; k < n; ++k)
1048 const double alpha = quadrature.
point(k)[0];
1049 points[k] = alpha * p2;
1050 points[k] += (1. - alpha) * p1;
1051 weights[k] = length * quadrature.
weight(k);
1061 const unsigned int face_no,
1062 const bool face_orientation,
1063 const bool face_flip,
1064 const bool face_rotation,
1065 const unsigned int n_quadrature_points)
1067 return face(reference_cell,
1072 n_quadrature_points);
1081 const unsigned int face_no,
1083 const unsigned int n_quadrature_points)
1087 reference_cell.n_face_orientations(face_no));
1091 return {(reference_cell.n_face_orientations(face_no) * face_no +
1092 combined_orientation) *
1093 n_quadrature_points};
1102 const unsigned int face_no,
1103 const bool face_orientation,
1104 const bool face_flip,
1105 const bool face_rotation,
1108 return face(reference_cell,
1122 const unsigned int face_no,
1128 reference_cell.n_face_orientations(face_no));
1131 unsigned int offset = 0;
1132 if (quadrature.
size() == 1)
1134 reference_cell.n_face_orientations(0) * quadrature[0].
size() * face_no;
1136 for (
unsigned int i = 0; i < face_no; ++i)
1137 offset += reference_cell.n_face_orientations(i) * quadrature[i].
size();
1139 return {offset + combined_orientation *
1140 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1149 const unsigned int face_no,
1150 const unsigned int subface_no,
1151 const bool face_orientation,
1152 const bool face_flip,
1153 const bool face_rotation,
1154 const unsigned int n_quadrature_points,
1164 n_quadrature_points,
1174 const unsigned int face_no,
1175 const unsigned int subface_no,
1177 const unsigned int n_quadrature_points,
1181 (void)reference_cell;
1188 n_quadrature_points);
1197 const unsigned int face_no,
1198 const unsigned int subface_no,
1200 const unsigned int n_quadrature_points,
1204 (void)reference_cell;
1211 n_quadrature_points);
1220 const unsigned int face_no,
1221 const unsigned int subface_no,
1223 const unsigned int n_quadrature_points,
1226 const unsigned int dim = 3;
1228 const auto [face_orientation, face_rotation, face_flip] =
1232 (void)reference_cell;
1246 const unsigned int total_subfaces_per_face = 8;
1263 static const unsigned int ref_case_offset[3] = {
1269 const std::pair<unsigned int, RefinementCase<2>>
1270 final_subface_no_and_ref_case =
1271 internal::QProjector::select_subface_no_and_refinement_case(
1272 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1274 return ((face_no * total_subfaces_per_face +
1275 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1276 final_subface_no_and_ref_case.first) +
1278 combined_orientation) *
1279 n_quadrature_points;
1288 const unsigned int face_no)
1290 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1292 (void)reference_cell;
1294 std::vector<Point<dim>> points(quadrature.
size());
1305 const unsigned int face_no,
1306 const unsigned int subface_no,
1309 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1311 (void)reference_cell;
1313 std::vector<Point<dim>> points(quadrature.
size());
1315 reference_cell, quadrature, face_no, subface_no, points, ref_case);
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &quadrature)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
static Quadrature< dim > project_to_oriented_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const internal::SubfaceCase< dim > ref_case)
static Quadrature< dim > project_to_oriented_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation)
static Quadrature< dim > project_to_line(const ReferenceCell &reference_cell, const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static Quadrature< dim > project_to_all_children(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) const
types::geometric_orientation get_inverse_combined_orientation(const types::geometric_orientation orientation) const
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const types::geometric_orientation orientation) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
constexpr ReferenceCell Line
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr unsigned int invalid_unsigned_int
constexpr types::geometric_orientation reverse_line_orientation
constexpr types::geometric_orientation default_geometric_orientation
unsigned char geometric_orientation
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)