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];
85 const unsigned int face_no,
97 const unsigned int face_no,
103 const unsigned int dim = 1;
107 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
115 const unsigned int face_no,
127 const unsigned int face_no,
130 const unsigned int dim = 2;
139 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
158 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
188 const unsigned int face_no,
200 const unsigned int face_no,
206 const unsigned int dim = 3;
211 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
249 const unsigned int face_no,
250 const unsigned int subface_no,
264 const unsigned int face_no,
272 const unsigned int dim = 1;
276 q_points[0] =
Point<dim>(
static_cast<double>(face_no));
284 const unsigned int face_no,
285 const unsigned int subface_no,
303 const unsigned int face_no,
304 const unsigned int subface_no,
308 const unsigned int dim = 2;
319 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
341 quadrature.
point(p)(0) / 2);
345 0.5 + quadrature.
point(p)(0) / 2);
371 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
446 const unsigned int face_no,
447 const unsigned int subface_no,
465 const unsigned int face_no,
466 const unsigned int subface_no,
473 const unsigned int dim = 3;
482 double const_value = face_no % 2;
491 const_index = face_no / 2;
497 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
498 eta_translation = 0.0;
522 xi_translation = subface_no % 2 * 0.5;
526 eta_translation = subface_no % 2 * 0.5;
531 xi_translation =
int(subface_no % 2) * 0.5;
532 eta_translation =
int(subface_no / 2) * 0.5;
540 for (
unsigned int p = 0; p < quadrature.
size(); ++p)
542 q_points[p][xi_index] =
543 xi_scale * quadrature.
point(p)(0) + xi_translation;
544 q_points[p][eta_index] =
545 eta_scale * quadrature.
point(p)(1) + eta_translation;
546 q_points[p][const_index] = const_value;
560 const unsigned int dim = 1;
565 std::vector<Point<dim>> q_points;
566 q_points.reserve(n_points * n_faces);
567 std::vector<Point<dim>> help(n_points);
572 for (
unsigned int face = 0; face < n_faces; ++face)
574 project_to_face(quadrature[quadrature.
size() == 1 ? 0 : face],
577 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
581 std::vector<double> weights;
582 weights.reserve(n_points * n_faces);
583 for (
unsigned int face = 0; face < n_faces; ++face)
585 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().begin(),
586 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().end(),
587 std::back_inserter(weights));
604 const auto support_points_line =
605 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<2>> {
607 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
611 return std::vector<Point<2>>(temp.begin(),
612 temp.begin() + face.first.size());
616 const std::array<std::pair<std::array<Point<2>, 2>,
double>, 3> faces = {
626 std::vector<Point<2>> points;
627 std::vector<double> weights;
630 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
632 for (
unsigned int orientation = 0; orientation < 2; ++orientation)
634 const auto &face = faces[face_no];
638 std::vector<Point<2>> support_points =
639 support_points_line(face, orientation);
642 const auto &sub_quadrature_points =
643 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
644 const auto &sub_quadrature_weights =
645 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
648 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
653 for (
unsigned int i = 0; i < 2; ++i)
656 poly.compute_value(i, sub_quadrature_points[j]);
658 points.emplace_back(mapped_point);
661 weights.emplace_back(sub_quadrature_weights[j] * face.second);
666 return {points, weights};
671 const unsigned int dim = 2;
675 unsigned int n_points_total = 0;
677 if (quadrature.
size() == 1)
682 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
683 n_points_total += quadrature[i].size();
687 std::vector<Point<dim>> q_points;
688 q_points.reserve(n_points_total);
689 std::vector<Point<dim>> help;
694 for (
unsigned int face = 0; face < n_faces; ++face)
696 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
697 project_to_face(quadrature[quadrature.
size() == 1 ? 0 : face],
700 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
704 std::vector<double> weights;
705 weights.reserve(n_points_total);
706 for (
unsigned int face = 0; face < n_faces; ++face)
708 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().begin(),
709 quadrature[quadrature.
size() == 1 ? 0 : face].get_weights().end(),
710 std::back_inserter(weights));
725 const auto support_points_tri =
726 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
728 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
732 return std::vector<Point<3>>(temp.begin(),
733 temp.begin() + face.first.size());
736 const auto support_points_quad =
737 [](
const auto &face,
const auto &orientation) -> std::vector<
Point<3>> {
739 std::copy_n(face.first.begin(), face.first.size(),
vertices.begin());
743 return std::vector<Point<3>>(temp.begin(),
744 temp.begin() + face.first.size());
747 const auto process = [&](
const auto &faces) {
749 std::vector<Point<3>> points;
750 std::vector<double> weights;
758 for (
unsigned int face_no = 0; face_no < faces.size(); ++face_no)
762 const unsigned int n_shape_functions = faces[face_no].first.size();
765 n_shape_functions == 3 ?
770 for (
unsigned int orientation = 0;
771 orientation < (n_shape_functions * 2);
774 const auto &face = faces[face_no];
776 const auto support_points =
777 n_shape_functions == 3 ? support_points_tri(face, orientation) :
778 support_points_quad(face, orientation);
781 const auto &sub_quadrature_points =
782 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_points();
783 const auto &sub_quadrature_weights =
784 quadrature[quadrature.
size() == 1 ? 0 : face_no].get_weights();
787 for (
unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
792 for (
unsigned int i = 0; i < n_shape_functions; ++i)
795 poly.compute_value(i, sub_quadrature_points[j]);
797 points.push_back(mapped_point);
800 const double scaling = [&]() {
801 const auto & supp_pts = support_points;
802 const unsigned int dim_ = 2;
803 const unsigned int spacedim = 3;
805 double result[spacedim][dim_];
807 std::vector<Tensor<1, dim_>> shape_derivatives(
810 for (
unsigned int i = 0; i < n_shape_functions; ++i)
811 shape_derivatives[i] =
812 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
814 for (
unsigned int i = 0; i < spacedim; ++i)
815 for (
unsigned int j = 0; j < dim_; ++j)
816 result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
817 for (
unsigned int k = 1; k < n_shape_functions; ++k)
818 for (
unsigned int i = 0; i < spacedim; ++i)
819 for (
unsigned int j = 0; j < dim_; ++j)
821 shape_derivatives[k][j] * supp_pts[k][i];
825 for (
unsigned int i = 0; i < spacedim; ++i)
826 for (
unsigned int j = 0; j < dim_; ++j)
827 contravariant[i][j] = result[i][j];
831 for (
unsigned int i = 0; i < spacedim; ++i)
832 for (
unsigned int j = 0; j < dim_; ++j)
833 DX_t[j][i] = contravariant[i][j];
836 for (
unsigned int i = 0; i < dim_; ++i)
837 for (
unsigned int j = 0; j < dim_; ++j)
838 G[i][j] = DX_t[i] * DX_t[j];
843 weights.push_back(sub_quadrature_weights[j] * scaling);
856 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
874 return process(faces);
878 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
903 return process(faces);
907 const std::vector<std::pair<std::vector<Point<3>>,
double>> faces = {
930 return process(faces);
936 const unsigned int dim = 3;
938 unsigned int n_points_total = 0;
940 if (quadrature.
size() == 1)
945 for (
unsigned int i = 0; i < quadrature.
size(); ++i)
946 n_points_total += quadrature[i].size();
952 std::vector<Point<dim>> q_points;
953 q_points.reserve(n_points_total);
954 std::vector<Point<dim>> help;
957 std::vector<double> weights;
958 weights.reserve(n_points_total);
964 for (
unsigned int i = 0; i < 8; ++i)
968 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
973 const auto quadrature_f =
974 quadrature[quadrature.
size() == 1 ? 0 : face];
978 mutation = quadrature_f;
990 mutation = internal::QProjector::reflect(quadrature_f);
994 internal::QProjector::reflect(quadrature_f), 3);
998 internal::QProjector::reflect(quadrature_f), 2);
1002 internal::QProjector::reflect(quadrature_f), 1);
1008 help.resize(quadrature[quadrature.
size() == 1 ? 0 : face].
size());
1009 project_to_face(mutation, face, help);
1010 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1014 std::back_inserter(weights));
1044 const unsigned int dim = 1;
1051 std::vector<Point<dim>> q_points;
1052 q_points.reserve(n_points * n_faces * subfaces_per_face);
1053 std::vector<Point<dim>> help(n_points);
1057 for (
unsigned int face = 0; face < n_faces; ++face)
1058 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1060 project_to_subface(quadrature, face, subface, help);
1061 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1065 std::vector<double> weights;
1066 weights.reserve(n_points * n_faces * subfaces_per_face);
1067 for (
unsigned int face = 0; face < n_faces; ++face)
1068 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1071 std::back_inserter(weights));
1073 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1075 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1094 const unsigned int dim = 2;
1096 const unsigned int n_points = quadrature.
size(),
1102 std::vector<Point<dim>> q_points;
1103 q_points.reserve(n_points * n_faces * subfaces_per_face);
1104 std::vector<Point<dim>> help(n_points);
1108 for (
unsigned int face = 0; face < n_faces; ++face)
1109 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1111 project_to_subface(quadrature, face, subface, help);
1112 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1116 std::vector<double> weights;
1117 weights.reserve(n_points * n_faces * subfaces_per_face);
1118 for (
unsigned int face = 0; face < n_faces; ++face)
1119 for (
unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1122 std::back_inserter(weights));
1124 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1126 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1154 const unsigned int dim = 3;
1155 SubQuadrature q_reflected = internal::QProjector::reflect(quadrature);
1165 const unsigned int n_points = quadrature.
size(),
1167 total_subfaces_per_face = 2 + 2 + 4;
1170 std::vector<Point<dim>> q_points;
1171 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1172 std::vector<Point<dim>> help(n_points);
1174 std::vector<double> weights;
1175 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1181 for (
const auto &mutation : q)
1185 for (
unsigned int face = 0; face < n_faces; ++face)
1189 for (
unsigned int subface = 0;
1194 project_to_subface(mutation,
1199 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1203 for (
unsigned int face = 0; face < n_faces; ++face)
1207 for (
unsigned int subface = 0;
1211 std::copy(mutation.get_weights().begin(),
1212 mutation.get_weights().end(),
1213 std::back_inserter(weights));
1216 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1218 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1239 const unsigned int child_no)
1241 return project_to_child(ReferenceCells::get_hypercube<dim>(),
1252 const unsigned int child_no)
1260 const unsigned int n_q_points = quadrature.
size();
1262 std::vector<Point<dim>> q_points(n_q_points);
1263 for (
unsigned int i = 0; i < n_q_points; ++i)
1271 std::vector<double> weights = quadrature.
get_weights();
1272 for (
unsigned int i = 0; i < n_q_points; ++i)
1284 return project_to_all_children(ReferenceCells::get_hypercube<dim>(),
1299 const unsigned int n_points = quadrature.
size(),
1302 std::vector<Point<dim>> q_points(n_points * n_children);
1303 std::vector<double> weights(n_points * n_children);
1307 for (
unsigned int child = 0; child < n_children; ++child)
1310 for (
unsigned int i = 0; i < n_points; ++i)
1312 q_points[child * n_points + i] = help.
point(i);
1313 weights[child * n_points + i] = help.
weight(i);
1327 return project_to_line(ReferenceCells::get_hypercube<dim>(),
1346 const unsigned int n = quadrature.
size();
1347 std::vector<Point<dim>> points(n);
1348 std::vector<double> weights(n);
1349 const double length = p1.
distance(p2);
1351 for (
unsigned int k = 0; k < n; ++k)
1353 const double alpha = quadrature.
point(k)(0);
1354 points[k] = alpha * p2;
1355 points[k] += (1. - alpha) * p1;
1356 weights[k] = length * quadrature.
weight(k);
1366 const bool face_orientation,
1367 const bool face_flip,
1368 const bool face_rotation,
1369 const unsigned int n_quadrature_points)
1371 return face(ReferenceCells::get_hypercube<dim>(),
1376 n_quadrature_points);
1384 const unsigned int face_no,
1385 const bool face_orientation,
1386 const bool face_flip,
1387 const bool face_rotation,
1388 const unsigned int n_quadrature_points)
1394 return {(2 * face_no + face_orientation) * n_quadrature_points};
1397 const unsigned int orientation =
1398 (face_flip * 2 + face_rotation) * 2 + face_orientation;
1399 return {(6 * face_no + orientation) * n_quadrature_points};
1412 return face_no * n_quadrature_points;
1437 static const unsigned int offset[2][2][2] = {
1456 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1457 n_quadrature_points);
1472 const unsigned int face_no,
1473 const bool face_orientation,
1474 const bool face_flip,
1475 const bool face_rotation,
1483 unsigned int offset = 0;
1486 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1487 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1488 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1489 static const std::array<unsigned int, 5> scale_pyramid = {
1500 if (quadrature.
size() == 1)
1501 offset =
scale[0] * quadrature[0].size() * face_no;
1503 for (
unsigned int i = 0; i < face_no; ++i)
1504 offset +=
scale[i] * quadrature[i].size();
1509 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1512 const unsigned int orientation =
1513 (face_flip * 2 + face_rotation) * 2 + face_orientation;
1517 quadrature[quadrature.
size() == 1 ? 0 : face_no].
size()};
1531 if (quadrature.
size() == 1)
1532 return quadrature[0].
size() * face_no;
1535 unsigned int result = 0;
1536 for (
unsigned int i = 0; i < face_no; ++i)
1537 result += quadrature[i].size();
1563 static const unsigned int offset[2][2][2] = {
1574 if (quadrature.
size() == 1)
1576 offset[face_orientation][face_flip][face_rotation] *
1578 quadrature[0].
size();
1581 unsigned int n_points_i = 0;
1582 for (
unsigned int i = 0; i < face_no; ++i)
1583 n_points_i += quadrature[i].size();
1585 unsigned int n_points = 0;
1586 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1588 n_points += quadrature[i].size();
1590 return (n_points_i +
1591 offset[face_orientation][face_flip][face_rotation] *
1608 const unsigned int face_no,
1609 const unsigned int subface_no,
1613 const unsigned int n_quadrature_points,
1624 n_quadrature_points);
1632 const unsigned int face_no,
1633 const unsigned int subface_no,
1634 const bool face_orientation,
1635 const bool face_flip,
1636 const bool face_rotation,
1637 const unsigned int n_quadrature_points,
1646 n_quadrature_points,
1656 const unsigned int face_no,
1657 const unsigned int subface_no,
1661 const unsigned int n_quadrature_points,
1672 n_quadrature_points);
1680 const unsigned int face_no,
1681 const unsigned int subface_no,
1682 const bool face_orientation,
1683 const bool face_flip,
1684 const bool face_rotation,
1685 const unsigned int n_quadrature_points,
1694 n_quadrature_points,
1703 const unsigned int face_no,
1704 const unsigned int subface_no,
1705 const bool face_orientation,
1706 const bool face_flip,
1707 const bool face_rotation,
1708 const unsigned int n_quadrature_points,
1711 const unsigned int dim = 3;
1744 const unsigned int total_subfaces_per_face = 8;
1761 static const unsigned int orientation_offset[2][2][2] = {
1790 static const unsigned int ref_case_offset[3] = {
1866 static const unsigned int
1902 equivalent_refine_case[ref_case][subface_no];
1903 const unsigned int equ_subface_no =
1904 equivalent_subface_number[ref_case][subface_no];
1911 (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1940 const unsigned int final_subface_no =
1949 return (((face_no * total_subfaces_per_face +
1950 ref_case_offset[final_ref_case - 1] + final_subface_no) +
1951 orientation_offset[face_orientation][face_flip][face_rotation]) *
1952 n_quadrature_points);
1959 const unsigned int face_no,
1960 const unsigned int subface_no,
1961 const bool face_orientation,
1962 const bool face_flip,
1963 const bool face_rotation,
1964 const unsigned int n_quadrature_points,
1973 n_quadrature_points,
1982 const unsigned int face_no)
1995 const unsigned int face_no)
2001 std::vector<Point<dim>> points(quadrature.
size());
2011 const unsigned int face_no,
2012 const unsigned int subface_no,
2028 const unsigned int face_no,
2029 const unsigned int subface_no,
2036 std::vector<Point<dim>> points(quadrature.
size());
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 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 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 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_faces(const Quadrature< dim - 1 > &quadrature)
static Quadrature< dim > project_to_line(const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
static void project_to_face(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::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)
unsigned int size() 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)
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
::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)