34 reflect(
const std::vector<
Point<2>> &points)
37 std::vector<Point<2>> q_points;
38 q_points.reserve(points.size());
40 q_points.emplace_back(p[1], p[0]);
47 rotate(
const std::vector<
Point<2>> &points,
const unsigned int n_times)
49 std::vector<Point<2>> q_points;
50 q_points.reserve(points.size());
56 q_points.push_back(p);
61 q_points.emplace_back(1.0 - p[1], p[0]);
66 q_points.emplace_back(1.0 - p[0], 1.0 - p[1]);
71 q_points.emplace_back(p[1], 1.0 - p[0]);
85 project_to_hex_face_and_append(
87 const unsigned int face_no,
90 const unsigned int subface_no = 0)
94 const double const_value = face_no % 2;
99 const unsigned int xi_index = (1 + face_no / 2) % 3,
100 eta_index = (2 + face_no / 2) % 3,
101 const_index = face_no / 2;
108 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
109 eta_translation = 0.0;
118 xi_translation = subface_no % 2 * 0.5;
122 eta_translation = subface_no % 2 * 0.5;
127 xi_translation =
int(subface_no % 2) * 0.5;
128 eta_translation =
int(subface_no / 2) * 0.5;
140 cell_point[xi_index] = xi_scale * p(0) + xi_translation;
141 cell_point[eta_index] = eta_scale * p(1) + eta_translation;
142 cell_point[const_index] = const_value;
143 q_points.push_back(cell_point);
147 std::vector<Point<2>>
148 mutate_points_with_offset(
const std::vector<
Point<2>> &points,
149 const unsigned int offset)
158 return rotate(points, offset);
160 return reflect(points);
164 return rotate(reflect(points), 8 - offset);
173 const bool face_orientation,
174 const bool face_flip,
175 const bool face_rotation)
177 static const unsigned int offset[2][2][2] = {
188 mutate_points_with_offset(
190 offset[face_orientation][face_flip][face_rotation]),
194 std::pair<unsigned int, RefinementCase<2>>
195 select_subface_no_and_refinement_case(
196 const unsigned int subface_no,
197 const bool face_orientation,
198 const bool face_flip,
199 const bool face_rotation,
202 constexpr int dim = 3;
272 static const unsigned int
308 equivalent_refine_case[ref_case][subface_no];
309 const unsigned int equ_subface_no =
310 equivalent_subface_number[ref_case][subface_no];
317 (face_orientation == face_rotation ?
318 ref_case_permutation[equ_ref_case] :
321 const unsigned int final_subface_no =
331 return std::make_pair(final_subface_no, final_ref_case);
343 const unsigned int face_no,
347 (void)reference_cell;
349 const unsigned int dim = 1;
353 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
362 const unsigned int face_no,
365 const unsigned int dim = 2;
374 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
393 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
424 const unsigned int face_no,
428 (void)reference_cell;
434 internal::QProjector::project_to_hex_face_and_append(quadrature.
get_points(),
444 const unsigned int face_no,
458 const unsigned int face_no,
459 const bool face_orientation,
460 const bool face_flip,
461 const bool face_rotation)
465 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
466 quadrature, face_orientation, face_flip, face_rotation);
477 const unsigned int face_no,
483 (void)reference_cell;
485 const unsigned int dim = 1;
489 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
498 const unsigned int face_no,
499 const unsigned int subface_no,
503 const unsigned int dim = 2;
514 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
536 quadrature.
point(p)(0) / 2);
540 0.5 + quadrature.
point(p)(0) / 2);
566 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
642 const unsigned int face_no,
643 const unsigned int subface_no,
648 (void)reference_cell;
656 internal::QProjector::project_to_hex_face_and_append(
657 quadrature.
get_points(), face_no, q_points, ref_case, subface_no);
667 const unsigned int face_no,
668 const unsigned int subface_no,
689 const unsigned int face_no,
690 const unsigned int subface_no,
691 const bool face_orientation,
692 const bool face_flip,
693 const bool face_rotation,
698 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
699 quadrature, face_orientation, face_flip, face_rotation);
701 const std::pair<unsigned int, RefinementCase<2>>
702 final_subface_no_and_ref_case =
703 internal::QProjector::select_subface_no_and_refinement_case(
704 subface_no, face_orientation, face_flip, face_rotation, ref_case);
710 final_subface_no_and_ref_case.first,
711 final_subface_no_and_ref_case.second);
723 (void)reference_cell;
725 const unsigned int dim = 1;
730 std::vector<Point<dim>> q_points;
731 q_points.reserve(n_points * n_faces);
732 std::vector<Point<dim>> help(n_points);
737 for (
unsigned int face = 0; face < n_faces; ++face)
739 project_to_face(reference_cell,
740 quadrature[quadrature.
size() == 1 ? 0 : face],
743 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
747 std::vector<double> weights;
748 weights.reserve(n_points * n_faces);
749 for (
unsigned int face = 0; face < n_faces; ++face)
751 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
752 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
753 std::back_inserter(weights));
770 const auto support_points_line =
771 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<2>> {
778 return std::vector<Point<2>>(temp.begin(),
779 temp.begin() + face.first.size());
783 const std::array<std::pair<std::array<Point<2>, 2>,
double>, 3> faces = {
793 std::vector<Point<2>> points;
794 std::vector<double> weights;
797 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
799 for (
unsigned int orientation = 0; orientation < 2; ++orientation)
801 const auto &face = faces[face_no];
805 std::vector<Point<2>> support_points =
806 support_points_line(face, orientation);
809 const auto &sub_quadrature_points =
810 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
811 const auto &sub_quadrature_weights =
812 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
815 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
820 for (
unsigned int i = 0; i < 2; ++i)
823 poly.compute_value(i, sub_quadrature_points[j]);
825 points.emplace_back(mapped_point);
828 weights.emplace_back(sub_quadrature_weights[j] * face.second);
838 const unsigned int dim = 2;
842 unsigned int n_points_total = 0;
844 if (quadrature.
size() == 1)
849 for (
const auto &q : quadrature)
850 n_points_total += q.size();
854 std::vector<Point<dim>> q_points;
855 q_points.reserve(n_points_total);
856 std::vector<Point<dim>> help;
861 for (
unsigned int face = 0; face < n_faces; ++face)
863 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
864 project_to_face(reference_cell,
865 quadrature[quadrature.
size() == 1 ? 0 : face],
868 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
872 std::vector<double> weights;
873 weights.reserve(n_points_total);
874 for (
unsigned int face = 0; face < n_faces; ++face)
876 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
877 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
878 std::back_inserter(weights));
893 const auto process = [&](
const std::vector<std::vector<Point<3>>> &faces) {
895 std::vector<Point<3>> points;
896 std::vector<double> weights;
904 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
909 const unsigned int n_linear_shape_functions = faces[face_no].size();
910 std::vector<Tensor<1, 2>> shape_derivatives;
913 (n_linear_shape_functions == 3 ?
918 for (
unsigned char orientation = 0;
919 orientation < reference_cell.n_face_orientations(face_no);
922 const auto &face = faces[face_no];
924 const boost::container::small_vector<Point<3>, 8> support_points =
925 reference_cell.face_reference_cell(face_no)
926 .permute_by_combined_orientation<
Point<3>>(face, orientation);
929 const auto &sub_quadrature_points =
930 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
931 const auto &sub_quadrature_weights =
932 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
935 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
940 for (
unsigned int i = 0; i < n_linear_shape_functions; ++i)
943 poly.compute_value(i, sub_quadrature_points[j]);
945 points.push_back(mapped_point);
948 const double scaling = [&]() {
949 const unsigned int dim_ = 2;
950 const unsigned int spacedim = 3;
954 shape_derivatives.resize(n_linear_shape_functions);
956 for (
unsigned int i = 0; i < n_linear_shape_functions; ++i)
957 shape_derivatives[i] =
958 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
960 for (
unsigned int k = 0; k < n_linear_shape_functions; ++k)
961 for (
unsigned int i = 0; i < spacedim; ++i)
962 for (
unsigned int j = 0; j < dim_; ++j)
964 shape_derivatives[k][j] * support_points[k][i];
967 for (
unsigned int i = 0; i < dim_; ++i)
968 for (
unsigned int j = 0; j < dim_; ++j)
969 G[i][j] = DX_t[i] * DX_t[j];
974 weights.push_back(sub_quadrature_weights[j] * scaling);
987 std::vector<std::vector<Point<3>>> face_vertex_locations(
988 reference_cell.n_faces());
989 for (
const unsigned int f : reference_cell.face_indices())
991 face_vertex_locations[f].resize(
992 reference_cell.face_reference_cell(f).n_vertices());
993 for (
const unsigned int v :
994 reference_cell.face_reference_cell(f).vertex_indices())
995 face_vertex_locations[f][v] =
996 reference_cell.face_vertex_location<3>(f, v);
999 return process(face_vertex_locations);
1005 const unsigned int dim = 3;
1007 unsigned int n_points_total = 0;
1009 if (quadrature.
size() == 1)
1015 for (
const auto &q : quadrature)
1016 n_points_total += q.size();
1019 n_points_total *= 8;
1022 std::vector<Point<dim>> q_points;
1023 q_points.reserve(n_points_total);
1025 std::vector<double> weights;
1026 weights.reserve(n_points_total);
1028 for (
unsigned int offset = 0; offset < 8; ++offset)
1030 const auto mutation = internal::QProjector::mutate_points_with_offset(
1031 quadrature[0].get_points(), offset);
1033 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
1036 const unsigned int q_index = quadrature.
size() == 1 ? 0 : face;
1038 internal::QProjector::project_to_hex_face_and_append(
1039 q_index > 0 ? internal::QProjector::mutate_points_with_offset(
1040 quadrature[face].get_points(), offset) :
1045 std::copy(quadrature[q_index].get_weights().begin(),
1046 quadrature[q_index].get_weights().end(),
1047 std::back_inserter(weights));
1066 (void)reference_cell;
1068 const unsigned int dim = 1;
1075 std::vector<Point<dim>> q_points;
1076 q_points.reserve(n_points * n_faces * subfaces_per_face);
1077 std::vector<Point<dim>> help(n_points);
1081 for (
unsigned int face = 0; face < n_faces; ++face)
1082 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1084 project_to_subface(reference_cell, quadrature, face, subface, help);
1085 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1089 std::vector<double> weights;
1090 weights.reserve(n_points * n_faces * subfaces_per_face);
1091 for (
unsigned int face = 0; face < n_faces; ++face)
1092 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1095 std::back_inserter(weights));
1097 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1099 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1118 const unsigned int dim = 2;
1120 const unsigned int n_points = quadrature.
size(),
1126 std::vector<Point<dim>> q_points;
1127 q_points.reserve(n_points * n_faces * subfaces_per_face);
1128 std::vector<Point<dim>> help(n_points);
1132 for (
unsigned int face = 0; face < n_faces; ++face)
1133 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1135 project_to_subface(reference_cell, quadrature, face, subface, help);
1136 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1140 std::vector<double> weights;
1141 weights.reserve(n_points * n_faces * subfaces_per_face);
1142 for (
unsigned int face = 0; face < n_faces; ++face)
1143 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1146 std::back_inserter(weights));
1148 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1150 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1169 const unsigned int dim = 3;
1171 const unsigned int n_points = quadrature.
size(),
1173 total_subfaces_per_face = 2 + 2 + 4;
1176 std::vector<Point<dim>> q_points;
1177 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1179 std::vector<double> weights;
1180 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1184 for (
unsigned int offset = 0; offset < 8; ++offset)
1186 const auto mutation =
1187 internal::QProjector::mutate_points_with_offset(quadrature.
get_points(),
1191 for (
unsigned int face = 0; face < n_faces; ++face)
1195 for (
unsigned int subface = 0;
1200 internal::QProjector::project_to_hex_face_and_append(
1210 std::back_inserter(weights));
1214 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1216 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1228 const unsigned int child_no)
1232 (void)reference_cell;
1236 const unsigned int n_q_points = quadrature.
size();
1238 std::vector<Point<dim>> q_points(n_q_points);
1239 for (
unsigned int i = 0; i < n_q_points; ++i)
1247 std::vector<double> weights = quadrature.
get_weights();
1248 for (
unsigned int i = 0; i < n_q_points; ++i)
1263 (void)reference_cell;
1265 const unsigned int n_points = quadrature.
size(),
1268 std::vector<Point<dim>> q_points(n_points * n_children);
1269 std::vector<double> weights(n_points * n_children);
1273 for (
unsigned int child = 0; child < n_children; ++child)
1276 project_to_child(reference_cell, quadrature, child);
1277 for (
unsigned int i = 0; i < n_points; ++i)
1279 q_points[child * n_points + i] = help.
point(i);
1280 weights[child * n_points + i] = help.
weight(i);
1297 (void)reference_cell;
1299 const unsigned int n = quadrature.
size();
1300 std::vector<Point<dim>> points(n);
1301 std::vector<double> weights(n);
1302 const double length = p1.
distance(p2);
1304 for (
unsigned int k = 0; k < n; ++k)
1306 const double alpha = quadrature.
point(k)(0);
1307 points[k] = alpha * p2;
1308 points[k] += (1. - alpha) * p1;
1309 weights[k] = length * quadrature.
weight(k);
1319 const unsigned int face_no,
1320 const bool face_orientation,
1321 const bool face_flip,
1322 const bool face_rotation,
1323 const unsigned int n_quadrature_points)
1329 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1330 n_quadrature_points};
1333 const unsigned int orientation = (face_flip ? 4 : 0) +
1334 (face_rotation ? 2 : 0) +
1335 (face_orientation ? 1 : 0);
1336 return {(6 * face_no + orientation) * n_quadrature_points};
1349 return face_no * n_quadrature_points;
1374 static const unsigned int offset[2][2][2] = {
1393 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1394 n_quadrature_points);
1409 const unsigned int face_no,
1410 const bool face_orientation,
1411 const bool face_flip,
1412 const bool face_rotation,
1420 unsigned int offset = 0;
1423 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1424 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1425 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1426 static const std::array<unsigned int, 5> scale_pyramid = {
1437 if (quadrature.
size() == 1)
1438 offset = scale[0] * quadrature[0].
size() * face_no;
1440 for (
unsigned int i = 0; i < face_no; ++i)
1441 offset += scale[i] * quadrature[i].size();
1446 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1449 const unsigned int orientation = (face_flip ? 4 : 0) +
1450 (face_rotation ? 2 : 0) +
1451 (face_orientation ? 1 : 0);
1455 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1469 if (quadrature.
size() == 1)
1470 return quadrature[0].
size() * face_no;
1473 unsigned int result = 0;
1474 for (
unsigned int i = 0; i < face_no; ++i)
1475 result += quadrature[i].size();
1501 static const unsigned int offset[2][2][2] = {
1512 if (quadrature.
size() == 1)
1514 offset[face_orientation][face_flip][face_rotation] *
1516 quadrature[0].
size();
1519 unsigned int n_points_i = 0;
1520 for (
unsigned int i = 0; i < face_no; ++i)
1521 n_points_i += quadrature[i].size();
1523 unsigned int n_points = 0;
1524 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1526 n_points += quadrature[i].size();
1528 return (n_points_i +
1529 offset[face_orientation][face_flip][face_rotation] *
1546 const unsigned int face_no,
1547 const unsigned int subface_no,
1551 const unsigned int n_quadrature_points,
1555 (void)reference_cell;
1562 n_quadrature_points);
1571 const unsigned int face_no,
1572 const unsigned int subface_no,
1576 const unsigned int n_quadrature_points,
1580 (void)reference_cell;
1587 n_quadrature_points);
1596 const unsigned int face_no,
1597 const unsigned int subface_no,
1598 const bool face_orientation,
1599 const bool face_flip,
1600 const bool face_rotation,
1601 const unsigned int n_quadrature_points,
1604 const unsigned int dim = 3;
1607 (void)reference_cell;
1637 const unsigned int total_subfaces_per_face = 8;
1654 static const unsigned int orientation_offset[2][2][2] = {
1683 static const unsigned int ref_case_offset[3] = {
1689 const std::pair<unsigned int, RefinementCase<2>>
1690 final_subface_no_and_ref_case =
1691 internal::QProjector::select_subface_no_and_refinement_case(
1692 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1694 return (((face_no * total_subfaces_per_face +
1695 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1696 final_subface_no_and_ref_case.first) +
1697 orientation_offset[face_orientation][face_flip][face_rotation]) *
1698 n_quadrature_points);
1707 const unsigned int face_no)
1711 (void)reference_cell;
1713 std::vector<Point<dim>> points(quadrature.
size());
1724 const unsigned int face_no,
1725 const unsigned int subface_no,
1730 (void)reference_cell;
1732 std::vector<Point<dim>> points(quadrature.
size());
1734 reference_cell, quadrature, face_no, subface_no, points, ref_case);
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
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)
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &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_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
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
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const unsigned char orientation) const
CollectionIterator< T > begin() const
unsigned int size() const
CollectionIterator< T > end() const
unsigned int max_n_quadrature_points() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell & get_hypercube()
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
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)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)