35 std::vector<Point<2>> q_points(q.
get_points());
46 std::vector<Point<2>> q_points(q.
size());
47 for (
unsigned int i = 0; i < q.
size(); ++i)
53 q_points[i] = q.
point(i);
58 q_points[i][0] = 1.0 - q.
point(i)[1];
59 q_points[i][1] = q.
point(i)[0];
63 q_points[i][0] = 1.0 - q.
point(i)[0];
64 q_points[i][1] = 1.0 - q.
point(i)[1];
68 q_points[i][0] = q.
point(i)[1];
69 q_points[i][1] = 1.0 - q.
point(i)[0];
86 const unsigned int face_no,
92 const unsigned int dim = 1;
96 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
105 const unsigned int face_no,
108 const unsigned int dim = 2;
117 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
136 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
167 const unsigned int face_no,
173 const unsigned int dim = 3;
178 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
217 const unsigned int face_no,
225 const unsigned int dim = 1;
229 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
238 const unsigned int face_no,
239 const unsigned int subface_no,
243 const unsigned int dim = 2;
254 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
276 quadrature.
point(p)(0) / 2);
280 0.5 + quadrature.
point(p)(0) / 2);
306 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
382 const unsigned int face_no,
383 const unsigned int subface_no,
390 const unsigned int dim = 3;
399 double const_value = face_no % 2;
408 const_index = face_no / 2;
414 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
415 eta_translation = 0.0;
439 xi_translation = subface_no % 2 * 0.5;
443 eta_translation = subface_no % 2 * 0.5;
448 xi_translation = int(subface_no % 2) * 0.5;
449 eta_translation = int(subface_no / 2) * 0.5;
457 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
459 q_points[p][xi_index] =
460 xi_scale * quadrature.
point(p)(0) + xi_translation;
461 q_points[p][eta_index] =
462 eta_scale * quadrature.
point(p)(1) + eta_translation;
463 q_points[p][const_index] = const_value;
477 const unsigned int dim = 1;
482 std::vector<Point<dim>> q_points;
483 q_points.reserve(n_points * n_faces);
484 std::vector<Point<dim>> help(n_points);
489 for (
unsigned int face = 0; face < n_faces; ++face)
492 quadrature[quadrature.
size() == 1 ? 0 : face],
495 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
499 std::vector<double> weights;
500 weights.reserve(n_points * n_faces);
501 for (
unsigned int face = 0; face < n_faces; ++face)
503 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
504 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
505 std::back_inserter(weights));
522 const auto support_points_line =
523 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<2>> {
525 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
529 return std::vector<Point<2>>(temp.begin(),
530 temp.begin() + face.first.size());
534 const std::array<std::pair<std::array<Point<2>, 2>,
double>, 3> faces = {
544 std::vector<Point<2>> points;
545 std::vector<double> weights;
548 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
550 for (
unsigned int orientation = 0; orientation < 2; ++orientation)
552 const auto &face = faces[face_no];
556 std::vector<Point<2>> support_points =
557 support_points_line(face, orientation);
560 const auto &sub_quadrature_points =
561 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
562 const auto &sub_quadrature_weights =
563 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
566 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
571 for (
unsigned int i = 0; i < 2; ++i)
574 poly.compute_value(i, sub_quadrature_points[j]);
576 points.emplace_back(mapped_point);
579 weights.emplace_back(sub_quadrature_weights[j] * face.second);
584 return {points, weights};
589 const unsigned int dim = 2;
593 unsigned int n_points_total = 0;
595 if (quadrature.
size() == 1)
600 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
601 n_points_total += quadrature[i].size();
605 std::vector<Point<dim>> q_points;
606 q_points.reserve(n_points_total);
607 std::vector<Point<dim>> help;
612 for (
unsigned int face = 0; face < n_faces; ++face)
614 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
616 quadrature[quadrature.
size() == 1 ? 0 : face],
619 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
623 std::vector<double> weights;
624 weights.reserve(n_points_total);
625 for (
unsigned int face = 0; face < n_faces; ++face)
627 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
628 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
629 std::back_inserter(weights));
644 const auto support_points_tri =
645 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
647 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
651 return std::vector<Point<3>>(temp.begin(),
652 temp.begin() + face.first.size());
655 const auto support_points_quad =
656 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
658 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
662 return std::vector<Point<3>>(temp.begin(),
663 temp.begin() + face.first.size());
666 const auto process = [&](
const auto &faces) {
668 std::vector<Point<3>> points;
669 std::vector<double> weights;
677 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
681 const unsigned int n_shape_functions = faces[face_no].first.size();
684 n_shape_functions == 3 ?
689 for (
unsigned int orientation = 0;
690 orientation < (n_shape_functions * 2);
693 const auto &face = faces[face_no];
695 const auto support_points =
696 n_shape_functions == 3 ? support_points_tri(face, orientation) :
697 support_points_quad(face, orientation);
700 const auto &sub_quadrature_points =
701 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
702 const auto &sub_quadrature_weights =
703 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
706 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
711 for (
unsigned int i = 0; i < n_shape_functions; ++i)
714 poly.compute_value(i, sub_quadrature_points[j]);
716 points.push_back(mapped_point);
719 const double scaling = [&]() {
720 const auto & supp_pts = support_points;
721 const unsigned int dim_ = 2;
722 const unsigned int spacedim = 3;
724 double result[spacedim][dim_];
726 std::vector<Tensor<1, dim_>> shape_derivatives(
729 for (
unsigned int i = 0; i < n_shape_functions; ++i)
730 shape_derivatives[i] =
731 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
733 for (
unsigned int i = 0; i < spacedim; ++i)
734 for (
unsigned int j = 0; j < dim_; ++j)
735 result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
736 for (
unsigned int k = 1; k < n_shape_functions; ++k)
737 for (
unsigned int i = 0; i < spacedim; ++i)
738 for (
unsigned int j = 0; j < dim_; ++j)
740 shape_derivatives[k][j] * supp_pts[k][i];
744 for (
unsigned int i = 0; i < spacedim; ++i)
745 for (
unsigned int j = 0; j < dim_; ++j)
746 contravariant[i][j] = result[i][j];
750 for (
unsigned int i = 0; i < spacedim; ++i)
751 for (
unsigned int j = 0; j < dim_; ++j)
752 DX_t[j][i] = contravariant[i][j];
755 for (
unsigned int i = 0; i < dim_; ++i)
756 for (
unsigned int j = 0; j < dim_; ++j)
757 G[i][j] = DX_t[i] * DX_t[j];
762 weights.push_back(sub_quadrature_weights[j] * scaling);
775 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
793 return process(faces);
797 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
822 return process(faces);
826 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
849 return process(faces);
855 const unsigned int dim = 3;
857 unsigned int n_points_total = 0;
859 if (quadrature.
size() == 1)
864 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
865 n_points_total += quadrature[i].size();
871 std::vector<Point<dim>> q_points;
872 q_points.reserve(n_points_total);
873 std::vector<Point<dim>> help;
876 std::vector<double> weights;
877 weights.reserve(n_points_total);
883 for (
unsigned int i = 0; i < 8; ++i)
887 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
892 const auto quadrature_f =
893 quadrature[quadrature.
size() == 1 ? 0 : face];
897 mutation = quadrature_f;
909 mutation = internal::QProjector::reflect(quadrature_f);
913 internal::QProjector::reflect(quadrature_f), 3);
917 internal::QProjector::reflect(quadrature_f), 2);
921 internal::QProjector::reflect(quadrature_f), 1);
927 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
929 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
933 std::back_inserter(weights));
954 const unsigned int dim = 1;
961 std::vector<Point<dim>> q_points;
962 q_points.reserve(n_points * n_faces * subfaces_per_face);
963 std::vector<Point<dim>> help(n_points);
967 for (
unsigned int face = 0; face < n_faces; ++face)
968 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
970 project_to_subface(
reference_cell, quadrature, face, subface, help);
971 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
975 std::vector<double> weights;
976 weights.reserve(n_points * n_faces * subfaces_per_face);
977 for (
unsigned int face = 0; face < n_faces; ++face)
978 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
981 std::back_inserter(weights));
983 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
985 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1004 const unsigned int dim = 2;
1006 const unsigned int n_points = quadrature.
size(),
1012 std::vector<Point<dim>> q_points;
1013 q_points.reserve(n_points * n_faces * subfaces_per_face);
1014 std::vector<Point<dim>> help(n_points);
1018 for (
unsigned int face = 0; face < n_faces; ++face)
1019 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1021 project_to_subface(
reference_cell, quadrature, face, subface, help);
1022 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1026 std::vector<double> weights;
1027 weights.reserve(n_points * n_faces * subfaces_per_face);
1028 for (
unsigned int face = 0; face < n_faces; ++face)
1029 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1032 std::back_inserter(weights));
1034 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1036 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1055 const unsigned int dim = 3;
1056 SubQuadrature q_reflected = internal::QProjector::reflect(quadrature);
1066 const unsigned int n_points = quadrature.
size(),
1068 total_subfaces_per_face = 2 + 2 + 4;
1071 std::vector<Point<dim>> q_points;
1072 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1073 std::vector<Point<dim>> help(n_points);
1075 std::vector<double> weights;
1076 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1082 for (
const auto &mutation : q)
1086 for (
unsigned int face = 0; face < n_faces; ++face)
1090 for (
unsigned int subface = 0;
1101 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1105 for (
unsigned int face = 0; face < n_faces; ++face)
1109 for (
unsigned int subface = 0;
1113 std::copy(mutation.get_weights().begin(),
1114 mutation.get_weights().end(),
1115 std::back_inserter(weights));
1118 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1120 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1132 const unsigned int child_no)
1140 const unsigned int n_q_points = quadrature.
size();
1142 std::vector<Point<dim>> q_points(n_q_points);
1143 for (
unsigned int i = 0; i < n_q_points; ++i)
1151 std::vector<double> weights = quadrature.
get_weights();
1152 for (
unsigned int i = 0; i < n_q_points; ++i)
1169 const unsigned int n_points = quadrature.
size(),
1172 std::vector<Point<dim>> q_points(n_points * n_children);
1173 std::vector<double> weights(n_points * n_children);
1177 for (
unsigned int child = 0; child < n_children; ++child)
1181 for (
unsigned int i = 0; i < n_points; ++i)
1183 q_points[child * n_points + i] = help.
point(i);
1184 weights[child * n_points + i] = help.
weight(i);
1203 const unsigned int n = quadrature.
size();
1204 std::vector<Point<dim>> points(n);
1205 std::vector<double> weights(n);
1206 const double length = p1.
distance(p2);
1208 for (
unsigned int k = 0; k < n; ++k)
1210 const double alpha = quadrature.
point(k)(0);
1211 points[k] = alpha * p2;
1212 points[k] += (1. - alpha) * p1;
1213 weights[k] = length * quadrature.
weight(k);
1223 const unsigned int face_no,
1224 const bool face_orientation,
1225 const bool face_flip,
1226 const bool face_rotation,
1227 const unsigned int n_quadrature_points)
1233 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1234 n_quadrature_points};
1237 const unsigned int orientation = (face_flip ? 4 : 0) +
1238 (face_rotation ? 2 : 0) +
1239 (face_orientation ? 1 : 0);
1240 return {(6 * face_no + orientation) * n_quadrature_points};
1253 return face_no * n_quadrature_points;
1278 static const unsigned int offset[2][2][2] = {
1297 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1298 n_quadrature_points);
1313 const unsigned int face_no,
1314 const bool face_orientation,
1315 const bool face_flip,
1316 const bool face_rotation,
1324 unsigned int offset = 0;
1327 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1328 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1329 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1330 static const std::array<unsigned int, 5> scale_pyramid = {
1341 if (quadrature.
size() == 1)
1342 offset =
scale[0] * quadrature[0].size() * face_no;
1344 for (
unsigned int i = 0; i < face_no; ++i)
1345 offset +=
scale[i] * quadrature[i].size();
1350 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1353 const unsigned int orientation = (face_flip ? 4 : 0) +
1354 (face_rotation ? 2 : 0) +
1355 (face_orientation ? 1 : 0);
1359 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1373 if (quadrature.
size() == 1)
1374 return quadrature[0].
size() * face_no;
1377 unsigned int result = 0;
1378 for (
unsigned int i = 0; i < face_no; ++i)
1379 result += quadrature[i].size();
1405 static const unsigned int offset[2][2][2] = {
1416 if (quadrature.
size() == 1)
1418 offset[face_orientation][face_flip][face_rotation] *
1420 quadrature[0].
size();
1423 unsigned int n_points_i = 0;
1424 for (
unsigned int i = 0; i < face_no; ++i)
1425 n_points_i += quadrature[i].size();
1427 unsigned int n_points = 0;
1428 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1430 n_points += quadrature[i].size();
1432 return (n_points_i +
1433 offset[face_orientation][face_flip][face_rotation] *
1450 const unsigned int face_no,
1451 const unsigned int subface_no,
1455 const unsigned int n_quadrature_points,
1466 n_quadrature_points);
1475 const unsigned int face_no,
1476 const unsigned int subface_no,
1480 const unsigned int n_quadrature_points,
1491 n_quadrature_points);
1500 const unsigned int face_no,
1501 const unsigned int subface_no,
1502 const bool face_orientation,
1503 const bool face_flip,
1504 const bool face_rotation,
1505 const unsigned int n_quadrature_points,
1508 const unsigned int dim = 3;
1541 const unsigned int total_subfaces_per_face = 8;
1558 static const unsigned int orientation_offset[2][2][2] = {
1587 static const unsigned int ref_case_offset[3] = {
1663 static const unsigned int
1699 equivalent_refine_case[ref_case][subface_no];
1700 const unsigned int equ_subface_no =
1701 equivalent_subface_number[ref_case][subface_no];
1708 (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1737 const unsigned int final_subface_no =
1746 return (((face_no * total_subfaces_per_face +
1747 ref_case_offset[final_ref_case - 1] + final_subface_no) +
1748 orientation_offset[face_orientation][face_flip][face_rotation]) *
1749 n_quadrature_points);
1758 const unsigned int face_no)
1764 std::vector<Point<dim>> points(quadrature.
size());
1775 const unsigned int face_no,
1776 const unsigned int subface_no,
1783 std::vector<Point<dim>> points(quadrature.
size());
1785 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_line(const ReferenceCell reference_cell, const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
static Quadrature< dim > project_to_child(const ReferenceCell reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
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_subfaces(const ReferenceCell reference_cell, const SubQuadrature &quadrature)
static void project_to_face(const ReferenceCell reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
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)
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
std::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
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)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
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
void copy(const T *begin, const T *end, U *dest)
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)