From d36d2f01d9555c7b9df89f1d3e463c5cc4e8b2c3 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 18 Jul 2020 22:20:46 +0200 Subject: [PATCH] Intoduce ReferenceCell in QProjector --- include/deal.II/base/qprojector.h | 296 ++++- include/deal.II/fe/fe_poly.templates.h | 10 +- include/deal.II/fe/fe_poly_face.h | 2 +- include/deal.II/fe/fe_tools.templates.h | 10 +- include/deal.II/grid/reference_cell.h | 46 + source/base/CMakeLists.txt | 1 + source/base/qprojector.cc | 1509 +++++++++++++++++++++++ source/base/quadrature.cc | 1102 ----------------- source/fe/fe.cc | 6 +- source/fe/fe_abf.cc | 24 +- source/fe/fe_bdm.cc | 14 +- source/fe/fe_bernstein.cc | 9 +- source/fe/fe_nedelec.cc | 45 +- source/fe/fe_nedelec_sz.cc | 11 +- source/fe/fe_p1nc.cc | 8 +- source/fe/fe_poly_tensor.cc | 10 +- source/fe/fe_q_base.cc | 31 +- source/fe/fe_q_bubbles.cc | 1 - source/fe/fe_q_dg0.cc | 1 - source/fe/fe_raviart_thomas.cc | 24 +- source/fe/fe_raviart_thomas_nodal.cc | 23 +- source/fe/fe_rt_bubbles.cc | 14 +- source/fe/fe_system.cc | 1 - source/fe/fe_trace.cc | 1 - source/fe/mapping_cartesian.cc | 42 +- source/fe/mapping_fe_field.cc | 29 +- source/fe/mapping_manifold.cc | 13 +- source/fe/mapping_q1.cc | 1 - source/fe/mapping_q_generic.cc | 13 +- source/grid/grid_out.cc | 4 +- 30 files changed, 2073 insertions(+), 1228 deletions(-) create mode 100644 source/base/qprojector.cc diff --git a/include/deal.II/base/qprojector.h b/include/deal.II/base/qprojector.h index af3818566e..7e231242f7 100644 --- a/include/deal.II/base/qprojector.h +++ b/include/deal.II/base/qprojector.h @@ -22,6 +22,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN @@ -86,20 +88,49 @@ public: * Compute the quadrature points on the cell if the given quadrature formula * is used on face face_no. For further details, see the general * doc for this class. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static void + DEAL_II_DEPRECATED static void project_to_face(const SubQuadrature & quadrature, const unsigned int face_no, std::vector> &q_points); + /** + * Compute the quadrature points on the cell if the given quadrature formula + * is used on face face_no. For further details, see the general + * doc for this class. + */ + static void + project_to_face(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no, + std::vector> & q_points); + /** * Compute the cell quadrature formula corresponding to using * quadrature on face face_no. For further details, see * the general doc for this class. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_face(const SubQuadrature &quadrature, const unsigned int face_no); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on face face_no. For further details, see + * the general doc for this class. + */ + static Quadrature + project_to_face(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no); + /** * Compute the quadrature points on the cell if the given quadrature formula * is used on face face_no, subface number subface_no @@ -108,8 +139,12 @@ public: * * @note Only the points are transformed. The quadrature weights are the * same as those of the original rule. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static void + DEAL_II_DEPRECATED static void project_to_subface(const SubQuadrature & quadrature, const unsigned int face_no, const unsigned int subface_no, @@ -117,6 +152,24 @@ public: const RefinementCase &ref_case = RefinementCase::isotropic_refinement); + /** + * Compute the quadrature points on the cell if the given quadrature formula + * is used on face face_no, subface number subface_no + * corresponding to RefineCase::Type ref_case. The last argument is + * only used in 3D. + * + * @note Only the points are transformed. The quadrature weights are the + * same as those of the original rule. + */ + static void + project_to_subface(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase &ref_case = + RefinementCase::isotropic_refinement); + /** * Compute the cell quadrature formula corresponding to using * quadrature on subface subface_no of face @@ -125,14 +178,39 @@ public: * * @note Only the points are transformed. The quadrature weights are the * same as those of the original rule. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_subface(const SubQuadrature & quadrature, const unsigned int face_no, const unsigned int subface_no, const RefinementCase &ref_case = RefinementCase::isotropic_refinement); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on subface subface_no of face + * face_no with RefinementCase ref_case. The last + * argument is only used in 3D. + * + * @note Only the points are transformed. The quadrature weights are the + * same as those of the original rule. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. + */ + static Quadrature + project_to_subface(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const RefinementCase &ref_case = + RefinementCase::isotropic_refinement); + /** * Take a face quadrature formula and generate a cell quadrature formula * from it where the quadrature points of the given argument are projected @@ -148,10 +226,34 @@ public: * * @note In 3D, this function produces eight sets of quadrature points for * each face, in order to cope possibly different orientations of the mesh. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_all_faces(const SubQuadrature &quadrature); + /** + * Take a face quadrature formula and generate a cell quadrature formula + * from it where the quadrature points of the given argument are projected + * on all faces. + * + * The weights of the new rule are replications of the original weights. + * Thus, the sum of the weights is not one, but the number of faces, which + * is the surface of the reference cell. + * + * This in particular allows us to extract a subset of points corresponding + * to a single face and use it as a quadrature on this face, as is done in + * FEFaceValues. + * + * @note In 3D, this function produces eight sets of quadrature points for + * each face, in order to cope possibly different orientations of the mesh. + */ + static Quadrature + project_to_all_faces(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature); + /** * Take a face quadrature formula and generate a cell quadrature formula * from it where the quadrature points of the given argument are projected @@ -164,10 +266,31 @@ public: * This in particular allows us to extract a subset of points corresponding * to a single subface and use it as a quadrature on this face, as is done * in FESubfaceValues. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_all_subfaces(const SubQuadrature &quadrature); + /** + * Take a face quadrature formula and generate a cell quadrature formula + * from it where the quadrature points of the given argument are projected + * on all subfaces. + * + * Like in project_to_all_faces(), the weights of the new rule sum up to the + * number of faces (not subfaces), which is the surface of the reference + * cell. + * + * This in particular allows us to extract a subset of points corresponding + * to a single subface and use it as a quadrature on this face, as is done + * in FESubfaceValues. + */ + static Quadrature + project_to_all_subfaces(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature); + /** * Project a given quadrature formula to a child of a cell. You may want to * use this function in case you want to extend an integral only over the @@ -177,11 +300,30 @@ public: * As integration using this quadrature formula now only extends over a * fraction of the cell, the weights of the resulting object are divided by * GeometryInfo::children_per_cell. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_child(const Quadrature &quadrature, const unsigned int child_no); + /** + * Project a given quadrature formula to a child of a cell. You may want to + * use this function in case you want to extend an integral only over the + * area which a potential child would occupy. The child numbering is the + * same as the children would be numbered upon refinement of the cell. + * + * As integration using this quadrature formula now only extends over a + * fraction of the cell, the weights of the resulting object are divided by + * GeometryInfo::children_per_cell. + */ + static Quadrature + project_to_child(const ReferenceCell::Type reference_cell_type, + const Quadrature & quadrature, + const unsigned int child_no); + /** * Project a quadrature rule to all children of a cell. Similarly to * project_to_all_subfaces(), this function replicates the formula generated @@ -190,19 +332,50 @@ public: * * The child numbering is the same as the children would be numbered upon * refinement of the cell. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_all_children(const Quadrature &quadrature); + /** + * Project a quadrature rule to all children of a cell. Similarly to + * project_to_all_subfaces(), this function replicates the formula generated + * by project_to_child() for all children, such that the weights sum up to + * one, the volume of the total cell again. + * + * The child numbering is the same as the children would be numbered upon + * refinement of the cell. + */ + static Quadrature + project_to_all_children(const ReferenceCell::Type reference_cell_type, + const Quadrature & quadrature); + /** * Project the one dimensional rule quadrature to the straight line * connecting the points p1 and p2. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static Quadrature + DEAL_II_DEPRECATED static Quadrature project_to_line(const Quadrature<1> &quadrature, const Point & p1, const Point & p2); + /** + * Project the one dimensional rule quadrature to the straight line + * connecting the points p1 and p2. + */ + static Quadrature + project_to_line(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const Point & p1, + const Point & p2); + /** * Since the project_to_all_faces() and project_to_all_subfaces() functions * chain together the quadrature points and weights of all projections of a @@ -241,14 +414,36 @@ public: * The last argument denotes the number of quadrature points the lower- * dimensional face quadrature formula (the one that has been projected * onto the faces) has. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static DataSetDescriptor + DEAL_II_DEPRECATED 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 function to generate an offset object for a given face of a cell + * with the given face orientation, flip and rotation. This function of + * course is only allowed if dim>=2, and the face orientation, + * flip and rotation are ignored if the space dimension equals 2. + * + * The last argument denotes the number of quadrature points the lower- + * dimensional face quadrature formula (the one that has been projected + * onto the faces) has. + */ + static DataSetDescriptor + face(const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const unsigned int n_quadrature_points); + /** * Static function to generate an offset object for a given subface of a * cell with the given face orientation, flip and rotation. This function @@ -260,8 +455,12 @@ public: * projected onto the faces) has. * * Through the last argument anisotropic refinement can be respected. + * + * @deprecated This function makes an implicit assumption that the cell is + * a line (1D), a quad (2D), or a hex (3D). Use the other version of + * this function that takes the reference cell type instead. */ - static DataSetDescriptor + DEAL_II_DEPRECATED static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, @@ -271,6 +470,29 @@ public: const internal::SubfaceCase ref_case = internal::SubfaceCase::case_isotropic); + /** + * Static function to generate an offset object for a given subface of a + * cell with the given face orientation, flip and rotation. This function + * of course is only allowed if dim>=2, and the face orientation, + * flip and rotation are ignored if the space dimension equals 2. + * + * The last but one argument denotes the number of quadrature points the + * lower-dimensional face quadrature formula (the one that has been + * projected onto the faces) has. + * + * Through the last argument anisotropic refinement can be respected. + */ + static DataSetDescriptor + subface(const ReferenceCell::Type reference_cell_type, + 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 ref_case = + internal::SubfaceCase::case_isotropic); + /** * Conversion operator to an integer denoting the offset of the first * element of this dataset in the set of quadrature formulas all projected @@ -365,18 +587,41 @@ QProjector<1>::project_to_face(const Quadrature<0> &, std::vector> &); template <> void +QProjector<1>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<0> &, + const unsigned int, + std::vector> &); +template <> +void QProjector<2>::project_to_face(const Quadrature<1> & quadrature, const unsigned int face_no, std::vector> &q_points); template <> void +QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const unsigned int face_no, + std::vector> & q_points); +template <> +void QProjector<3>::project_to_face(const Quadrature<2> & quadrature, const unsigned int face_no, std::vector> &q_points); +template <> +void +QProjector<3>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<2> & quadrature, + const unsigned int face_no, + std::vector> & q_points); template <> Quadrature<1> QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature); +template <> +Quadrature<1> +QProjector<1>::project_to_all_faces( + const ReferenceCell::Type reference_cell_type, + const Quadrature<0> & quadrature); template <> @@ -388,6 +633,14 @@ QProjector<1>::project_to_subface(const Quadrature<0> &, const RefinementCase<0> &); template <> void +QProjector<1>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<0> &, + const unsigned int, + const unsigned int, + std::vector> &, + const RefinementCase<0> &); +template <> +void QProjector<2>::project_to_subface(const Quadrature<1> & quadrature, const unsigned int face_no, const unsigned int subface_no, @@ -395,15 +648,36 @@ QProjector<2>::project_to_subface(const Quadrature<1> & quadrature, const RefinementCase<1> &); template <> void +QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<1> &); +template <> +void QProjector<3>::project_to_subface(const Quadrature<2> & quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector> & q_points, const RefinementCase<2> &face_ref_case); +template <> +void +QProjector<3>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<2> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<2> & face_ref_case); template <> Quadrature<1> QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature); +template <> +Quadrature<1> +QProjector<1>::project_to_all_subfaces( + const ReferenceCell::Type reference_cell_type, + const Quadrature<0> & quadrature); #endif // DOXYGEN diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index c4bcb2c100..cf88f802a7 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -358,8 +358,9 @@ FE_Poly::fill_fe_face_values( // to take (all data sets for all // faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::face(face_no, + const auto offset = + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), @@ -450,8 +451,9 @@ FE_Poly::fill_fe_subface_values( // to take (all data sets for all // sub-faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::subface(face_no, + const auto offset = + QProjector::DataSetDescriptor::subface(this->reference_cell_type(), + face_no, sub_no, cell->face_orientation(face_no), cell->face_flip(face_no), diff --git a/include/deal.II/fe/fe_poly_face.h b/include/deal.II/fe/fe_poly_face.h index d2b3834c33..cd2be76ff1 100644 --- a/include/deal.II/fe/fe_poly_face.h +++ b/include/deal.II/fe/fe_poly_face.h @@ -167,7 +167,7 @@ protected: return get_face_data(update_flags, mapping, QProjector::project_to_all_children( - quadrature), + ReferenceCell::get_hypercube(dim - 1), quadrature), output_data); } diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index 02f472e41d..34c7d92789 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -2035,7 +2035,9 @@ namespace FETools // hating the anisotropic implementation QGauss q_gauss(degree + 1); const Quadrature q_fine = - QProjector::project_to_face(q_gauss, face_fine); + QProjector::project_to_face(fe.reference_cell_type(), + q_gauss, + face_fine); const unsigned int nq = q_fine.size(); FEValues fine(mapping, @@ -2084,10 +2086,8 @@ namespace FETools cell_number < GeometryInfo::max_children_per_face; ++cell_number) { - const Quadrature q_coarse = - QProjector::project_to_subface(q_gauss, - face_coarse, - cell_number); + const Quadrature q_coarse = QProjector::project_to_subface( + fe.reference_cell_type(), q_gauss, face_coarse, cell_number); FEValues coarse(mapping, fe, q_coarse, update_values); typename Triangulation::active_cell_iterator fine_cell = diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 8832704e63..d678ce0f61 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -45,6 +45,52 @@ namespace ReferenceCell Invalid = static_cast(-1) }; + /** + * Return the correct simplex reference cell type for the given dimension + * @p dim. + */ + inline Type + get_simplex(const unsigned int dim) + { + switch (dim) + { + case 0: + return Type::Vertex; + case 1: + return Type::Line; + case 2: + return Type::Tri; + case 3: + return Type::Tet; + default: + Assert(false, ExcNotImplemented()); + return Type::Invalid; + } + } + + /** + * Return the correct hypercube reference cell type for the given dimension + * @p dim. + */ + inline Type + get_hypercube(const unsigned int dim) + { + switch (dim) + { + case 0: + return Type::Vertex; + case 1: + return Type::Line; + case 2: + return Type::Quad; + case 3: + return Type::Hex; + default: + Assert(false, ExcNotImplemented()); + return Type::Invalid; + } + } + namespace internal { /** diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index a046beecc8..5a56310fe7 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -76,6 +76,7 @@ SET(_unity_include_src polynomials_raviart_thomas.cc polynomials_rt_bubbles.cc process_grid.cc + qprojector.cc quadrature.cc quadrature_lib.cc quadrature_selector.cc diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc new file mode 100644 index 0000000000..275a9356be --- /dev/null +++ b/source/base/qprojector.cc @@ -0,0 +1,1509 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +template +Quadrature<2> +QProjector::reflect(const Quadrature<2> &q) +{ + std::vector> q_points(q.size()); + std::vector weights(q.size()); + for (unsigned int i = 0; i < q.size(); ++i) + { + q_points[i][0] = q.point(i)[1]; + q_points[i][1] = q.point(i)[0]; + + weights[i] = q.weight(i); + } + + return Quadrature<2>(q_points, weights); +} + + +template +Quadrature<2> +QProjector::rotate(const Quadrature<2> &q, const unsigned int n_times) +{ + std::vector> q_points(q.size()); + std::vector weights(q.size()); + for (unsigned int i = 0; i < q.size(); ++i) + { + switch (n_times % 4) + { + case 0: + // 0 degree + q_points[i][0] = q.point(i)[0]; + q_points[i][1] = q.point(i)[1]; + break; + case 1: + // 90 degree counterclockwise + q_points[i][0] = 1.0 - q.point(i)[1]; + q_points[i][1] = q.point(i)[0]; + break; + case 2: + // 180 degree counterclockwise + q_points[i][0] = 1.0 - q.point(i)[0]; + q_points[i][1] = 1.0 - q.point(i)[1]; + break; + case 3: + // 270 degree counterclockwise + q_points[i][0] = q.point(i)[1]; + q_points[i][1] = 1.0 - q.point(i)[0]; + break; + } + + weights[i] = q.weight(i); + } + + return Quadrature<2>(q_points, weights); +} + + +template <> +void +QProjector<1>::project_to_face(const Quadrature<0> & quadrature, + const unsigned int face_no, + std::vector> &q_points) +{ + project_to_face(ReferenceCell::Type::Line, quadrature, face_no, q_points); +} + + + +template <> +void +QProjector<1>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<0> &, + const unsigned int face_no, + std::vector> &q_points) +{ + Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 1; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertDimension(q_points.size(), 1); + + q_points[0] = Point(static_cast(face_no)); +} + + + +template <> +void +QProjector<2>::project_to_face(const Quadrature<1> & quadrature, + const unsigned int face_no, + std::vector> &q_points) +{ + project_to_face(ReferenceCell::Type::Quad, quadrature, face_no, q_points); +} + + + +template <> +void +QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const unsigned int face_no, + std::vector> & q_points) +{ + Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 2; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + Assert(q_points.size() == quadrature.size(), + ExcDimensionMismatch(q_points.size(), quadrature.size())); + + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + q_points[p] = Point(0, quadrature.point(p)(0)); + break; + case 1: + q_points[p] = Point(1, quadrature.point(p)(0)); + break; + case 2: + q_points[p] = Point(quadrature.point(p)(0), 0); + break; + case 3: + q_points[p] = Point(quadrature.point(p)(0), 1); + break; + default: + Assert(false, ExcInternalError()); + } +} + + + +template <> +void +QProjector<3>::project_to_face(const Quadrature<2> & quadrature, + const unsigned int face_no, + std::vector> &q_points) +{ + project_to_face(ReferenceCell::Type::Hex, quadrature, face_no, q_points); +} + + + +template <> +void +QProjector<3>::project_to_face(const ReferenceCell::Type reference_cell_type, + const Quadrature<2> & quadrature, + const unsigned int face_no, + std::vector> & q_points) +{ + Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 3; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + Assert(q_points.size() == quadrature.size(), + ExcDimensionMismatch(q_points.size(), quadrature.size())); + + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + q_points[p] = + Point(0, quadrature.point(p)(0), quadrature.point(p)(1)); + break; + case 1: + q_points[p] = + Point(1, quadrature.point(p)(0), quadrature.point(p)(1)); + break; + case 2: + q_points[p] = + Point(quadrature.point(p)(1), 0, quadrature.point(p)(0)); + break; + case 3: + q_points[p] = + Point(quadrature.point(p)(1), 1, quadrature.point(p)(0)); + break; + case 4: + q_points[p] = + Point(quadrature.point(p)(0), quadrature.point(p)(1), 0); + break; + case 5: + q_points[p] = + Point(quadrature.point(p)(0), quadrature.point(p)(1), 1); + break; + + default: + Assert(false, ExcInternalError()); + } +} + + + +template <> +void +QProjector<1>::project_to_subface(const Quadrature<0> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<0> &ref_case) +{ + project_to_subface(ReferenceCell::Type::Line, + quadrature, + face_no, + subface_no, + q_points, + ref_case); +} + + + +template <> +void +QProjector<1>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<0> &, + const unsigned int face_no, + const unsigned int, + std::vector> &q_points, + const RefinementCase<0> &) +{ + Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 1; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertDimension(q_points.size(), 1); + + q_points[0] = Point(static_cast(face_no)); +} + + + +template <> +void +QProjector<2>::project_to_subface(const Quadrature<1> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<1> &ref_case) +{ + project_to_subface(ReferenceCell::Type::Quad, + quadrature, + face_no, + subface_no, + q_points, + ref_case); +} + + + +template <> +void +QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<1> &) +{ + Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 2; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertIndexRange(subface_no, GeometryInfo::max_children_per_face); + + Assert(q_points.size() == quadrature.size(), + ExcDimensionMismatch(q_points.size(), quadrature.size())); + + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + switch (subface_no) + { + case 0: + q_points[p] = Point(0, quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = Point(0, quadrature.point(p)(0) / 2 + 0.5); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 1: + switch (subface_no) + { + case 0: + q_points[p] = Point(1, quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = Point(1, quadrature.point(p)(0) / 2 + 0.5); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 2: + switch (subface_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0) / 2, 0); + break; + case 1: + q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 0); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 3: + switch (subface_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0) / 2, 1); + break; + case 1: + q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 1); + break; + default: + Assert(false, ExcInternalError()); + } + break; + + default: + Assert(false, ExcInternalError()); + } +} + + + +template <> +void +QProjector<3>::project_to_subface(const Quadrature<2> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<2> &ref_case) +{ + project_to_subface(ReferenceCell::Type::Hex, + quadrature, + face_no, + subface_no, + q_points, + ref_case); +} + + + +template <> +void +QProjector<3>::project_to_subface(const ReferenceCell::Type reference_cell_type, + const Quadrature<2> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + std::vector> & q_points, + const RefinementCase<2> & ref_case) +{ + Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 3; + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertIndexRange(subface_no, GeometryInfo::max_children_per_face); + Assert(q_points.size() == quadrature.size(), + ExcDimensionMismatch(q_points.size(), quadrature.size())); + + // one coordinate is at a const value. for + // faces 0, 2 and 4 this value is 0.0, for + // faces 1, 3 and 5 it is 1.0 + double const_value = face_no % 2; + // local 2d coordinates are xi and eta, + // global 3d coordinates are x, y and + // z. those have to be mapped. the following + // indices tell, which global coordinate + // (0->x, 1->y, 2->z) corresponds to which + // local one + unsigned int xi_index = numbers::invalid_unsigned_int, + eta_index = numbers::invalid_unsigned_int, + const_index = face_no / 2; + // the xi and eta values have to be scaled + // (by factor 0.5 or factor 1.0) depending on + // the refinement case and translated (by 0.0 + // or 0.5) depending on the refinement case + // and subface_no. + double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0, + eta_translation = 0.0; + // set the index mapping between local and + // global coordinates + switch (face_no / 2) + { + case 0: + xi_index = 1; + eta_index = 2; + break; + case 1: + xi_index = 2; + eta_index = 0; + break; + case 2: + xi_index = 0; + eta_index = 1; + break; + } + // set the scale and translation parameter + // for individual subfaces + switch (ref_case) + { + case RefinementCase::cut_x: + xi_scale = 0.5; + xi_translation = subface_no % 2 * 0.5; + break; + case RefinementCase::cut_y: + eta_scale = 0.5; + eta_translation = subface_no % 2 * 0.5; + break; + case RefinementCase::cut_xy: + xi_scale = 0.5; + eta_scale = 0.5; + xi_translation = int(subface_no % 2) * 0.5; + eta_translation = int(subface_no / 2) * 0.5; + break; + default: + Assert(false, ExcInternalError()); + break; + } + // finally, compute the scaled, translated, + // projected quadrature points + for (unsigned int p = 0; p < quadrature.size(); ++p) + { + q_points[p][xi_index] = + xi_scale * quadrature.point(p)(0) + xi_translation; + q_points[p][eta_index] = + eta_scale * quadrature.point(p)(1) + eta_translation; + q_points[p][const_index] = const_value; + } +} + + +template <> +Quadrature<1> +QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature) +{ + return project_to_all_faces(ReferenceCell::Type::Line, quadrature); +} + + +template <> +Quadrature<1> +QProjector<1>::project_to_all_faces( + const ReferenceCell::Type reference_cell_type, + const Quadrature<0> & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 1; + + const unsigned int n_points = 1, n_faces = GeometryInfo::faces_per_cell; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces); + std::vector> help(n_points); + + + // project to each face and append + // results + for (unsigned int face = 0; face < n_faces; ++face) + { + project_to_face(quadrature, face, help); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + std::vector weights; + weights.reserve(n_points * n_faces); + for (unsigned int face = 0; face < n_faces; ++face) + std::copy(quadrature.get_weights().begin(), + quadrature.get_weights().end(), + std::back_inserter(weights)); + + Assert(q_points.size() == n_points * n_faces, ExcInternalError()); + Assert(weights.size() == n_points * n_faces, ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<2> +QProjector<2>::project_to_all_faces( + const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 2; + + const unsigned int n_points = quadrature.size(), + n_faces = GeometryInfo::faces_per_cell; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces); + std::vector> help(n_points); + + // project to each face and append + // results + for (unsigned int face = 0; face < n_faces; ++face) + { + project_to_face(quadrature, face, help); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + std::vector weights; + weights.reserve(n_points * n_faces); + for (unsigned int face = 0; face < n_faces; ++face) + std::copy(quadrature.get_weights().begin(), + quadrature.get_weights().end(), + std::back_inserter(weights)); + + Assert(q_points.size() == n_points * n_faces, ExcInternalError()); + Assert(weights.size() == n_points * n_faces, ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<2> +QProjector<2>::project_to_all_faces(const SubQuadrature &quadrature) +{ + return project_to_all_faces(ReferenceCell::Type::Quad, quadrature); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_all_faces( + const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 3; + + SubQuadrature q_reflected = reflect(quadrature); + SubQuadrature q[8] = {quadrature, + rotate(quadrature, 1), + rotate(quadrature, 2), + rotate(quadrature, 3), + q_reflected, + rotate(q_reflected, 3), + rotate(q_reflected, 2), + rotate(q_reflected, 1)}; + + + + const unsigned int n_points = quadrature.size(), + n_faces = GeometryInfo::faces_per_cell; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces * 8); + std::vector> help(n_points); + + std::vector weights; + weights.reserve(n_points * n_faces * 8); + + // do the following for all possible + // mutations of a face (mutation==0 + // corresponds to a face with standard + // orientation, no flip and no rotation) + for (const auto &mutation : q) + { + // project to each face and append + // results + for (unsigned int face = 0; face < n_faces; ++face) + { + project_to_face(mutation, face, help); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + for (unsigned int face = 0; face < n_faces; ++face) + std::copy(mutation.get_weights().begin(), + mutation.get_weights().end(), + std::back_inserter(weights)); + } + + + Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError()); + Assert(weights.size() == n_points * n_faces * 8, ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_all_faces(const SubQuadrature &quadrature) +{ + return project_to_all_faces(ReferenceCell::Type::Hex, quadrature); +} + + + +template <> +Quadrature<1> +QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature) +{ + return project_to_all_subfaces(ReferenceCell::Type::Line, quadrature); +} + + + +template <> +Quadrature<1> +QProjector<1>::project_to_all_subfaces( + const ReferenceCell::Type reference_cell_type, + const Quadrature<0> & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 1; + + const unsigned int n_points = 1, n_faces = GeometryInfo::faces_per_cell, + subfaces_per_face = + GeometryInfo::max_children_per_face; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces * subfaces_per_face); + std::vector> help(n_points); + + // project to each face and copy + // results + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) + { + project_to_subface(quadrature, face, subface, help); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + std::vector weights; + weights.reserve(n_points * n_faces * subfaces_per_face); + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) + std::copy(quadrature.get_weights().begin(), + quadrature.get_weights().end(), + std::back_inserter(weights)); + + Assert(q_points.size() == n_points * n_faces * subfaces_per_face, + ExcInternalError()); + Assert(weights.size() == n_points * n_faces * subfaces_per_face, + ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<2> +QProjector<2>::project_to_all_subfaces( + const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 2; + + const unsigned int n_points = quadrature.size(), + n_faces = GeometryInfo::faces_per_cell, + subfaces_per_face = + GeometryInfo::max_children_per_face; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces * subfaces_per_face); + std::vector> help(n_points); + + // project to each face and copy + // results + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) + { + project_to_subface(quadrature, face, subface, help); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + std::vector weights; + weights.reserve(n_points * n_faces * subfaces_per_face); + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) + std::copy(quadrature.get_weights().begin(), + quadrature.get_weights().end(), + std::back_inserter(weights)); + + Assert(q_points.size() == n_points * n_faces * subfaces_per_face, + ExcInternalError()); + Assert(weights.size() == n_points * n_faces * subfaces_per_face, + ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<2> +QProjector<2>::project_to_all_subfaces(const SubQuadrature &quadrature) +{ + return project_to_all_subfaces(ReferenceCell::Type::Quad, quadrature); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_all_subfaces( + const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int dim = 3; + SubQuadrature q_reflected = reflect(quadrature); + SubQuadrature q[8] = {quadrature, + rotate(quadrature, 1), + rotate(quadrature, 2), + rotate(quadrature, 3), + q_reflected, + rotate(q_reflected, 3), + rotate(q_reflected, 2), + rotate(q_reflected, 1)}; + + const unsigned int n_points = quadrature.size(), + n_faces = GeometryInfo::faces_per_cell, + total_subfaces_per_face = 2 + 2 + 4; + + // first fix quadrature points + std::vector> q_points; + q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8); + std::vector> help(n_points); + + std::vector weights; + weights.reserve(n_points * n_faces * total_subfaces_per_face * 8); + + // do the following for all possible + // mutations of a face (mutation==0 + // corresponds to a face with standard + // orientation, no flip and no rotation) + for (const auto &mutation : q) + { + // project to each face and copy + // results + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int ref_case = RefinementCase::cut_xy; + ref_case >= RefinementCase::cut_x; + --ref_case) + for (unsigned int subface = 0; + subface < GeometryInfo::n_children( + RefinementCase(ref_case)); + ++subface) + { + project_to_subface(mutation, + face, + subface, + help, + RefinementCase(ref_case)); + std::copy(help.begin(), help.end(), std::back_inserter(q_points)); + } + + // next copy over weights + for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int ref_case = RefinementCase::cut_xy; + ref_case >= RefinementCase::cut_x; + --ref_case) + for (unsigned int subface = 0; + subface < GeometryInfo::n_children( + RefinementCase(ref_case)); + ++subface) + std::copy(mutation.get_weights().begin(), + mutation.get_weights().end(), + std::back_inserter(weights)); + } + + Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8, + ExcInternalError()); + Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8, + ExcInternalError()); + + return Quadrature(q_points, weights); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_all_subfaces(const SubQuadrature &quadrature) +{ + return project_to_all_subfaces(ReferenceCell::Type::Hex, quadrature); +} + + + +// This function is not used in the library +template +Quadrature +QProjector::project_to_child(const Quadrature &quadrature, + const unsigned int child_no) +{ + return project_to_child(ReferenceCell::get_hypercube(dim), + quadrature, + child_no); +} + + + +template +Quadrature +QProjector::project_to_child(const ReferenceCell::Type reference_cell_type, + const Quadrature & quadrature, + const unsigned int child_no) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + AssertIndexRange(child_no, GeometryInfo::max_children_per_cell); + + const unsigned int n_q_points = quadrature.size(); + + std::vector> q_points(n_q_points); + for (unsigned int i = 0; i < n_q_points; ++i) + q_points[i] = + GeometryInfo::child_to_cell_coordinates(quadrature.point(i), + child_no); + + // for the weights, things are + // equally simple: copy them and + // scale them + std::vector weights = quadrature.get_weights(); + for (unsigned int i = 0; i < n_q_points; ++i) + weights[i] *= (1. / GeometryInfo::max_children_per_cell); + + return Quadrature(q_points, weights); +} + + + +template +Quadrature +QProjector::project_to_all_children(const Quadrature &quadrature) +{ + return project_to_all_children(ReferenceCell::get_hypercube(dim), quadrature); +} + + + +template +Quadrature +QProjector::project_to_all_children( + const ReferenceCell::Type reference_cell_type, + const Quadrature & quadrature) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int n_points = quadrature.size(), + n_children = GeometryInfo::max_children_per_cell; + + std::vector> q_points(n_points * n_children); + std::vector weights(n_points * n_children); + + // project to each child and copy + // results + for (unsigned int child = 0; child < n_children; ++child) + { + Quadrature help = project_to_child(quadrature, child); + for (unsigned int i = 0; i < n_points; ++i) + { + q_points[child * n_points + i] = help.point(i); + weights[child * n_points + i] = help.weight(i); + } + } + return Quadrature(q_points, weights); +} + + + +template +Quadrature +QProjector::project_to_line(const Quadrature<1> &quadrature, + const Point & p1, + const Point & p2) +{ + return project_to_line(ReferenceCell::get_hypercube(dim), quadrature, p1, p2); +} + + + +template +Quadrature +QProjector::project_to_line(const ReferenceCell::Type reference_cell_type, + const Quadrature<1> & quadrature, + const Point & p1, + const Point & p2) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + const unsigned int n = quadrature.size(); + std::vector> points(n); + std::vector weights(n); + const double length = p1.distance(p2); + + for (unsigned int k = 0; k < n; ++k) + { + const double alpha = quadrature.point(k)(0); + points[k] = alpha * p2; + points[k] += (1. - alpha) * p1; + weights[k] = length * quadrature.weight(k); + } + return Quadrature(points, weights); +} + + + +template +typename QProjector::DataSetDescriptor +QProjector::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) +{ + return face(ReferenceCell::get_hypercube(dim), + face_no, + face_orientation, + face_flip, + face_rotation, + n_quadrature_points); +} + + + +template +typename QProjector::DataSetDescriptor +QProjector::DataSetDescriptor::face( + const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const unsigned int n_quadrature_points) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + + switch (dim) + { + case 1: + case 2: + return face_no * n_quadrature_points; + + + case 3: + { + // in 3d, we have to account for faces that + // have non-standard face orientation, flip + // and rotation. thus, we have to store + // _eight_ data sets per face or subface + + // set up a table with the according offsets + // for non-standard orientation, first index: + // face_orientation (standard true=1), second + // index: face_flip (standard false=0), third + // index: face_rotation (standard false=0) + // + // note, that normally we should use the + // obvious offsets 0,1,2,3,4,5,6,7. However, + // prior to the changes enabling flipped and + // rotated faces, in many places of the + // library the convention was used, that the + // first dataset with offset 0 corresponds to + // a face in standard orientation. therefore + // we use the offsets 4,5,6,7,0,1,2,3 here to + // stick to that (implicit) convention + static const unsigned int offset[2][2][2] = { + {{4 * GeometryInfo::faces_per_cell, + 5 * GeometryInfo:: + faces_per_cell}, // face_orientation=false; face_flip=false; + // face_rotation=false and true + {6 * GeometryInfo::faces_per_cell, + 7 * GeometryInfo:: + faces_per_cell}}, // face_orientation=false; face_flip=true; + // face_rotation=false and true + {{0 * GeometryInfo::faces_per_cell, + 1 * GeometryInfo:: + faces_per_cell}, // face_orientation=true; face_flip=false; + // face_rotation=false and true + {2 * GeometryInfo::faces_per_cell, + 3 * GeometryInfo:: + faces_per_cell}}}; // face_orientation=true; face_flip=true; + // face_rotation=false and true + + return ( + (face_no + offset[face_orientation][face_flip][face_rotation]) * + n_quadrature_points); + } + + default: + Assert(false, ExcInternalError()); + } + return numbers::invalid_unsigned_int; +} + + + +template <> +QProjector<1>::DataSetDescriptor +QProjector<1>::DataSetDescriptor::subface( + const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const unsigned int subface_no, + const bool, + const bool, + const bool, + const unsigned int n_quadrature_points, + const internal::SubfaceCase<1>) +{ + Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); + (void)reference_cell_type; + + Assert(face_no < GeometryInfo<1>::faces_per_cell, ExcInternalError()); + Assert(subface_no < GeometryInfo<1>::max_children_per_face, + ExcInternalError()); + + return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) * + n_quadrature_points); +} + + + +template <> +QProjector<1>::DataSetDescriptor +QProjector<1>::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<1> ref_case) +{ + return subface(ReferenceCell::Type::Line, + face_no, + subface_no, + face_orientation, + face_flip, + face_rotation, + n_quadrature_points, + ref_case); +} + + + +template <> +QProjector<2>::DataSetDescriptor +QProjector<2>::DataSetDescriptor::subface( + const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const unsigned int subface_no, + const bool, + const bool, + const bool, + const unsigned int n_quadrature_points, + const internal::SubfaceCase<2>) +{ + Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); + (void)reference_cell_type; + + Assert(face_no < GeometryInfo<2>::faces_per_cell, ExcInternalError()); + Assert(subface_no < GeometryInfo<2>::max_children_per_face, + ExcInternalError()); + + return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) * + n_quadrature_points); +} + + + +template <> +QProjector<2>::DataSetDescriptor +QProjector<2>::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<2> ref_case) +{ + return subface(ReferenceCell::Type::Quad, + face_no, + subface_no, + face_orientation, + face_flip, + face_rotation, + n_quadrature_points, + ref_case); +} + + +template <> +QProjector<3>::DataSetDescriptor +QProjector<3>::DataSetDescriptor::subface( + const ReferenceCell::Type reference_cell_type, + 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<3> ref_case) +{ + const unsigned int dim = 3; + + Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented()); + (void)reference_cell_type; + + Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + Assert(subface_no < GeometryInfo::max_children_per_face, + ExcInternalError()); + + // As the quadrature points created by + // QProjector are on subfaces in their + // "standard location" we have to use a + // permutation of the equivalent subface + // number in order to respect face + // orientation, flip and rotation. The + // information we need here is exactly the + // same as the + // GeometryInfo<3>::child_cell_on_face info + // for the bottom face (face 4) of a hex, as + // on this the RefineCase of the cell matches + // that of the face and the subfaces are + // numbered in the same way as the child + // cells. + + // in 3d, we have to account for faces that + // have non-standard face orientation, flip + // and rotation. thus, we have to store + // _eight_ data sets per face or subface + // already for the isotropic + // case. Additionally, we have three + // different refinement cases, resulting in + // 4 + 2 + 2 = 8 different subfaces + // for each face. + const unsigned int total_subfaces_per_face = 8; + + // set up a table with the according offsets + // for non-standard orientation, first index: + // face_orientation (standard true=1), second + // index: face_flip (standard false=0), third + // index: face_rotation (standard false=0) + // + // note, that normally we should use the + // obvious offsets 0,1,2,3,4,5,6,7. However, + // prior to the changes enabling flipped and + // rotated faces, in many places of the + // library the convention was used, that the + // first dataset with offset 0 corresponds to + // a face in standard orientation. therefore + // we use the offsets 4,5,6,7,0,1,2,3 here to + // stick to that (implicit) convention + static const unsigned int orientation_offset[2][2][2] = { + {// face_orientation=false; face_flip=false; face_rotation=false and true + {4 * GeometryInfo::faces_per_cell * total_subfaces_per_face, + 5 * GeometryInfo::faces_per_cell * total_subfaces_per_face}, + // face_orientation=false; face_flip=true; face_rotation=false and true + {6 * GeometryInfo::faces_per_cell * total_subfaces_per_face, + 7 * GeometryInfo::faces_per_cell * total_subfaces_per_face}}, + {// face_orientation=true; face_flip=false; face_rotation=false and true + {0 * GeometryInfo::faces_per_cell * total_subfaces_per_face, + 1 * GeometryInfo::faces_per_cell * total_subfaces_per_face}, + // face_orientation=true; face_flip=true; face_rotation=false and true + {2 * GeometryInfo::faces_per_cell * total_subfaces_per_face, + 3 * GeometryInfo::faces_per_cell * total_subfaces_per_face}}}; + + // set up a table with the offsets for a + // given refinement case respecting the + // corresponding number of subfaces. the + // index corresponds to (RefineCase::Type - 1) + + // note, that normally we should use the + // obvious offsets 0,2,6. However, prior to + // the implementation of anisotropic + // refinement, in many places of the library + // the convention was used, that the first + // dataset with offset 0 corresponds to a + // standard (isotropic) face + // refinement. therefore we use the offsets + // 6,4,0 here to stick to that (implicit) + // convention + static const unsigned int ref_case_offset[3] = { + 6, // cut_x + 4, // cut_y + 0 // cut_xy + }; + + + // for each subface of a given FaceRefineCase + // there is a corresponding equivalent + // subface number of one of the "standard" + // RefineCases (cut_x, cut_y, cut_xy). Map + // the given values to those equivalent + // ones. + + // first, define an invalid number + static const unsigned int e = numbers::invalid_unsigned_int; + + static const RefinementCase + equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] + [GeometryInfo<3>::max_children_per_face] = { + // case_none. there should be only + // invalid values here. However, as + // this function is also called (in + // tests) for cells which have no + // refined faces, use isotropic + // refinement instead + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_x + {RefinementCase::cut_x, + RefinementCase::cut_x, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_x1y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_x, + RefinementCase::no_refinement}, + // case_x2y + {RefinementCase::cut_x, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_x1y2y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_y + {RefinementCase::cut_y, + RefinementCase::cut_y, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_y1x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_y, + RefinementCase::no_refinement}, + // case_y2x + {RefinementCase::cut_y, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_y1x2x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_xy (case_isotropic) + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}}; + + static const unsigned int + equivalent_subface_number[internal::SubfaceCase::case_isotropic + 1] + [GeometryInfo<3>::max_children_per_face] = { + // case_none, see above + {0, 1, 2, 3}, + // case_x + {0, 1, e, e}, + // case_x1y + {0, 2, 1, e}, + // case_x2y + {0, 1, 3, e}, + // case_x1y2y + {0, 2, 1, 3}, + // case_y + {0, 1, e, e}, + // case_y1x + {0, 1, 1, e}, + // case_y2x + {0, 2, 3, e}, + // case_y1x2x + {0, 1, 2, 3}, + // case_xy (case_isotropic) + {0, 1, 2, 3}}; + + // If face-orientation or face_rotation are + // non-standard, cut_x and cut_y have to be + // exchanged. + static const RefinementCase ref_case_permutation[4] = { + RefinementCase::no_refinement, + RefinementCase::cut_y, + RefinementCase::cut_x, + RefinementCase::cut_xy}; + + // set a corresponding (equivalent) + // RefineCase and subface number + const RefinementCase equ_ref_case = + equivalent_refine_case[ref_case][subface_no]; + const unsigned int equ_subface_no = + equivalent_subface_number[ref_case][subface_no]; + // make sure, that we got a valid subface and RefineCase + Assert(equ_ref_case != RefinementCase::no_refinement, + ExcInternalError()); + Assert(equ_subface_no != e, ExcInternalError()); + // now, finally respect non-standard faces + const RefinementCase final_ref_case = + (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] : + equ_ref_case); + + // what we have now is the number of + // the subface in the natural + // orientation of the *face*. what we + // need to know is the number of the + // subface concerning the standard face + // orientation as seen from the *cell*. + + // this mapping is not trivial, but we + // have done exactly this stuff in the + // child_cell_on_face function. in + // order to reduce the amount of code + // as well as to make maintaining the + // functionality easier we want to + // reuse that information. So we note + // that on the bottom face (face 4) of + // a hex cell the local x and y + // coordinates of the face and the cell + // coincide, thus also the refinement + // case of the face corresponds to the + // refinement case of the cell + // (ignoring cell refinement along the + // z direction). Using this knowledge + // we can (ab)use the + // child_cell_on_face function to do + // exactly the transformation we are in + // need of now + const unsigned int final_subface_no = + GeometryInfo::child_cell_on_face(RefinementCase(final_ref_case), + 4, + equ_subface_no, + face_orientation, + face_flip, + face_rotation, + equ_ref_case); + + return (((face_no * total_subfaces_per_face + + ref_case_offset[final_ref_case - 1] + final_subface_no) + + orientation_offset[face_orientation][face_flip][face_rotation]) * + n_quadrature_points); +} + + +template <> +QProjector<3>::DataSetDescriptor +QProjector<3>::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<3> ref_case) +{ + return subface(ReferenceCell::Type::Hex, + face_no, + subface_no, + face_orientation, + face_flip, + face_rotation, + n_quadrature_points, + ref_case); +} + + + +template +Quadrature +QProjector::project_to_face(const SubQuadrature &quadrature, + const unsigned int face_no) +{ + return project_to_face(ReferenceCell::get_hypercube(dim), + quadrature, + face_no); +} + + + +template +Quadrature +QProjector::project_to_face(const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + std::vector> points(quadrature.size()); + project_to_face(quadrature, face_no, points); + return Quadrature(points, quadrature.get_weights()); +} + + + +template +Quadrature +QProjector::project_to_subface(const SubQuadrature & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const RefinementCase &ref_case) +{ + return project_to_subface(ReferenceCell::get_hypercube(dim), + quadrature, + face_no, + subface_no, + ref_case); +} + + + +template +Quadrature +QProjector::project_to_subface( + const ReferenceCell::Type reference_cell_type, + const SubQuadrature & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const RefinementCase &ref_case) +{ + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + (void)reference_cell_type; + + std::vector> points(quadrature.size()); + project_to_subface(quadrature, face_no, subface_no, points, ref_case); + return Quadrature(points, quadrature.get_weights()); +} + + +// explicit instantiations; note: we need them all for all dimensions +template class QProjector<1>; +template class QProjector<2>; +template class QProjector<3>; + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/quadrature.cc b/source/base/quadrature.cc index 25829085e2..9e79b4f987 100644 --- a/source/base/quadrature.cc +++ b/source/base/quadrature.cc @@ -424,1105 +424,6 @@ QAnisotropic<3>::QAnisotropic(const Quadrature<1> &qx, -//--------------------------------------------------------------------------- - - - -template -Quadrature<2> -QProjector::reflect(const Quadrature<2> &q) -{ - std::vector> q_points(q.size()); - std::vector weights(q.size()); - for (unsigned int i = 0; i < q.size(); ++i) - { - q_points[i][0] = q.point(i)[1]; - q_points[i][1] = q.point(i)[0]; - - weights[i] = q.weight(i); - } - - return Quadrature<2>(q_points, weights); -} - - -template -Quadrature<2> -QProjector::rotate(const Quadrature<2> &q, const unsigned int n_times) -{ - std::vector> q_points(q.size()); - std::vector weights(q.size()); - for (unsigned int i = 0; i < q.size(); ++i) - { - switch (n_times % 4) - { - case 0: - // 0 degree - q_points[i][0] = q.point(i)[0]; - q_points[i][1] = q.point(i)[1]; - break; - case 1: - // 90 degree counterclockwise - q_points[i][0] = 1.0 - q.point(i)[1]; - q_points[i][1] = q.point(i)[0]; - break; - case 2: - // 180 degree counterclockwise - q_points[i][0] = 1.0 - q.point(i)[0]; - q_points[i][1] = 1.0 - q.point(i)[1]; - break; - case 3: - // 270 degree counterclockwise - q_points[i][0] = q.point(i)[1]; - q_points[i][1] = 1.0 - q.point(i)[0]; - break; - } - - weights[i] = q.weight(i); - } - - return Quadrature<2>(q_points, weights); -} - - -template <> -void -QProjector<1>::project_to_face(const Quadrature<0> &, - const unsigned int face_no, - std::vector> &q_points) -{ - const unsigned int dim = 1; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - AssertDimension(q_points.size(), 1); - - q_points[0] = Point(static_cast(face_no)); -} - - - -template <> -void -QProjector<2>::project_to_face(const Quadrature<1> & quadrature, - const unsigned int face_no, - std::vector> &q_points) -{ - const unsigned int dim = 2; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - Assert(q_points.size() == quadrature.size(), - ExcDimensionMismatch(q_points.size(), quadrature.size())); - - for (unsigned int p = 0; p < quadrature.size(); ++p) - switch (face_no) - { - case 0: - q_points[p] = Point(0, quadrature.point(p)(0)); - break; - case 1: - q_points[p] = Point(1, quadrature.point(p)(0)); - break; - case 2: - q_points[p] = Point(quadrature.point(p)(0), 0); - break; - case 3: - q_points[p] = Point(quadrature.point(p)(0), 1); - break; - default: - Assert(false, ExcInternalError()); - } -} - - - -template <> -void -QProjector<3>::project_to_face(const Quadrature<2> & quadrature, - const unsigned int face_no, - std::vector> &q_points) -{ - const unsigned int dim = 3; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - Assert(q_points.size() == quadrature.size(), - ExcDimensionMismatch(q_points.size(), quadrature.size())); - - for (unsigned int p = 0; p < quadrature.size(); ++p) - switch (face_no) - { - case 0: - q_points[p] = - Point(0, quadrature.point(p)(0), quadrature.point(p)(1)); - break; - case 1: - q_points[p] = - Point(1, quadrature.point(p)(0), quadrature.point(p)(1)); - break; - case 2: - q_points[p] = - Point(quadrature.point(p)(1), 0, quadrature.point(p)(0)); - break; - case 3: - q_points[p] = - Point(quadrature.point(p)(1), 1, quadrature.point(p)(0)); - break; - case 4: - q_points[p] = - Point(quadrature.point(p)(0), quadrature.point(p)(1), 0); - break; - case 5: - q_points[p] = - Point(quadrature.point(p)(0), quadrature.point(p)(1), 1); - break; - - default: - Assert(false, ExcInternalError()); - } -} - - - -template <> -void -QProjector<1>::project_to_subface(const Quadrature<0> &, - const unsigned int face_no, - const unsigned int, - std::vector> &q_points, - const RefinementCase<0> &) -{ - const unsigned int dim = 1; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - AssertDimension(q_points.size(), 1); - - q_points[0] = Point(static_cast(face_no)); -} - - - -template <> -void -QProjector<2>::project_to_subface(const Quadrature<1> & quadrature, - const unsigned int face_no, - const unsigned int subface_no, - std::vector> &q_points, - const RefinementCase<1> &) -{ - const unsigned int dim = 2; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - AssertIndexRange(subface_no, GeometryInfo::max_children_per_face); - - Assert(q_points.size() == quadrature.size(), - ExcDimensionMismatch(q_points.size(), quadrature.size())); - - for (unsigned int p = 0; p < quadrature.size(); ++p) - switch (face_no) - { - case 0: - switch (subface_no) - { - case 0: - q_points[p] = Point(0, quadrature.point(p)(0) / 2); - break; - case 1: - q_points[p] = Point(0, quadrature.point(p)(0) / 2 + 0.5); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 1: - switch (subface_no) - { - case 0: - q_points[p] = Point(1, quadrature.point(p)(0) / 2); - break; - case 1: - q_points[p] = Point(1, quadrature.point(p)(0) / 2 + 0.5); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 2: - switch (subface_no) - { - case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 0); - break; - case 1: - q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 0); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 3: - switch (subface_no) - { - case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 1); - break; - case 1: - q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 1); - break; - default: - Assert(false, ExcInternalError()); - } - break; - - default: - Assert(false, ExcInternalError()); - } -} - - - -template <> -void -QProjector<3>::project_to_subface(const Quadrature<2> & quadrature, - const unsigned int face_no, - const unsigned int subface_no, - std::vector> & q_points, - const RefinementCase<2> &ref_case) -{ - const unsigned int dim = 3; - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); - AssertIndexRange(subface_no, GeometryInfo::max_children_per_face); - Assert(q_points.size() == quadrature.size(), - ExcDimensionMismatch(q_points.size(), quadrature.size())); - - // one coordinate is at a const value. for - // faces 0, 2 and 4 this value is 0.0, for - // faces 1, 3 and 5 it is 1.0 - double const_value = face_no % 2; - // local 2d coordinates are xi and eta, - // global 3d coordinates are x, y and - // z. those have to be mapped. the following - // indices tell, which global coordinate - // (0->x, 1->y, 2->z) corresponds to which - // local one - unsigned int xi_index = numbers::invalid_unsigned_int, - eta_index = numbers::invalid_unsigned_int, - const_index = face_no / 2; - // the xi and eta values have to be scaled - // (by factor 0.5 or factor 1.0) depending on - // the refinement case and translated (by 0.0 - // or 0.5) depending on the refinement case - // and subface_no. - double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0, - eta_translation = 0.0; - // set the index mapping between local and - // global coordinates - switch (face_no / 2) - { - case 0: - xi_index = 1; - eta_index = 2; - break; - case 1: - xi_index = 2; - eta_index = 0; - break; - case 2: - xi_index = 0; - eta_index = 1; - break; - } - // set the scale and translation parameter - // for individual subfaces - switch (ref_case) - { - case RefinementCase::cut_x: - xi_scale = 0.5; - xi_translation = subface_no % 2 * 0.5; - break; - case RefinementCase::cut_y: - eta_scale = 0.5; - eta_translation = subface_no % 2 * 0.5; - break; - case RefinementCase::cut_xy: - xi_scale = 0.5; - eta_scale = 0.5; - xi_translation = int(subface_no % 2) * 0.5; - eta_translation = int(subface_no / 2) * 0.5; - break; - default: - Assert(false, ExcInternalError()); - break; - } - // finally, compute the scaled, translated, - // projected quadrature points - for (unsigned int p = 0; p < quadrature.size(); ++p) - { - q_points[p][xi_index] = - xi_scale * quadrature.point(p)(0) + xi_translation; - q_points[p][eta_index] = - eta_scale * quadrature.point(p)(1) + eta_translation; - q_points[p][const_index] = const_value; - } -} - - -template <> -Quadrature<1> -QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature) -{ - const unsigned int dim = 1; - - const unsigned int n_points = 1, n_faces = GeometryInfo::faces_per_cell; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces); - std::vector> help(n_points); - - - // project to each face and append - // results - for (unsigned int face = 0; face < n_faces; ++face) - { - project_to_face(quadrature, face, help); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - std::vector weights; - weights.reserve(n_points * n_faces); - for (unsigned int face = 0; face < n_faces; ++face) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); - - Assert(q_points.size() == n_points * n_faces, ExcInternalError()); - Assert(weights.size() == n_points * n_faces, ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -template <> -Quadrature<2> -QProjector<2>::project_to_all_faces(const SubQuadrature &quadrature) -{ - const unsigned int dim = 2; - - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces); - std::vector> help(n_points); - - // project to each face and append - // results - for (unsigned int face = 0; face < n_faces; ++face) - { - project_to_face(quadrature, face, help); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - std::vector weights; - weights.reserve(n_points * n_faces); - for (unsigned int face = 0; face < n_faces; ++face) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); - - Assert(q_points.size() == n_points * n_faces, ExcInternalError()); - Assert(weights.size() == n_points * n_faces, ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -template <> -Quadrature<3> -QProjector<3>::project_to_all_faces(const SubQuadrature &quadrature) -{ - const unsigned int dim = 3; - - SubQuadrature q_reflected = reflect(quadrature); - SubQuadrature q[8] = {quadrature, - rotate(quadrature, 1), - rotate(quadrature, 2), - rotate(quadrature, 3), - q_reflected, - rotate(q_reflected, 3), - rotate(q_reflected, 2), - rotate(q_reflected, 1)}; - - - - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces * 8); - std::vector> help(n_points); - - std::vector weights; - weights.reserve(n_points * n_faces * 8); - - // do the following for all possible - // mutations of a face (mutation==0 - // corresponds to a face with standard - // orientation, no flip and no rotation) - for (const auto &mutation : q) - { - // project to each face and append - // results - for (unsigned int face = 0; face < n_faces; ++face) - { - project_to_face(mutation, face, help); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - for (unsigned int face = 0; face < n_faces; ++face) - std::copy(mutation.get_weights().begin(), - mutation.get_weights().end(), - std::back_inserter(weights)); - } - - - Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError()); - Assert(weights.size() == n_points * n_faces * 8, ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -template <> -Quadrature<1> -QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature) -{ - const unsigned int dim = 1; - - const unsigned int n_points = 1, n_faces = GeometryInfo::faces_per_cell, - subfaces_per_face = - GeometryInfo::max_children_per_face; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces * subfaces_per_face); - std::vector> help(n_points); - - // project to each face and copy - // results - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) - { - project_to_subface(quadrature, face, subface, help); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - std::vector weights; - weights.reserve(n_points * n_faces * subfaces_per_face); - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); - - Assert(q_points.size() == n_points * n_faces * subfaces_per_face, - ExcInternalError()); - Assert(weights.size() == n_points * n_faces * subfaces_per_face, - ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -template <> -Quadrature<2> -QProjector<2>::project_to_all_subfaces(const SubQuadrature &quadrature) -{ - const unsigned int dim = 2; - - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell, - subfaces_per_face = - GeometryInfo::max_children_per_face; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces * subfaces_per_face); - std::vector> help(n_points); - - // project to each face and copy - // results - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) - { - project_to_subface(quadrature, face, subface, help); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - std::vector weights; - weights.reserve(n_points * n_faces * subfaces_per_face); - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int subface = 0; subface < subfaces_per_face; ++subface) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); - - Assert(q_points.size() == n_points * n_faces * subfaces_per_face, - ExcInternalError()); - Assert(weights.size() == n_points * n_faces * subfaces_per_face, - ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -template <> -Quadrature<3> -QProjector<3>::project_to_all_subfaces(const SubQuadrature &quadrature) -{ - const unsigned int dim = 3; - SubQuadrature q_reflected = reflect(quadrature); - SubQuadrature q[8] = {quadrature, - rotate(quadrature, 1), - rotate(quadrature, 2), - rotate(quadrature, 3), - q_reflected, - rotate(q_reflected, 3), - rotate(q_reflected, 2), - rotate(q_reflected, 1)}; - - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell, - total_subfaces_per_face = 2 + 2 + 4; - - // first fix quadrature points - std::vector> q_points; - q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8); - std::vector> help(n_points); - - std::vector weights; - weights.reserve(n_points * n_faces * total_subfaces_per_face * 8); - - // do the following for all possible - // mutations of a face (mutation==0 - // corresponds to a face with standard - // orientation, no flip and no rotation) - for (const auto &mutation : q) - { - // project to each face and copy - // results - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int ref_case = RefinementCase::cut_xy; - ref_case >= RefinementCase::cut_x; - --ref_case) - for (unsigned int subface = 0; - subface < GeometryInfo::n_children( - RefinementCase(ref_case)); - ++subface) - { - project_to_subface(mutation, - face, - subface, - help, - RefinementCase(ref_case)); - std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - - // next copy over weights - for (unsigned int face = 0; face < n_faces; ++face) - for (unsigned int ref_case = RefinementCase::cut_xy; - ref_case >= RefinementCase::cut_x; - --ref_case) - for (unsigned int subface = 0; - subface < GeometryInfo::n_children( - RefinementCase(ref_case)); - ++subface) - std::copy(mutation.get_weights().begin(), - mutation.get_weights().end(), - std::back_inserter(weights)); - } - - Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8, - ExcInternalError()); - Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8, - ExcInternalError()); - - return Quadrature(q_points, weights); -} - - - -// This function is not used in the library -template -Quadrature -QProjector::project_to_child(const Quadrature &quadrature, - const unsigned int child_no) -{ - AssertIndexRange(child_no, GeometryInfo::max_children_per_cell); - - const unsigned int n_q_points = quadrature.size(); - - std::vector> q_points(n_q_points); - for (unsigned int i = 0; i < n_q_points; ++i) - q_points[i] = - GeometryInfo::child_to_cell_coordinates(quadrature.point(i), - child_no); - - // for the weights, things are - // equally simple: copy them and - // scale them - std::vector weights = quadrature.get_weights(); - for (unsigned int i = 0; i < n_q_points; ++i) - weights[i] *= (1. / GeometryInfo::max_children_per_cell); - - return Quadrature(q_points, weights); -} - - -template -Quadrature -QProjector::project_to_all_children(const Quadrature &quadrature) -{ - const unsigned int n_points = quadrature.size(), - n_children = GeometryInfo::max_children_per_cell; - - std::vector> q_points(n_points * n_children); - std::vector weights(n_points * n_children); - - // project to each child and copy - // results - for (unsigned int child = 0; child < n_children; ++child) - { - Quadrature help = project_to_child(quadrature, child); - for (unsigned int i = 0; i < n_points; ++i) - { - q_points[child * n_points + i] = help.point(i); - weights[child * n_points + i] = help.weight(i); - } - } - return Quadrature(q_points, weights); -} - - - -template -Quadrature -QProjector::project_to_line(const Quadrature<1> &quadrature, - const Point & p1, - const Point & p2) -{ - const unsigned int n = quadrature.size(); - std::vector> points(n); - std::vector weights(n); - const double length = p1.distance(p2); - - for (unsigned int k = 0; k < n; ++k) - { - const double alpha = quadrature.point(k)(0); - points[k] = alpha * p2; - points[k] += (1. - alpha) * p1; - weights[k] = length * quadrature.weight(k); - } - return Quadrature(points, weights); -} - - - -template -typename QProjector::DataSetDescriptor -QProjector::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) -{ - Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); - - switch (dim) - { - case 1: - case 2: - return face_no * n_quadrature_points; - - - case 3: - { - // in 3d, we have to account for faces that - // have non-standard face orientation, flip - // and rotation. thus, we have to store - // _eight_ data sets per face or subface - - // set up a table with the according offsets - // for non-standard orientation, first index: - // face_orientation (standard true=1), second - // index: face_flip (standard false=0), third - // index: face_rotation (standard false=0) - // - // note, that normally we should use the - // obvious offsets 0,1,2,3,4,5,6,7. However, - // prior to the changes enabling flipped and - // rotated faces, in many places of the - // library the convention was used, that the - // first dataset with offset 0 corresponds to - // a face in standard orientation. therefore - // we use the offsets 4,5,6,7,0,1,2,3 here to - // stick to that (implicit) convention - static const unsigned int offset[2][2][2] = { - {{4 * GeometryInfo::faces_per_cell, - 5 * GeometryInfo:: - faces_per_cell}, // face_orientation=false; face_flip=false; - // face_rotation=false and true - {6 * GeometryInfo::faces_per_cell, - 7 * GeometryInfo:: - faces_per_cell}}, // face_orientation=false; face_flip=true; - // face_rotation=false and true - {{0 * GeometryInfo::faces_per_cell, - 1 * GeometryInfo:: - faces_per_cell}, // face_orientation=true; face_flip=false; - // face_rotation=false and true - {2 * GeometryInfo::faces_per_cell, - 3 * GeometryInfo:: - faces_per_cell}}}; // face_orientation=true; face_flip=true; - // face_rotation=false and true - - return ( - (face_no + offset[face_orientation][face_flip][face_rotation]) * - n_quadrature_points); - } - - default: - Assert(false, ExcInternalError()); - } - return numbers::invalid_unsigned_int; -} - - - -template <> -QProjector<1>::DataSetDescriptor -QProjector<1>::DataSetDescriptor::subface( - const unsigned int face_no, - const unsigned int subface_no, - const bool, - const bool, - const bool, - const unsigned int n_quadrature_points, - const internal::SubfaceCase<1>) -{ - Assert(face_no < GeometryInfo<1>::faces_per_cell, ExcInternalError()); - Assert(subface_no < GeometryInfo<1>::max_children_per_face, - ExcInternalError()); - - return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) * - n_quadrature_points); -} - - - -template <> -QProjector<2>::DataSetDescriptor -QProjector<2>::DataSetDescriptor::subface( - const unsigned int face_no, - const unsigned int subface_no, - const bool, - const bool, - const bool, - const unsigned int n_quadrature_points, - const internal::SubfaceCase<2>) -{ - Assert(face_no < GeometryInfo<2>::faces_per_cell, ExcInternalError()); - Assert(subface_no < GeometryInfo<2>::max_children_per_face, - ExcInternalError()); - - return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) * - n_quadrature_points); -} - - -template <> -QProjector<3>::DataSetDescriptor -QProjector<3>::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<3> ref_case) -{ - const unsigned int dim = 3; - - Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); - Assert(subface_no < GeometryInfo::max_children_per_face, - ExcInternalError()); - - // As the quadrature points created by - // QProjector are on subfaces in their - // "standard location" we have to use a - // permutation of the equivalent subface - // number in order to respect face - // orientation, flip and rotation. The - // information we need here is exactly the - // same as the - // GeometryInfo<3>::child_cell_on_face info - // for the bottom face (face 4) of a hex, as - // on this the RefineCase of the cell matches - // that of the face and the subfaces are - // numbered in the same way as the child - // cells. - - // in 3d, we have to account for faces that - // have non-standard face orientation, flip - // and rotation. thus, we have to store - // _eight_ data sets per face or subface - // already for the isotropic - // case. Additionally, we have three - // different refinement cases, resulting in - // 4 + 2 + 2 = 8 different subfaces - // for each face. - const unsigned int total_subfaces_per_face = 8; - - // set up a table with the according offsets - // for non-standard orientation, first index: - // face_orientation (standard true=1), second - // index: face_flip (standard false=0), third - // index: face_rotation (standard false=0) - // - // note, that normally we should use the - // obvious offsets 0,1,2,3,4,5,6,7. However, - // prior to the changes enabling flipped and - // rotated faces, in many places of the - // library the convention was used, that the - // first dataset with offset 0 corresponds to - // a face in standard orientation. therefore - // we use the offsets 4,5,6,7,0,1,2,3 here to - // stick to that (implicit) convention - static const unsigned int orientation_offset[2][2][2] = { - {// face_orientation=false; face_flip=false; face_rotation=false and true - {4 * GeometryInfo::faces_per_cell * total_subfaces_per_face, - 5 * GeometryInfo::faces_per_cell * total_subfaces_per_face}, - // face_orientation=false; face_flip=true; face_rotation=false and true - {6 * GeometryInfo::faces_per_cell * total_subfaces_per_face, - 7 * GeometryInfo::faces_per_cell * total_subfaces_per_face}}, - {// face_orientation=true; face_flip=false; face_rotation=false and true - {0 * GeometryInfo::faces_per_cell * total_subfaces_per_face, - 1 * GeometryInfo::faces_per_cell * total_subfaces_per_face}, - // face_orientation=true; face_flip=true; face_rotation=false and true - {2 * GeometryInfo::faces_per_cell * total_subfaces_per_face, - 3 * GeometryInfo::faces_per_cell * total_subfaces_per_face}}}; - - // set up a table with the offsets for a - // given refinement case respecting the - // corresponding number of subfaces. the - // index corresponds to (RefineCase::Type - 1) - - // note, that normally we should use the - // obvious offsets 0,2,6. However, prior to - // the implementation of anisotropic - // refinement, in many places of the library - // the convention was used, that the first - // dataset with offset 0 corresponds to a - // standard (isotropic) face - // refinement. therefore we use the offsets - // 6,4,0 here to stick to that (implicit) - // convention - static const unsigned int ref_case_offset[3] = { - 6, // cut_x - 4, // cut_y - 0 // cut_xy - }; - - - // for each subface of a given FaceRefineCase - // there is a corresponding equivalent - // subface number of one of the "standard" - // RefineCases (cut_x, cut_y, cut_xy). Map - // the given values to those equivalent - // ones. - - // first, define an invalid number - static const unsigned int e = numbers::invalid_unsigned_int; - - static const RefinementCase - equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] - [GeometryInfo<3>::max_children_per_face] = { - // case_none. there should be only - // invalid values here. However, as - // this function is also called (in - // tests) for cells which have no - // refined faces, use isotropic - // refinement instead - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_x - {RefinementCase::cut_x, - RefinementCase::cut_x, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_x1y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_x, - RefinementCase::no_refinement}, - // case_x2y - {RefinementCase::cut_x, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_x1y2y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_y - {RefinementCase::cut_y, - RefinementCase::cut_y, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_y1x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_y, - RefinementCase::no_refinement}, - // case_y2x - {RefinementCase::cut_y, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_y1x2x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_xy (case_isotropic) - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}}; - - static const unsigned int - equivalent_subface_number[internal::SubfaceCase::case_isotropic + 1] - [GeometryInfo<3>::max_children_per_face] = { - // case_none, see above - {0, 1, 2, 3}, - // case_x - {0, 1, e, e}, - // case_x1y - {0, 2, 1, e}, - // case_x2y - {0, 1, 3, e}, - // case_x1y2y - {0, 2, 1, 3}, - // case_y - {0, 1, e, e}, - // case_y1x - {0, 1, 1, e}, - // case_y2x - {0, 2, 3, e}, - // case_y1x2x - {0, 1, 2, 3}, - // case_xy (case_isotropic) - {0, 1, 2, 3}}; - - // If face-orientation or face_rotation are - // non-standard, cut_x and cut_y have to be - // exchanged. - static const RefinementCase ref_case_permutation[4] = { - RefinementCase::no_refinement, - RefinementCase::cut_y, - RefinementCase::cut_x, - RefinementCase::cut_xy}; - - // set a corresponding (equivalent) - // RefineCase and subface number - const RefinementCase equ_ref_case = - equivalent_refine_case[ref_case][subface_no]; - const unsigned int equ_subface_no = - equivalent_subface_number[ref_case][subface_no]; - // make sure, that we got a valid subface and RefineCase - Assert(equ_ref_case != RefinementCase::no_refinement, - ExcInternalError()); - Assert(equ_subface_no != e, ExcInternalError()); - // now, finally respect non-standard faces - const RefinementCase final_ref_case = - (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] : - equ_ref_case); - - // what we have now is the number of - // the subface in the natural - // orientation of the *face*. what we - // need to know is the number of the - // subface concerning the standard face - // orientation as seen from the *cell*. - - // this mapping is not trivial, but we - // have done exactly this stuff in the - // child_cell_on_face function. in - // order to reduce the amount of code - // as well as to make maintaining the - // functionality easier we want to - // reuse that information. So we note - // that on the bottom face (face 4) of - // a hex cell the local x and y - // coordinates of the face and the cell - // coincide, thus also the refinement - // case of the face corresponds to the - // refinement case of the cell - // (ignoring cell refinement along the - // z direction). Using this knowledge - // we can (ab)use the - // child_cell_on_face function to do - // exactly the transformation we are in - // need of now - const unsigned int final_subface_no = - GeometryInfo::child_cell_on_face(RefinementCase(final_ref_case), - 4, - equ_subface_no, - face_orientation, - face_flip, - face_rotation, - equ_ref_case); - - return (((face_no * total_subfaces_per_face + - ref_case_offset[final_ref_case - 1] + final_subface_no) + - orientation_offset[face_orientation][face_flip][face_rotation]) * - n_quadrature_points); -} - - - -template -Quadrature -QProjector::project_to_face(const SubQuadrature &quadrature, - const unsigned int face_no) -{ - std::vector> points(quadrature.size()); - project_to_face(quadrature, face_no, points); - return Quadrature(points, quadrature.get_weights()); -} - - -template -Quadrature -QProjector::project_to_subface(const SubQuadrature & quadrature, - const unsigned int face_no, - const unsigned int subface_no, - const RefinementCase &ref_case) -{ - std::vector> points(quadrature.size()); - project_to_subface(quadrature, face_no, subface_no, points, ref_case); - return Quadrature(points, quadrature.get_weights()); -} - - // ------------------------------------------------------------ // namespace internal @@ -1792,8 +693,5 @@ template class QAnisotropic<3>; template class QIterated<1>; template class QIterated<2>; template class QIterated<3>; -template class QProjector<1>; -template class QProjector<2>; -template class QProjector<3>; DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe.cc b/source/fe/fe.cc index f95944c77a..abf4015204 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -1218,7 +1218,8 @@ FiniteElement::get_face_data( { return get_data(flags, mapping, - QProjector::project_to_all_faces(quadrature), + QProjector::project_to_all_faces( + this->reference_cell_type(), quadrature), output_data); } @@ -1236,7 +1237,8 @@ FiniteElement::get_subface_data( { return get_data(flags, mapping, - QProjector::project_to_all_subfaces(quadrature), + QProjector::project_to_all_subfaces( + this->reference_cell_type(), quadrature), output_data); } diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 432f891806..83732630e5 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -207,7 +207,8 @@ FE_ABF::initialize_support_points(const unsigned int deg) } Quadrature faces = - QProjector::project_to_all_faces(face_points); + QProjector::project_to_all_faces(this->reference_cell_type(), + face_points); for (; current < GeometryInfo::faces_per_cell * n_face_points; ++current) { @@ -339,7 +340,10 @@ FE_ABF::initialize_restriction() // child cell are evaluated // in the quadrature points // of a full face. - Quadrature q_face = QProjector::project_to_face(q_base, face); + Quadrature q_face = + QProjector::project_to_face(this->reference_cell_type(), + q_base, + face); // Store shape values, since the // evaluation suffers if not // ordered by point @@ -357,8 +361,8 @@ FE_ABF::initialize_restriction() // the coarse face are // evaluated on the subface // only. - Quadrature q_sub = - QProjector::project_to_subface(q_base, face, sub); + Quadrature q_sub = QProjector::project_to_subface( + this->reference_cell_type(), q_base, face, sub); const unsigned int child = GeometryInfo::child_cell_on_face( RefinementCase::isotropic_refinement, face, sub); @@ -434,7 +438,10 @@ FE_ABF::initialize_restriction() for (unsigned int child = 0; child < GeometryInfo::max_children_per_cell; ++child) { - Quadrature q_sub = QProjector::project_to_child(q_cell, child); + Quadrature q_sub = + QProjector::project_to_child(this->reference_cell_type(), + q_cell, + child); for (unsigned int k = 0; k < q_sub.size(); ++k) for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell(); @@ -592,7 +599,12 @@ FE_ABF::convert_generalized_support_point_values_to_dof_values( // TODO: Check what the face_orientation, face_flip and face_rotation // have to be in 3D unsigned int k = QProjector::DataSetDescriptor::face( - face, false, false, false, n_face_points); + this->reference_cell_type(), + face, + false, + false, + false, + n_face_points); for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i) nodal_values[start_abf_dofs + i] += n_orient * boundary_weights_abf(k + fp, i) * diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 264203cfb3..a85ede4162 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -349,13 +349,19 @@ FE_BDM::initialize_support_points(const unsigned int deg) this->generalized_support_points.resize(npoints); - Quadrature faces = QProjector::project_to_all_faces(face_points); + Quadrature faces = + QProjector::project_to_all_faces(this->reference_cell_type(), + face_points); for (unsigned int k = 0; k < face_points.size() * GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = - faces.point(k + QProjector::DataSetDescriptor::face( - 0, true, false, false, this->n_dofs_per_face())); + this->generalized_support_points[k] = faces.point( + k + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face())); // Currently, for backward compatibility, we do not use moments, but // point values on faces in 2D. In 3D, this is impossible, since the diff --git a/source/fe/fe_bernstein.cc b/source/fe/fe_bernstein.cc index f8078b5173..587301ecbb 100644 --- a/source/fe/fe_bernstein.cc +++ b/source/fe/fe_bernstein.cc @@ -152,8 +152,13 @@ FE_Bernstein::get_subface_interpolation_matrix( // be done for the face orientation flag in 3D. const Quadrature subface_quadrature = subface == numbers::invalid_unsigned_int ? - QProjector::project_to_face(quad_face_support, 0) : - QProjector::project_to_subface(quad_face_support, 0, subface); + QProjector::project_to_face(this->reference_cell_type(), + quad_face_support, + 0) : + QProjector::project_to_subface(this->reference_cell_type(), + quad_face_support, + 0, + subface); for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i) { diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index a6b9717f79..9645e6a0c9 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -270,7 +270,8 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) const unsigned int n_boundary_points = GeometryInfo::lines_per_cell * n_edge_points; const Quadrature edge_quadrature = - QProjector::project_to_all_faces(reference_edge_quadrature); + QProjector::project_to_all_faces(this->reference_cell_type(), + reference_edge_quadrature); this->generalized_face_support_points.resize(n_edge_points); @@ -297,7 +298,12 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) ++line) this->generalized_support_points[line * n_edge_points + q_point] = edge_quadrature.point(QProjector::DataSetDescriptor::face( - line, true, false, false, n_edge_points) + + this->reference_cell_type(), + line, + true, + false, + false, + n_edge_points) + q_point); for (unsigned int i = 0; i < order; ++i) @@ -323,7 +329,12 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) this->generalized_support_points[line * n_edge_points + q_point] = edge_quadrature.point(QProjector::DataSetDescriptor::face( - line, true, false, false, n_edge_points) + + this->reference_cell_type(), + line, + true, + false, + false, + n_edge_points) + q_point); } } @@ -350,7 +361,9 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) const QGauss<1> reference_edge_quadrature(order + 1); const unsigned int n_edge_points = reference_edge_quadrature.size(); const Quadrature &edge_quadrature = - QProjector::project_to_all_faces(reference_edge_quadrature); + QProjector::project_to_all_faces(ReferenceCell::get_hypercube(dim - + 1), + reference_edge_quadrature); if (order > 0) { @@ -382,7 +395,12 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( - line, true, false, false, n_edge_points) + + ReferenceCell::get_hypercube(dim - 1), + line, + true, + false, + false, + n_edge_points) + q_point); for (unsigned int i = 0; i < 2; ++i) @@ -440,7 +458,8 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) } const Quadrature &face_quadrature = - QProjector::project_to_all_faces(reference_face_quadrature); + QProjector::project_to_all_faces(this->reference_cell_type(), + reference_face_quadrature); for (const unsigned int face : GeometryInfo::face_indices()) for (unsigned int q_point = 0; q_point < n_face_points; ++q_point) @@ -449,7 +468,12 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) GeometryInfo::lines_per_cell * n_edge_points] = face_quadrature.point(QProjector::DataSetDescriptor::face( - face, true, false, false, n_face_points) + + this->reference_cell_type(), + face, + true, + false, + false, + n_face_points) + q_point); } @@ -474,7 +498,12 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( - line, true, false, false, n_edge_points) + + ReferenceCell::get_hypercube(dim - 1), + line, + true, + false, + false, + n_edge_points) + q_point); for (unsigned int i = 0; i < 2; ++i) diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc index ebf33c332b..7a9d9bb646 100644 --- a/source/fe/fe_nedelec_sz.cc +++ b/source/fe/fe_nedelec_sz.cc @@ -2056,19 +2056,22 @@ FE_NedelecSZ::fill_fe_face_values( // This will fill in the missing items in the InternalData // (fe_internal/fe_data) which was not filled in by get_data. fill_edge_values(cell, - QProjector::project_to_all_faces(quadrature), + QProjector::project_to_all_faces( + this->reference_cell_type(), quadrature), fe_data); if (dim == 3 && this->degree > 1) { fill_face_values(cell, - QProjector::project_to_all_faces(quadrature), + QProjector::project_to_all_faces( + this->reference_cell_type(), quadrature), fe_data); } const UpdateFlags flags(fe_data.update_each); const unsigned int n_q_points = quadrature.size(); - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::face(face_no, + const auto offset = + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc index e66652bebc..e86950d1c7 100644 --- a/source/fe/fe_p1nc.cc +++ b/source/fe/fe_p1nc.cc @@ -262,7 +262,9 @@ FE_P1NC::fill_fe_face_values( // compute on the face const Quadrature<2> quadrature_on_face = - QProjector<2>::project_to_face(quadrature, face_no); + QProjector<2>::project_to_face(this->reference_cell_type(), + quadrature, + face_no); if (flags & update_values) for (unsigned int i = 0; i < quadrature_on_face.size(); ++i) @@ -306,8 +308,8 @@ FE_P1NC::fill_fe_subface_values( get_linear_shape_coefficients(cell); // compute on the subface - const Quadrature<2> quadrature_on_subface = - QProjector<2>::project_to_subface(quadrature, face_no, sub_no); + const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface( + this->reference_cell_type(), quadrature, face_no, sub_no); if (flags & update_values) for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i) diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 0f7abca45d..d6f1794bf3 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -969,8 +969,9 @@ FE_PolyTensor::fill_fe_face_values( // to take (all data sets for all // faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::face(face_no, + const auto offset = + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), @@ -1595,8 +1596,9 @@ FE_PolyTensor::fill_fe_subface_values( // offset determines which data set // to take (all data sets for all // sub-faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::subface(face_no, + const auto offset = + QProjector::DataSetDescriptor::subface(this->reference_cell_type(), + face_no, sub_no, cell->face_orientation(face_no), cell->face_flip(face_no), diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 03e69ae8a0..d6664d6ea9 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -260,19 +260,23 @@ struct FE_Q_Base::Implementation // Add nodes of lines interior in the "mother-face" // line 5: use line 9 - QProjector::project_to_subface(qline, 0, 0, p_line); + QProjector::project_to_subface( + ReferenceCell::get_hypercube(dim - 1), qline, 0, 0, p_line); for (unsigned int i = 0; i < n; ++i) constraint_points.push_back(p_line[i] + Point(0.5, 0)); // line 6: use line 10 - QProjector::project_to_subface(qline, 0, 1, p_line); + QProjector::project_to_subface( + ReferenceCell::get_hypercube(dim - 1), qline, 0, 1, p_line); for (unsigned int i = 0; i < n; ++i) constraint_points.push_back(p_line[i] + Point(0.5, 0)); // line 7: use line 13 - QProjector::project_to_subface(qline, 2, 0, p_line); + QProjector::project_to_subface( + ReferenceCell::get_hypercube(dim - 1), qline, 2, 0, p_line); for (unsigned int i = 0; i < n; ++i) constraint_points.push_back(p_line[i] + Point(0, 0.5)); // line 8: use line 14 - QProjector::project_to_subface(qline, 2, 1, p_line); + QProjector::project_to_subface( + ReferenceCell::get_hypercube(dim - 1), qline, 2, 1, p_line); for (unsigned int i = 0; i < n; ++i) constraint_points.push_back(p_line[i] + Point(0, 0.5)); @@ -284,10 +288,12 @@ struct FE_Q_Base::Implementation subface < GeometryInfo::max_children_per_face; ++subface) { - QProjector::project_to_subface(qline, - face, - subface, - p_line); + QProjector::project_to_subface( + ReferenceCell::get_hypercube(dim - 1), + qline, + face, + subface, + p_line); constraint_points.insert(constraint_points.end(), p_line.begin(), p_line.end()); @@ -658,8 +664,13 @@ FE_Q_Base::get_subface_interpolation_matrix( // be done for the face orientation flag in 3D. const Quadrature subface_quadrature = subface == numbers::invalid_unsigned_int ? - QProjector::project_to_face(quad_face_support, 0) : - QProjector::project_to_subface(quad_face_support, 0, subface); + QProjector::project_to_face(this->reference_cell_type(), + quad_face_support, + 0) : + QProjector::project_to_subface(this->reference_cell_type(), + quad_face_support, + 0, + subface); for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i) { const Point &p = subface_quadrature.point(i); diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index dc38e9d687..ce1b950f14 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -14,7 +14,6 @@ // --------------------------------------------------------------------- -#include #include #include #include diff --git a/source/fe/fe_q_dg0.cc b/source/fe/fe_q_dg0.cc index 3f1b209c35..e85acf4dc5 100644 --- a/source/fe/fe_q_dg0.cc +++ b/source/fe/fe_q_dg0.cc @@ -14,7 +14,6 @@ // --------------------------------------------------------------------- -#include #include #include #include diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 19e0582061..310806e5bc 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -191,7 +191,8 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) } Quadrature faces = - QProjector::project_to_all_faces(face_points); + QProjector::project_to_all_faces(this->reference_cell_type(), + face_points); for (; current < GeometryInfo::faces_per_cell * n_face_points; ++current) { @@ -199,7 +200,12 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) // into the vector this->generalized_support_points[current] = faces.point(current + QProjector::DataSetDescriptor::face( - 0, true, false, false, n_face_points)); + this->reference_cell_type(), + 0, + true, + false, + false, + n_face_points)); } } @@ -276,7 +282,10 @@ FE_RaviartThomas::initialize_restriction() // child cell are evaluated // in the quadrature points // of a full face. - Quadrature q_face = QProjector::project_to_face(q_base, face); + Quadrature q_face = + QProjector::project_to_face(this->reference_cell_type(), + q_base, + face); // Store shape values, since the // evaluation suffers if not // ordered by point @@ -294,8 +303,8 @@ FE_RaviartThomas::initialize_restriction() // the coarse face are // evaluated on the subface // only. - Quadrature q_sub = - QProjector::project_to_subface(q_base, face, sub); + Quadrature q_sub = QProjector::project_to_subface( + this->reference_cell_type(), q_base, face, sub); const unsigned int child = GeometryInfo::child_cell_on_face( RefinementCase::isotropic_refinement, face, sub); @@ -372,7 +381,10 @@ FE_RaviartThomas::initialize_restriction() for (unsigned int child = 0; child < GeometryInfo::max_children_per_cell; ++child) { - Quadrature q_sub = QProjector::project_to_child(q_cell, child); + Quadrature q_sub = + QProjector::project_to_child(this->reference_cell_type(), + q_cell, + child); for (unsigned int k = 0; k < q_sub.size(); ++k) for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell(); diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 19b150ef5d..0ddeb83502 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -166,13 +166,19 @@ FE_RaviartThomasNodal::initialize_support_points(const unsigned int deg) for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) this->generalized_face_support_points[k] = face_points.point(k); Quadrature faces = - QProjector::project_to_all_faces(face_points); + QProjector::project_to_all_faces(this->reference_cell_type(), + face_points); for (unsigned int k = 0; k < this->n_dofs_per_face() * GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = - faces.point(k + QProjector::DataSetDescriptor::face( - 0, true, false, false, this->n_dofs_per_face())); + this->generalized_support_points[k] = faces.point( + k + + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face())); current = this->n_dofs_per_face() * GeometryInfo::faces_per_cell; } @@ -621,7 +627,9 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( // matrix by simply taking the // value at the support points. const Quadrature face_projection = - QProjector::project_to_face(quad_face_support, 0); + QProjector::project_to_face(this->reference_cell_type(), + quad_face_support, + 0); for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i) { @@ -727,7 +735,10 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( // value at the support points. const Quadrature subface_projection = - QProjector::project_to_subface(quad_face_support, 0, subface); + QProjector::project_to_subface(this->reference_cell_type(), + quad_face_support, + 0, + subface); for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i) { diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc index aacdccfcfd..b7c66bdf63 100644 --- a/source/fe/fe_rt_bubbles.cc +++ b/source/fe/fe_rt_bubbles.cc @@ -147,13 +147,19 @@ FE_RT_Bubbles::initialize_support_points(const unsigned int deg) for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) this->generalized_face_support_points[k] = face_points.point(k); Quadrature faces = - QProjector::project_to_all_faces(face_points); + QProjector::project_to_all_faces(this->reference_cell_type(), + face_points); for (unsigned int k = 0; k < this->n_dofs_per_face() * GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = - faces.point(k + QProjector::DataSetDescriptor::face( - 0, true, false, false, this->n_dofs_per_face())); + this->generalized_support_points[k] = faces.point( + k + + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face())); current = this->n_dofs_per_face() * GeometryInfo::faces_per_cell; } diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 2c169eb84e..ef668c0171 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -14,7 +14,6 @@ // --------------------------------------------------------------------- #include -#include #include #include diff --git a/source/fe/fe_trace.cc b/source/fe/fe_trace.cc index 97cb68e8b8..577c3ae778 100644 --- a/source/fe/fe_trace.cc +++ b/source/fe/fe_trace.cc @@ -15,7 +15,6 @@ #include -#include #include #include diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index c6e5918588..2affce2eb5 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -116,7 +116,8 @@ MappingCartesian::get_face_data( { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique( - QProjector::project_to_all_faces(quadrature)); + QProjector::project_to_all_faces(ReferenceCell::get_hypercube(dim), + quadrature)); auto &data = dynamic_cast(*data_ptr); // verify that we have computed the transitive hull of the required @@ -140,8 +141,8 @@ MappingCartesian::get_subface_data( const Quadrature &quadrature) const { std::unique_ptr::InternalDataBase> data_ptr = - std::make_unique( - QProjector::project_to_all_subfaces(quadrature)); + std::make_unique(QProjector::project_to_all_subfaces( + ReferenceCell::get_hypercube(dim), quadrature)); auto &data = dynamic_cast(*data_ptr); // verify that we have computed the transitive hull of the required @@ -203,8 +204,7 @@ MappingCartesian::maybe_update_cell_quadrature_points( { if (data.update_each & update_quadrature_points) { - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::cell(); + const auto offset = QProjector::DataSetDescriptor::cell(); transform_quadrature_points(cell, data, offset, quadrature_points); } @@ -224,13 +224,13 @@ MappingCartesian::maybe_update_face_quadrature_points( if (data.update_each & update_quadrature_points) { - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::face(face_no, - cell->face_orientation( - face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature_points.size()); + const auto offset = QProjector::DataSetDescriptor::face( + ReferenceCell::get_hypercube(dim), + face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature_points.size()); transform_quadrature_points(cell, data, offset, quadrature_points); @@ -257,15 +257,15 @@ MappingCartesian::maybe_update_subface_quadrature_points( if (data.update_each & update_quadrature_points) { - const typename QProjector::DataSetDescriptor offset = - QProjector::DataSetDescriptor::subface( - face_no, - sub_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature_points.size(), - cell->subface_case(face_no)); + const auto offset = QProjector::DataSetDescriptor::subface( + ReferenceCell::get_hypercube(dim), + face_no, + sub_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature_points.size(), + cell->subface_case(face_no)); transform_quadrature_points(cell, data, offset, quadrature_points); } diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 62824c3f64..163719f12e 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -628,7 +628,9 @@ MappingFEField::get_face_data( std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); auto & data = dynamic_cast(*data_ptr); - const Quadrature q(QProjector::project_to_all_faces(quadrature)); + const Quadrature q( + QProjector::project_to_all_faces(ReferenceCell::get_hypercube(dim), + quadrature)); this->compute_face_data(update_flags, q, quadrature.size(), data); return data_ptr; @@ -644,7 +646,9 @@ MappingFEField::get_subface_data( std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); auto & data = dynamic_cast(*data_ptr); - const Quadrature q(QProjector::project_to_all_subfaces(quadrature)); + const Quadrature q( + QProjector::project_to_all_subfaces(ReferenceCell::get_hypercube(dim), + quadrature)); this->compute_face_data(update_flags, q, quadrature.size(), data); return data_ptr; @@ -1763,7 +1767,9 @@ MappingFEField::fill_fe_face_values( cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::face(face_no, + QProjector::DataSetDescriptor::face(ReferenceCell::get_hypercube( + dim), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), @@ -1803,14 +1809,15 @@ MappingFEField:: cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::subface(face_no, - subface_no, - cell->face_orientation( - face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size(), - cell->subface_case(face_no)), + QProjector::DataSetDescriptor::subface( + ReferenceCell::get_hypercube(dim), + face_no, + subface_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)), quadrature, data, euler_dof_handler->get_fe(), diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 2a22a445b8..b8dc584573 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -289,7 +289,8 @@ MappingManifold::get_face_data( std::make_unique(); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), - QProjector::project_to_all_faces(quadrature), + QProjector::project_to_all_faces( + ReferenceCell::get_hypercube(dim), quadrature), quadrature.size()); return data_ptr; @@ -307,7 +308,8 @@ MappingManifold::get_subface_data( std::make_unique(); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), - QProjector::project_to_all_subfaces(quadrature), + QProjector::project_to_all_subfaces( + ReferenceCell::get_hypercube(dim), quadrature), quadrature.size()); return data_ptr; @@ -1221,7 +1223,8 @@ MappingManifold::fill_fe_face_values( cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::face(face_no, + QProjector::DataSetDescriptor::face(ReferenceCell::get_hypercube(dim), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), @@ -1254,7 +1257,9 @@ MappingManifold::fill_fe_subface_values( cell, face_no, subface_no, - QProjector::DataSetDescriptor::subface(face_no, + QProjector::DataSetDescriptor::subface(ReferenceCell::get_hypercube( + dim), + face_no, subface_no, cell->face_orientation(face_no), cell->face_flip(face_no), diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index e812cb30c8..090d9a23b5 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -16,7 +16,6 @@ #include #include -#include #include #include #include diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 014403e94b..49a0eea8e4 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -2695,7 +2695,8 @@ MappingQGeneric::get_face_data( std::make_unique(polynomial_degree); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), - QProjector::project_to_all_faces(quadrature), + QProjector::project_to_all_faces( + ReferenceCell::get_hypercube(dim), quadrature), quadrature.size()); return data_ptr; @@ -2713,7 +2714,8 @@ MappingQGeneric::get_subface_data( std::make_unique(polynomial_degree); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), - QProjector::project_to_all_subfaces(quadrature), + QProjector::project_to_all_subfaces( + ReferenceCell::get_hypercube(dim), quadrature), quadrature.size()); return data_ptr; @@ -3235,7 +3237,8 @@ MappingQGeneric::fill_fe_face_values( cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::face(face_no, + QProjector::DataSetDescriptor::face(ReferenceCell::get_hypercube(dim), + face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), @@ -3281,7 +3284,9 @@ MappingQGeneric::fill_fe_subface_values( cell, face_no, subface_no, - QProjector::DataSetDescriptor::subface(face_no, + QProjector::DataSetDescriptor::subface(ReferenceCell::get_hypercube( + dim), + face_no, subface_no, cell->face_orientation(face_no), cell->face_flip(face_no), diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc index 330655c47a..2462f1ab7d 100644 --- a/source/grid/grid_out.cc +++ b/source/grid/grid_out.cc @@ -4031,7 +4031,9 @@ namespace internal std::vector dummy_weights(n_points, 1. / n_points); Quadrature quadrature(boundary_points, dummy_weights); - q_projector = QProjector::project_to_all_faces(quadrature); + q_projector = + QProjector::project_to_all_faces(ReferenceCell::Type::Quad, + quadrature); } for (const auto &cell : tria.active_cell_iterators()) -- 2.39.5