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);
191 std::pair<unsigned int, RefinementCase<2>>
192 select_subface_no_and_refinement_case(
193 const unsigned int subface_no,
194 const bool face_orientation,
195 const bool face_flip,
196 const bool face_rotation,
199 constexpr int dim = 3;
269 static const unsigned int
305 equivalent_refine_case[ref_case][subface_no];
306 const unsigned int equ_subface_no =
307 equivalent_subface_number[ref_case][subface_no];
314 (face_orientation == face_rotation ?
315 ref_case_permutation[equ_ref_case] :
318 const unsigned int final_subface_no =
328 return std::make_pair(final_subface_no, final_ref_case);
359 append_subobject_rule(
363 const double measure,
366 std::vector<double> &weights)
368 const auto support_points =
372 combined_orientation));
374 for (
unsigned int j = 0; j < quadrature.
size(); ++j)
379 for (
const unsigned int vertex_no :
382 support_points[vertex_no] *
386 points.push_back(mapped_point);
390 weights.push_back(quadrature.
weight(j) * measure /
391 face_reference_cell.
volume());
404 const unsigned int face_no,
408 const auto face_quadrature =
413 q_points = face_quadrature.get_points();
422 const unsigned int face_no,
426 const auto face_quadrature =
431 q_points = face_quadrature.get_points();
440 const unsigned int face_no,
444 (void)reference_cell;
450 internal::QProjector::project_to_hex_face_and_append(quadrature.
get_points(),
461 const unsigned int face_no,
462 const bool face_orientation,
463 const bool face_flip,
464 const bool face_rotation)
482 const unsigned int face_no,
487 reference_cell.n_face_orientations(face_no));
490 std::vector<Point<dim>> points;
491 std::vector<double> weights;
494 reference_cell.face_reference_cell(face_no);
495 std::vector<Point<dim>> face_vertices(face_reference_cell.
n_vertices());
496 for (
const unsigned int vertex_no : face_reference_cell.
vertex_indices())
497 face_vertices[vertex_no] =
498 reference_cell.face_vertex_location<dim>(face_no, vertex_no);
499 internal::QProjector::append_subobject_rule(face_reference_cell,
502 reference_cell.face_measure(
504 combined_orientation,
517 const unsigned int face_no,
523 (void)reference_cell;
525 const unsigned int dim = 1;
529 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
538 const unsigned int face_no,
539 const unsigned int subface_no,
544 const auto face_quadrature =
545 project_to_subface(reference_cell,
551 q_points = face_quadrature.get_points();
560 const unsigned int face_no,
561 const unsigned int subface_no,
566 (void)reference_cell;
568 const auto face_quadrature =
569 project_to_subface(reference_cell,
575 q_points = face_quadrature.get_points();
585 const unsigned int face_no,
586 const unsigned int subface_no,
607 const unsigned int face_no,
608 const unsigned int subface_no,
614 reference_cell.n_face_orientations(face_no));
617 reference_cell.face_reference_cell(face_no)
618 .template n_children<dim - 1>(ref_case));
620 std::vector<Point<dim>> q_points;
621 std::vector<double> q_weights = quadrature.
get_weights();
622 q_points.reserve(quadrature.
size());
624 if constexpr (dim == 1)
627 q_points.emplace_back(
static_cast<double>(face_no));
629 else if constexpr (dim == 2)
634 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
639 q_points.emplace_back(quadrature.
point(p)[0] / 2, 0);
641 q_points.emplace_back(0.5 + quadrature.
point(p)[0] / 2, 0);
643 else if (face_no == 1)
646 q_points.emplace_back(1 - quadrature.
point(p)[0] / 2,
647 quadrature.
point(p)[0] / 2);
649 q_points.emplace_back(0.5 - quadrature.
point(p)[0] / 2,
650 0.5 + quadrature.
point(p)[0] / 2);
652 else if (face_no == 2)
655 q_points.emplace_back(0, 1 - quadrature.
point(p)[0] / 2);
657 q_points.emplace_back(0, 0.5 - quadrature.
point(p)[0] / 2);
663 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
668 q_points.emplace_back(0, quadrature.
point(p)[0] / 2);
670 q_points.emplace_back(0, quadrature.
point(p)[0] / 2 + 0.5);
672 else if (face_no == 1)
675 q_points.emplace_back(1, quadrature.
point(p)[0] / 2);
677 q_points.emplace_back(1, quadrature.
point(p)[0] / 2 + 0.5);
679 else if (face_no == 2)
682 q_points.emplace_back(quadrature.
point(p)[0] / 2, 0);
684 q_points.emplace_back(quadrature.
point(p)[0] / 2 + 0.5, 0);
686 else if (face_no == 3)
689 q_points.emplace_back(quadrature.
point(p)[0] / 2, 1);
691 q_points.emplace_back(quadrature.
point(p)[0] / 2 + 0.5, 1);
701 std::reverse(q_points.begin(), q_points.end());
702 std::reverse(q_weights.begin(), q_weights.end());
704 for (
auto &w : q_weights)
705 w *= reference_cell.face_measure(face_no);
707 else if constexpr (dim == 3)
710 internal::QProjector::project_to_hex_face_and_append(
711 quadrature.
get_points(), face_no, q_points, ref_case, subface_no);
729 std::vector<Point<dim>> points;
730 std::vector<double> weights;
732 for (
const unsigned int face_no : reference_cell.face_indices())
735 reference_cell.face_reference_cell(face_no);
736 std::vector<Point<dim>> face_vertices(face_reference_cell.
n_vertices());
737 for (
const unsigned int vertex_no : face_reference_cell.
vertex_indices())
738 face_vertices[vertex_no] =
739 reference_cell.face_vertex_location<dim>(face_no, vertex_no);
742 combined_orientation < reference_cell.n_face_orientations(face_no);
743 ++combined_orientation)
744 internal::QProjector::append_subobject_rule(
746 quadrature[quadrature.
size() == 1 ? 0 : face_no],
748 reference_cell.face_measure(face_no),
749 combined_orientation,
765 (void)reference_cell;
767 const unsigned int dim = 1;
774 std::vector<Point<dim>> q_points;
775 q_points.reserve(n_points * n_faces * subfaces_per_face);
776 std::vector<Point<dim>> help(n_points);
780 for (
unsigned int face = 0; face < n_faces; ++face)
781 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
783 project_to_subface(reference_cell, quadrature, face, subface, help);
784 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
788 std::vector<double> weights;
789 weights.reserve(n_points * n_faces * subfaces_per_face);
790 for (
unsigned int face = 0; face < n_faces; ++face)
791 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
794 std::back_inserter(weights));
796 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
798 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
815 const unsigned int dim = 2;
817 std::vector<Point<dim>> q_points;
818 std::vector<double> weights;
822 for (
unsigned int face = 0; face < reference_cell.n_faces(); ++face)
824 orientation < reference_cell.n_face_orientations(face);
826 for (
unsigned int subface = 0;
828 reference_cell.face_reference_cell(face).n_isotropic_children();
831 const unsigned int local_subface =
834 const auto sub_quadrature =
835 project_to_subface(reference_cell,
841 q_points.insert(q_points.end(),
842 sub_quadrature.get_points().begin(),
843 sub_quadrature.get_points().end());
844 weights.insert(weights.end(),
845 sub_quadrature.get_weights().begin(),
846 sub_quadrature.get_weights().end());
865 const unsigned int dim = 3;
867 const unsigned int n_points = quadrature.
size(),
869 total_subfaces_per_face = 2 + 2 + 4;
872 std::vector<Point<dim>> q_points;
873 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
875 std::vector<double> weights;
876 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
880 for (
unsigned char offset = 0; offset < 8; ++offset)
882 const auto mutation =
883 internal::QProjector::mutate_points_with_offset(quadrature.
get_points(),
887 for (
unsigned int face = 0; face < n_faces; ++face)
891 for (
unsigned int subface = 0;
896 internal::QProjector::project_to_hex_face_and_append(
906 std::back_inserter(weights));
910 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
912 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
924 const unsigned int child_no)
926 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
928 (void)reference_cell;
932 const unsigned int n_q_points = quadrature.
size();
934 std::vector<Point<dim>> q_points(n_q_points);
935 for (
unsigned int i = 0; i < n_q_points; ++i)
943 std::vector<double> weights = quadrature.
get_weights();
944 for (
unsigned int i = 0; i < n_q_points; ++i)
957 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
959 (void)reference_cell;
961 const unsigned int n_points = quadrature.
size(),
964 std::vector<Point<dim>> q_points(n_points * n_children);
965 std::vector<double> weights(n_points * n_children);
969 for (
unsigned int child = 0; child < n_children; ++child)
972 project_to_child(reference_cell, quadrature, child);
973 for (
unsigned int i = 0; i < n_points; ++i)
975 q_points[child * n_points + i] = help.
point(i);
976 weights[child * n_points + i] = help.
weight(i);
991 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
993 (void)reference_cell;
995 const unsigned int n = quadrature.
size();
996 std::vector<Point<dim>> points(n);
997 std::vector<double> weights(n);
998 const double length = p1.
distance(p2);
1000 for (
unsigned int k = 0; k < n; ++k)
1002 const double alpha = quadrature.
point(k)[0];
1003 points[k] = alpha * p2;
1004 points[k] += (1. - alpha) * p1;
1005 weights[k] = length * quadrature.
weight(k);
1015 const unsigned int face_no,
1016 const bool face_orientation,
1017 const bool face_flip,
1018 const bool face_rotation,
1019 const unsigned int n_quadrature_points)
1021 return face(reference_cell,
1026 n_quadrature_points);
1035 const unsigned int face_no,
1037 const unsigned int n_quadrature_points)
1041 reference_cell.n_face_orientations(face_no));
1045 return {(reference_cell.n_face_orientations(face_no) * face_no +
1046 combined_orientation) *
1047 n_quadrature_points};
1056 const unsigned int face_no,
1057 const bool face_orientation,
1058 const bool face_flip,
1059 const bool face_rotation,
1062 return face(reference_cell,
1076 const unsigned int face_no,
1082 reference_cell.n_face_orientations(face_no));
1085 unsigned int offset = 0;
1086 if (quadrature.
size() == 1)
1088 reference_cell.n_face_orientations(0) * quadrature[0].
size() * face_no;
1090 for (
unsigned int i = 0; i < face_no; ++i)
1091 offset += reference_cell.n_face_orientations(i) * quadrature[i].
size();
1093 return {offset + combined_orientation *
1094 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1103 const unsigned int face_no,
1104 const unsigned int subface_no,
1105 const bool face_orientation,
1106 const bool face_flip,
1107 const bool face_rotation,
1108 const unsigned int n_quadrature_points,
1118 n_quadrature_points,
1128 const unsigned int face_no,
1129 const unsigned int subface_no,
1131 const unsigned int n_quadrature_points,
1135 (void)reference_cell;
1142 n_quadrature_points);
1151 const unsigned int face_no,
1152 const unsigned int subface_no,
1154 const unsigned int n_quadrature_points,
1161 const unsigned int n_faces = reference_cell.n_faces();
1162 const unsigned int n_subfaces =
1163 reference_cell.face_reference_cell(face_no).n_isotropic_children();
1164 const unsigned int n_orientations =
1165 reference_cell.n_face_orientations(face_no);
1171 return ((face_no * n_orientations * n_subfaces +
1172 combined_orientation * n_subfaces + subface_no) *
1173 n_quadrature_points);
1182 const unsigned int face_no,
1183 const unsigned int subface_no,
1185 const unsigned int n_quadrature_points,
1188 const unsigned int dim = 3;
1190 const auto [face_orientation, face_rotation, face_flip] =
1194 (void)reference_cell;
1208 const unsigned int total_subfaces_per_face = 8;
1225 static const unsigned int ref_case_offset[3] = {
1231 const std::pair<unsigned int, RefinementCase<2>>
1232 final_subface_no_and_ref_case =
1233 internal::QProjector::select_subface_no_and_refinement_case(
1234 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1236 return ((face_no * total_subfaces_per_face +
1237 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1238 final_subface_no_and_ref_case.first) +
1240 combined_orientation) *
1241 n_quadrature_points;
1250 const unsigned int face_no)
1252 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1254 (void)reference_cell;
1256 std::vector<Point<dim>> points(quadrature.
size());
1267 const unsigned int face_no,
1268 const unsigned int subface_no,
1271 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1273 (void)reference_cell;
1275 std::vector<Point<dim>> points(quadrature.
size());
1277 reference_cell, quadrature, face_no, subface_no, points, ref_case);
auto make_const_array_view(const Container &container) -> decltype(make_array_view(container))
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
unsigned int n_vertices() 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)