35 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]);
48 rotate(
const std::vector<
Point<2>> &points,
const unsigned int n_times)
50 std::vector<Point<2>> q_points;
51 q_points.reserve(points.size());
57 q_points.push_back(p);
62 q_points.emplace_back(1.0 - p[1], p[0]);
67 q_points.emplace_back(1.0 - p[0], 1.0 - p[1]);
72 q_points.emplace_back(p[1], 1.0 - p[0]);
86 project_to_hex_face_and_append(
88 const unsigned int face_no,
91 const unsigned int subface_no = 0)
95 const double const_value = face_no % 2;
100 const unsigned int xi_index = (1 + face_no / 2) % 3,
101 eta_index = (2 + face_no / 2) % 3,
102 const_index = face_no / 2;
109 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
110 eta_translation = 0.0;
119 xi_translation = subface_no % 2 * 0.5;
123 eta_translation = subface_no % 2 * 0.5;
128 xi_translation = int(subface_no % 2) * 0.5;
129 eta_translation = int(subface_no / 2) * 0.5;
141 cell_point[xi_index] = xi_scale * p(0) + xi_translation;
142 cell_point[eta_index] = eta_scale * p(1) + eta_translation;
143 cell_point[const_index] = const_value;
144 q_points.push_back(cell_point);
148 std::vector<Point<2>>
149 mutate_points_with_offset(
const std::vector<
Point<2>> &points,
150 const unsigned int offset)
159 return rotate(points, offset);
161 return reflect(points);
165 return rotate(reflect(points), 8 - offset);
174 const bool face_orientation,
175 const bool face_flip,
176 const bool face_rotation)
178 static const unsigned int offset[2][2][2] = {
189 mutate_points_with_offset(
191 offset[face_orientation][face_flip][face_rotation]),
195 std::pair<unsigned int, RefinementCase<2>>
196 select_subface_no_and_refinement_case(
197 const unsigned int subface_no,
198 const bool face_orientation,
199 const bool face_flip,
200 const bool face_rotation,
203 constexpr
int dim = 3;
273 static const unsigned int
309 equivalent_refine_case[ref_case][subface_no];
310 const unsigned int equ_subface_no =
311 equivalent_subface_number[ref_case][subface_no];
318 (face_orientation == face_rotation ?
319 ref_case_permutation[equ_ref_case] :
322 const unsigned int final_subface_no =
332 return std::make_pair(final_subface_no, final_ref_case);
344 const unsigned int face_no,
350 const unsigned int dim = 1;
354 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
363 const unsigned int face_no,
366 const unsigned int dim = 2;
375 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
394 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
425 const unsigned int face_no,
435 internal::QProjector::project_to_hex_face_and_append(quadrature.
get_points(),
445 const unsigned int face_no,
459 const unsigned int face_no,
460 const bool face_orientation,
461 const bool face_flip,
462 const bool face_rotation)
466 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
467 quadrature, face_orientation, face_flip, face_rotation);
478 const unsigned int face_no,
486 const unsigned int dim = 1;
490 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
499 const unsigned int face_no,
500 const unsigned int subface_no,
504 const unsigned int dim = 2;
515 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
537 quadrature.
point(p)(0) / 2);
541 0.5 + quadrature.
point(p)(0) / 2);
567 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
643 const unsigned int face_no,
644 const unsigned int subface_no,
657 internal::QProjector::project_to_hex_face_and_append(
658 quadrature.
get_points(), face_no, q_points, ref_case, subface_no);
668 const unsigned int face_no,
669 const unsigned int subface_no,
690 const unsigned int face_no,
691 const unsigned int subface_no,
692 const bool face_orientation,
693 const bool face_flip,
694 const bool face_rotation,
699 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
700 quadrature, face_orientation, face_flip, face_rotation);
702 const std::pair<unsigned int, RefinementCase<2>>
703 final_subface_no_and_ref_case =
704 internal::QProjector::select_subface_no_and_refinement_case(
705 subface_no, face_orientation, face_flip, face_rotation, ref_case);
711 final_subface_no_and_ref_case.first,
712 final_subface_no_and_ref_case.second);
726 const unsigned int dim = 1;
731 std::vector<Point<dim>> q_points;
732 q_points.reserve(n_points * n_faces);
733 std::vector<Point<dim>> help(n_points);
738 for (
unsigned int face = 0; face < n_faces; ++face)
741 quadrature[quadrature.
size() == 1 ? 0 : face],
744 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
748 std::vector<double> weights;
749 weights.reserve(n_points * n_faces);
750 for (
unsigned int face = 0; face < n_faces; ++face)
752 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
753 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
754 std::back_inserter(weights));
771 const auto support_points_line =
772 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<2>> {
779 return std::vector<Point<2>>(temp.begin(),
780 temp.begin() + face.first.size());
784 const std::array<std::pair<std::array<Point<2>, 2>,
double>, 3> faces = {
794 std::vector<Point<2>> points;
795 std::vector<double> weights;
798 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
800 for (
unsigned int orientation = 0; orientation < 2; ++orientation)
802 const auto &face = faces[face_no];
806 std::vector<Point<2>> support_points =
807 support_points_line(face, orientation);
810 const auto &sub_quadrature_points =
811 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
812 const auto &sub_quadrature_weights =
813 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
816 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
821 for (
unsigned int i = 0; i < 2; ++i)
824 poly.compute_value(i, sub_quadrature_points[j]);
826 points.emplace_back(mapped_point);
829 weights.emplace_back(sub_quadrature_weights[j] * face.second);
839 const unsigned int dim = 2;
843 unsigned int n_points_total = 0;
845 if (quadrature.
size() == 1)
850 for (
const auto &q : quadrature)
851 n_points_total += q.size();
855 std::vector<Point<dim>> q_points;
856 q_points.reserve(n_points_total);
857 std::vector<Point<dim>> help;
862 for (
unsigned int face = 0; face < n_faces; ++face)
864 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
866 quadrature[quadrature.
size() == 1 ? 0 : face],
869 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
873 std::vector<double> weights;
874 weights.reserve(n_points_total);
875 for (
unsigned int face = 0; face < n_faces; ++face)
877 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
878 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
879 std::back_inserter(weights));
894 const auto process = [&](
const std::vector<std::vector<Point<3>>> &faces) {
896 std::vector<Point<3>> points;
897 std::vector<double> weights;
905 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
910 const unsigned int n_linear_shape_functions = faces[face_no].size();
911 std::vector<Tensor<1, 2>> shape_derivatives;
914 (n_linear_shape_functions == 3 ?
919 for (
unsigned char orientation = 0;
923 const auto &face = faces[face_no];
925 const boost::container::small_vector<Point<3>, 8> support_points =
927 .permute_by_combined_orientation<
Point<3>>(face, orientation);
930 const auto &sub_quadrature_points =
931 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
932 const auto &sub_quadrature_weights =
933 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
936 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
941 for (
unsigned int i = 0; i < n_linear_shape_functions; ++i)
944 poly.compute_value(i, sub_quadrature_points[j]);
946 points.push_back(mapped_point);
949 const double scaling = [&]() {
950 const unsigned int dim_ = 2;
951 const unsigned int spacedim = 3;
955 shape_derivatives.resize(n_linear_shape_functions);
957 for (
unsigned int i = 0; i < n_linear_shape_functions; ++i)
958 shape_derivatives[i] =
959 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
961 for (
unsigned int k = 0; k < n_linear_shape_functions; ++k)
962 for (
unsigned int i = 0; i < spacedim; ++i)
963 for (
unsigned int j = 0; j < dim_; ++j)
965 shape_derivatives[k][j] * support_points[k][i];
968 for (
unsigned int i = 0; i < dim_; ++i)
969 for (
unsigned int j = 0; j < dim_; ++j)
970 G[i][j] = DX_t[i] * DX_t[j];
975 weights.push_back(sub_quadrature_weights[j] * scaling);
988 std::vector<std::vector<Point<3>>> face_vertex_locations(
992 face_vertex_locations[f].resize(
994 for (
const unsigned int v :
996 face_vertex_locations[f][v] =
1000 return process(face_vertex_locations);
1006 const unsigned int dim = 3;
1008 unsigned int n_points_total = 0;
1010 if (quadrature.
size() == 1)
1016 for (
const auto &q : quadrature)
1017 n_points_total += q.size();
1020 n_points_total *= 8;
1023 std::vector<Point<dim>> q_points;
1024 q_points.reserve(n_points_total);
1026 std::vector<double> weights;
1027 weights.reserve(n_points_total);
1029 for (
unsigned int offset = 0; offset < 8; ++offset)
1031 const auto mutation = internal::QProjector::mutate_points_with_offset(
1032 quadrature[0].get_points(), offset);
1034 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
1037 const unsigned int q_index = quadrature.
size() == 1 ? 0 : face;
1039 internal::QProjector::project_to_hex_face_and_append(
1040 q_index > 0 ? internal::QProjector::mutate_points_with_offset(
1041 quadrature[face].get_points(), offset) :
1047 quadrature[q_index].get_weights().
end(),
1048 std::back_inserter(weights));
1069 const unsigned int dim = 1;
1076 std::vector<Point<dim>> q_points;
1077 q_points.reserve(n_points * n_faces * subfaces_per_face);
1078 std::vector<Point<dim>> help(n_points);
1082 for (
unsigned int face = 0; face < n_faces; ++face)
1083 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1085 project_to_subface(
reference_cell, quadrature, face, subface, help);
1086 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1090 std::vector<double> weights;
1091 weights.reserve(n_points * n_faces * subfaces_per_face);
1092 for (
unsigned int face = 0; face < n_faces; ++face)
1093 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1096 std::back_inserter(weights));
1098 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1100 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1119 const unsigned int dim = 2;
1121 const unsigned int n_points = quadrature.
size(),
1127 std::vector<Point<dim>> q_points;
1128 q_points.reserve(n_points * n_faces * subfaces_per_face);
1129 std::vector<Point<dim>> help(n_points);
1133 for (
unsigned int face = 0; face < n_faces; ++face)
1134 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1136 project_to_subface(
reference_cell, quadrature, face, subface, help);
1137 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1141 std::vector<double> weights;
1142 weights.reserve(n_points * n_faces * subfaces_per_face);
1143 for (
unsigned int face = 0; face < n_faces; ++face)
1144 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1147 std::back_inserter(weights));
1149 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1151 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1170 const unsigned int dim = 3;
1172 const unsigned int n_points = quadrature.
size(),
1174 total_subfaces_per_face = 2 + 2 + 4;
1177 std::vector<Point<dim>> q_points;
1178 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1180 std::vector<double> weights;
1181 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1185 for (
unsigned int offset = 0; offset < 8; ++offset)
1187 const auto mutation =
1188 internal::QProjector::mutate_points_with_offset(quadrature.
get_points(),
1192 for (
unsigned int face = 0; face < n_faces; ++face)
1196 for (
unsigned int subface = 0;
1201 internal::QProjector::project_to_hex_face_and_append(
1211 std::back_inserter(weights));
1215 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1217 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1229 const unsigned int child_no)
1237 const unsigned int n_q_points = quadrature.
size();
1239 std::vector<Point<dim>> q_points(n_q_points);
1240 for (
unsigned int i = 0; i < n_q_points; ++i)
1248 std::vector<double> weights = quadrature.
get_weights();
1249 for (
unsigned int i = 0; i < n_q_points; ++i)
1266 const unsigned int n_points = quadrature.
size(),
1269 std::vector<Point<dim>> q_points(n_points * n_children);
1270 std::vector<double> weights(n_points * n_children);
1274 for (
unsigned int child = 0; child < n_children; ++child)
1278 for (
unsigned int i = 0; i < n_points; ++i)
1280 q_points[child * n_points + i] = help.
point(i);
1281 weights[child * n_points + i] = help.
weight(i);
1300 const unsigned int n = quadrature.
size();
1301 std::vector<Point<dim>> points(n);
1302 std::vector<double> weights(n);
1303 const double length = p1.
distance(p2);
1305 for (
unsigned int k = 0; k < n; ++k)
1307 const double alpha = quadrature.
point(k)(0);
1308 points[k] = alpha * p2;
1309 points[k] += (1. - alpha) * p1;
1310 weights[k] = length * quadrature.
weight(k);
1320 const unsigned int face_no,
1321 const bool face_orientation,
1322 const bool face_flip,
1323 const bool face_rotation,
1324 const unsigned int n_quadrature_points)
1330 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1331 n_quadrature_points};
1334 const unsigned char orientation =
1341 n_quadrature_points};
1354 return face_no * n_quadrature_points;
1379 static const unsigned int offset[2][2][2] = {
1398 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1399 n_quadrature_points);
1414 const unsigned int face_no,
1415 const bool face_orientation,
1416 const bool face_flip,
1417 const bool face_rotation,
1425 unsigned int offset = 0;
1428 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1429 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1430 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1431 static const std::array<unsigned int, 5> scale_pyramid = {
1442 if (quadrature.
size() == 1)
1443 offset =
scale[0] * quadrature[0].size() * face_no;
1445 for (
unsigned int i = 0; i < face_no; ++i)
1446 offset +=
scale[i] * quadrature[i].size();
1451 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1458 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1472 if (quadrature.
size() == 1)
1473 return quadrature[0].
size() * face_no;
1476 unsigned int result = 0;
1477 for (
unsigned int i = 0; i < face_no; ++i)
1478 result += quadrature[i].size();
1504 static const unsigned int offset[2][2][2] = {
1515 if (quadrature.
size() == 1)
1517 offset[face_orientation][face_flip][face_rotation] *
1519 quadrature[0].
size();
1522 unsigned int n_points_i = 0;
1523 for (
unsigned int i = 0; i < face_no; ++i)
1524 n_points_i += quadrature[i].size();
1526 unsigned int n_points = 0;
1527 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1529 n_points += quadrature[i].size();
1531 return (n_points_i +
1532 offset[face_orientation][face_flip][face_rotation] *
1549 const unsigned int face_no,
1550 const unsigned int subface_no,
1554 const unsigned int n_quadrature_points,
1565 n_quadrature_points);
1574 const unsigned int face_no,
1575 const unsigned int subface_no,
1579 const unsigned int n_quadrature_points,
1590 n_quadrature_points);
1599 const unsigned int face_no,
1600 const unsigned int subface_no,
1601 const bool face_orientation,
1602 const bool face_flip,
1603 const bool face_rotation,
1604 const unsigned int n_quadrature_points,
1607 const unsigned int dim = 3;
1640 const unsigned int total_subfaces_per_face = 8;
1657 static const unsigned int orientation_offset[2][2][2] = {
1686 static const unsigned int ref_case_offset[3] = {
1692 const std::pair<unsigned int, RefinementCase<2>>
1693 final_subface_no_and_ref_case =
1694 internal::QProjector::select_subface_no_and_refinement_case(
1695 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1697 return (((face_no * total_subfaces_per_face +
1698 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1699 final_subface_no_and_ref_case.first) +
1700 orientation_offset[face_orientation][face_flip][face_rotation]) *
1701 n_quadrature_points);
1710 const unsigned int face_no)
1716 std::vector<Point<dim>> points(quadrature.
size());
1727 const unsigned int face_no,
1728 const unsigned int subface_no,
1735 std::vector<Point<dim>> points(quadrature.
size());
1737 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 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_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 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_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 std::vector< Point< dim > > & get_points() const
const std::vector< double > & get_weights() const
double weight(const unsigned int i) const
const Point< dim > & point(const unsigned int i) 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
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
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 & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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 Line
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
void copy(const T *begin, const T *end, U *dest)
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
static const unsigned int invalid_unsigned_int
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)