37 std::vector<Point<2>> q_points(q.
get_points());
48 std::vector<Point<2>> q_points(q.
size());
49 for (
unsigned int i = 0; i < q.
size(); ++i)
55 q_points[i] = q.
point(i);
60 q_points[i][0] = 1.0 - q.
point(i)[1];
61 q_points[i][1] = q.
point(i)[0];
65 q_points[i][0] = 1.0 - q.
point(i)[0];
66 q_points[i][1] = 1.0 - q.
point(i)[1];
70 q_points[i][0] = q.
point(i)[1];
71 q_points[i][1] = 1.0 - q.
point(i)[0];
88 const unsigned int face_no,
94 const unsigned int dim = 1;
98 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
107 const unsigned int face_no,
110 const unsigned int dim = 2;
119 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
138 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
169 const unsigned int face_no,
175 const unsigned int dim = 3;
180 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
219 const unsigned int face_no,
227 const unsigned int dim = 1;
231 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
240 const unsigned int face_no,
241 const unsigned int subface_no,
245 const unsigned int dim = 2;
256 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
278 quadrature.
point(p)(0) / 2);
282 0.5 + quadrature.
point(p)(0) / 2);
308 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
384 const unsigned int face_no,
385 const unsigned int subface_no,
392 const unsigned int dim = 3;
401 double const_value = face_no % 2;
410 const_index = face_no / 2;
416 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
417 eta_translation = 0.0;
441 xi_translation = subface_no % 2 * 0.5;
445 eta_translation = subface_no % 2 * 0.5;
450 xi_translation = int(subface_no % 2) * 0.5;
451 eta_translation = int(subface_no / 2) * 0.5;
459 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
461 q_points[p][xi_index] =
462 xi_scale * quadrature.
point(p)(0) + xi_translation;
463 q_points[p][eta_index] =
464 eta_scale * quadrature.
point(p)(1) + eta_translation;
465 q_points[p][const_index] = const_value;
479 const unsigned int dim = 1;
484 std::vector<Point<dim>> q_points;
485 q_points.reserve(n_points * n_faces);
486 std::vector<Point<dim>> help(n_points);
491 for (
unsigned int face = 0; face < n_faces; ++face)
494 quadrature[quadrature.
size() == 1 ? 0 : face],
497 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
501 std::vector<double> weights;
502 weights.reserve(n_points * n_faces);
503 for (
unsigned int face = 0; face < n_faces; ++face)
505 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
506 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
507 std::back_inserter(weights));
524 const auto support_points_line =
525 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<2>> {
527 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
531 return std::vector<Point<2>>(temp.begin(),
532 temp.begin() + face.first.size());
536 const std::array<std::pair<std::array<Point<2>, 2>,
double>, 3> faces = {
546 std::vector<Point<2>> points;
547 std::vector<double> weights;
550 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
552 for (
unsigned int orientation = 0; orientation < 2; ++orientation)
554 const auto &face = faces[face_no];
558 std::vector<Point<2>> support_points =
559 support_points_line(face, orientation);
562 const auto &sub_quadrature_points =
563 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
564 const auto &sub_quadrature_weights =
565 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
568 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
573 for (
unsigned int i = 0; i < 2; ++i)
576 poly.compute_value(i, sub_quadrature_points[j]);
578 points.emplace_back(mapped_point);
581 weights.emplace_back(sub_quadrature_weights[j] * face.second);
586 return {points, weights};
591 const unsigned int dim = 2;
595 unsigned int n_points_total = 0;
597 if (quadrature.
size() == 1)
602 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
603 n_points_total += quadrature[i].size();
607 std::vector<Point<dim>> q_points;
608 q_points.reserve(n_points_total);
609 std::vector<Point<dim>> help;
614 for (
unsigned int face = 0; face < n_faces; ++face)
616 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
618 quadrature[quadrature.
size() == 1 ? 0 : face],
621 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
625 std::vector<double> weights;
626 weights.reserve(n_points_total);
627 for (
unsigned int face = 0; face < n_faces; ++face)
629 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
begin(),
630 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().
end(),
631 std::back_inserter(weights));
646 const auto support_points_tri =
647 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
649 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
653 return std::vector<Point<3>>(temp.begin(),
654 temp.begin() + face.first.size());
657 const auto support_points_quad =
658 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
660 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
664 return std::vector<Point<3>>(temp.begin(),
665 temp.begin() + face.first.size());
668 const auto process = [&](
const auto &faces) {
670 std::vector<Point<3>> points;
671 std::vector<double> weights;
679 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
683 const unsigned int n_shape_functions = faces[face_no].first.size();
686 n_shape_functions == 3 ?
691 for (
unsigned int orientation = 0;
692 orientation < (n_shape_functions * 2);
695 const auto &face = faces[face_no];
697 const auto support_points =
698 n_shape_functions == 3 ? support_points_tri(face, orientation) :
699 support_points_quad(face, orientation);
702 const auto &sub_quadrature_points =
703 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
704 const auto &sub_quadrature_weights =
705 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
708 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
713 for (
unsigned int i = 0; i < n_shape_functions; ++i)
716 poly.compute_value(i, sub_quadrature_points[j]);
718 points.push_back(mapped_point);
721 const double scaling = [&]() {
722 const auto & supp_pts = support_points;
723 const unsigned int dim_ = 2;
724 const unsigned int spacedim = 3;
726 double result[spacedim][dim_];
728 std::vector<Tensor<1, dim_>> shape_derivatives(
731 for (
unsigned int i = 0; i < n_shape_functions; ++i)
732 shape_derivatives[i] =
733 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
735 for (
unsigned int i = 0; i < spacedim; ++i)
736 for (
unsigned int j = 0; j < dim_; ++j)
737 result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
738 for (
unsigned int k = 1; k < n_shape_functions; ++k)
739 for (
unsigned int i = 0; i < spacedim; ++i)
740 for (
unsigned int j = 0; j < dim_; ++j)
742 shape_derivatives[k][j] * supp_pts[k][i];
746 for (
unsigned int i = 0; i < spacedim; ++i)
747 for (
unsigned int j = 0; j < dim_; ++j)
748 contravariant[i][j] = result[i][j];
752 for (
unsigned int i = 0; i < spacedim; ++i)
753 for (
unsigned int j = 0; j < dim_; ++j)
754 DX_t[j][i] = contravariant[i][j];
757 for (
unsigned int i = 0; i < dim_; ++i)
758 for (
unsigned int j = 0; j < dim_; ++j)
759 G[i][j] = DX_t[i] * DX_t[j];
764 weights.push_back(sub_quadrature_weights[j] * scaling);
777 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
795 return process(faces);
799 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
824 return process(faces);
828 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
851 return process(faces);
857 const unsigned int dim = 3;
859 unsigned int n_points_total = 0;
861 if (quadrature.
size() == 1)
866 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
867 n_points_total += quadrature[i].size();
873 std::vector<Point<dim>> q_points;
874 q_points.reserve(n_points_total);
875 std::vector<Point<dim>> help;
878 std::vector<double> weights;
879 weights.reserve(n_points_total);
885 for (
unsigned int i = 0; i < 8; ++i)
889 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
894 const auto &quadrature_f =
895 quadrature[quadrature.
size() == 1 ? 0 : face];
899 mutation = quadrature_f;
911 mutation = internal::QProjector::reflect(quadrature_f);
915 internal::QProjector::reflect(quadrature_f), 3);
919 internal::QProjector::reflect(quadrature_f), 2);
923 internal::QProjector::reflect(quadrature_f), 1);
929 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
931 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
935 std::back_inserter(weights));
956 const unsigned int dim = 1;
963 std::vector<Point<dim>> q_points;
964 q_points.reserve(n_points * n_faces * subfaces_per_face);
965 std::vector<Point<dim>> help(n_points);
969 for (
unsigned int face = 0; face < n_faces; ++face)
970 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
972 project_to_subface(
reference_cell, quadrature, face, subface, help);
973 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
977 std::vector<double> weights;
978 weights.reserve(n_points * n_faces * subfaces_per_face);
979 for (
unsigned int face = 0; face < n_faces; ++face)
980 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
983 std::back_inserter(weights));
985 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
987 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1006 const unsigned int dim = 2;
1008 const unsigned int n_points = quadrature.
size(),
1014 std::vector<Point<dim>> q_points;
1015 q_points.reserve(n_points * n_faces * subfaces_per_face);
1016 std::vector<Point<dim>> help(n_points);
1020 for (
unsigned int face = 0; face < n_faces; ++face)
1021 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1023 project_to_subface(
reference_cell, quadrature, face, subface, help);
1024 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1028 std::vector<double> weights;
1029 weights.reserve(n_points * n_faces * subfaces_per_face);
1030 for (
unsigned int face = 0; face < n_faces; ++face)
1031 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1034 std::back_inserter(weights));
1036 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1038 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1057 const unsigned int dim = 3;
1058 SubQuadrature q_reflected = internal::QProjector::reflect(quadrature);
1068 const unsigned int n_points = quadrature.
size(),
1070 total_subfaces_per_face = 2 + 2 + 4;
1073 std::vector<Point<dim>> q_points;
1074 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1075 std::vector<Point<dim>> help(n_points);
1077 std::vector<double> weights;
1078 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1084 for (
const auto &mutation : q)
1088 for (
unsigned int face = 0; face < n_faces; ++face)
1092 for (
unsigned int subface = 0;
1103 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1107 for (
unsigned int face = 0; face < n_faces; ++face)
1111 for (
unsigned int subface = 0;
1115 std::copy(mutation.get_weights().begin(),
1116 mutation.get_weights().end(),
1117 std::back_inserter(weights));
1120 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1122 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1134 const unsigned int child_no)
1142 const unsigned int n_q_points = quadrature.
size();
1144 std::vector<Point<dim>> q_points(n_q_points);
1145 for (
unsigned int i = 0; i < n_q_points; ++i)
1153 std::vector<double> weights = quadrature.
get_weights();
1154 for (
unsigned int i = 0; i < n_q_points; ++i)
1171 const unsigned int n_points = quadrature.
size(),
1174 std::vector<Point<dim>> q_points(n_points * n_children);
1175 std::vector<double> weights(n_points * n_children);
1179 for (
unsigned int child = 0; child < n_children; ++child)
1183 for (
unsigned int i = 0; i < n_points; ++i)
1185 q_points[child * n_points + i] = help.
point(i);
1186 weights[child * n_points + i] = help.
weight(i);
1205 const unsigned int n = quadrature.
size();
1206 std::vector<Point<dim>> points(n);
1207 std::vector<double> weights(n);
1208 const double length = p1.
distance(p2);
1210 for (
unsigned int k = 0; k < n; ++k)
1212 const double alpha = quadrature.
point(k)(0);
1213 points[k] = alpha * p2;
1214 points[k] += (1. - alpha) * p1;
1215 weights[k] = length * quadrature.
weight(k);
1225 const unsigned int face_no,
1226 const bool face_orientation,
1227 const bool face_flip,
1228 const bool face_rotation,
1229 const unsigned int n_quadrature_points)
1235 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1236 n_quadrature_points};
1239 const unsigned int orientation = (face_flip ? 4 : 0) +
1240 (face_rotation ? 2 : 0) +
1241 (face_orientation ? 1 : 0);
1242 return {(6 * face_no + orientation) * n_quadrature_points};
1255 return face_no * n_quadrature_points;
1280 static const unsigned int offset[2][2][2] = {
1299 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1300 n_quadrature_points);
1315 const unsigned int face_no,
1316 const bool face_orientation,
1317 const bool face_flip,
1318 const bool face_rotation,
1326 unsigned int offset = 0;
1329 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1330 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1331 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1332 static const std::array<unsigned int, 5> scale_pyramid = {
1343 if (quadrature.
size() == 1)
1344 offset =
scale[0] * quadrature[0].size() * face_no;
1346 for (
unsigned int i = 0; i < face_no; ++i)
1347 offset +=
scale[i] * quadrature[i].size();
1352 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1355 const unsigned int orientation = (face_flip ? 4 : 0) +
1356 (face_rotation ? 2 : 0) +
1357 (face_orientation ? 1 : 0);
1361 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1375 if (quadrature.
size() == 1)
1376 return quadrature[0].
size() * face_no;
1379 unsigned int result = 0;
1380 for (
unsigned int i = 0; i < face_no; ++i)
1381 result += quadrature[i].size();
1407 static const unsigned int offset[2][2][2] = {
1418 if (quadrature.
size() == 1)
1420 offset[face_orientation][face_flip][face_rotation] *
1422 quadrature[0].
size();
1425 unsigned int n_points_i = 0;
1426 for (
unsigned int i = 0; i < face_no; ++i)
1427 n_points_i += quadrature[i].size();
1429 unsigned int n_points = 0;
1430 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1432 n_points += quadrature[i].size();
1434 return (n_points_i +
1435 offset[face_orientation][face_flip][face_rotation] *
1452 const unsigned int face_no,
1453 const unsigned int subface_no,
1457 const unsigned int n_quadrature_points,
1468 n_quadrature_points);
1477 const unsigned int face_no,
1478 const unsigned int subface_no,
1482 const unsigned int n_quadrature_points,
1493 n_quadrature_points);
1502 const unsigned int face_no,
1503 const unsigned int subface_no,
1504 const bool face_orientation,
1505 const bool face_flip,
1506 const bool face_rotation,
1507 const unsigned int n_quadrature_points,
1510 const unsigned int dim = 3;
1543 const unsigned int total_subfaces_per_face = 8;
1560 static const unsigned int orientation_offset[2][2][2] = {
1589 static const unsigned int ref_case_offset[3] = {
1665 static const unsigned int
1701 equivalent_refine_case[ref_case][subface_no];
1702 const unsigned int equ_subface_no =
1703 equivalent_subface_number[ref_case][subface_no];
1710 (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1739 const unsigned int final_subface_no =
1748 return (((face_no * total_subfaces_per_face +
1749 ref_case_offset[final_ref_case - 1] + final_subface_no) +
1750 orientation_offset[face_orientation][face_flip][face_rotation]) *
1751 n_quadrature_points);
1760 const unsigned int face_no)
1766 std::vector<Point<dim>> points(quadrature.
size());
1777 const unsigned int face_no,
1778 const unsigned int subface_no,
1785 std::vector<Point<dim>> points(quadrature.
size());
1787 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_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
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 > &)
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< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
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)