From: Peter Munch Date: Wed, 23 Jun 2021 08:56:19 +0000 (+0200) Subject: Work on get_face_quadrature_collection() and get_unique_face_quadratures() X-Git-Tag: v9.4.0-rc1~1199^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12491%2Fhead;p=dealii.git Work on get_face_quadrature_collection() and get_unique_face_quadratures() --- diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 898c7a192b..2a9b360878 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -328,9 +328,12 @@ namespace internal { this->n_q_points_face = quad_face[0].size(); - n_q_points_faces.resize(quad_face.size()); - for (unsigned int i = 0; i < quad_face.size(); ++i) - n_q_points_faces[i] = quad_face[i].size(); + const unsigned int n_faces = temp.first.n_faces(); + + n_q_points_faces.resize(n_faces); + for (unsigned int i = 0; i < n_faces; ++i) + n_q_points_faces[i] = + quad_face[quad_face.size() == 1 ? 0 : i].size(); unsigned int n_q_points_face_max = 0; @@ -340,8 +343,7 @@ namespace internal unsigned int n_max_vertices = 0; - for (unsigned int face_no = 0; face_no < quad_face.size(); - ++face_no) + for (unsigned int face_no = 0; face_no < n_faces; ++face_no) n_max_vertices = std::max( n_max_vertices, reference_cell.face_reference_cell(face_no).n_vertices()); @@ -353,16 +355,16 @@ namespace internal const unsigned int n_max_face_orientations = dim == 2 ? 2 : (2 * n_max_vertices); - shape_values_face.reinit({quad_face.size(), + shape_values_face.reinit({n_faces, n_max_face_orientations, n_dofs * n_q_points_face_max}); - shape_gradients_face.reinit({quad_face.size(), + shape_gradients_face.reinit({n_faces, n_max_face_orientations, dim, n_dofs * n_q_points_face_max}); - for (unsigned int f = 0; f < quad_face.size(); ++f) + for (unsigned int f = 0; f < n_faces; ++f) { const unsigned int n_face_orientations = dim == 2 ? @@ -370,7 +372,8 @@ namespace internal (2 * reference_cell.face_reference_cell(f).n_vertices()); - const unsigned int n_q_points_face = quad_face[f].size(); + const unsigned int n_q_points_face = + quad_face[quad_face.size() == 1 ? 0 : f].size(); for (unsigned int o = 0; o < n_face_orientations; ++o) { diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h index 6b43808cfd..323f2093e0 100644 --- a/include/deal.II/matrix_free/util.h +++ b/include/deal.II/matrix_free/util.h @@ -34,6 +34,11 @@ namespace internal { namespace MatrixFreeFunctions { + /** + * Given a quadrature rule @p quad defined on a cell, return the type of + * the cell and a collection of lower-dimensional quadrature rules that + * are defined on each face. + */ template inline std::pair> get_face_quadrature_collection(const Quadrature &quad, @@ -43,29 +48,15 @@ namespace internal { for (unsigned int i = 1; i <= 4; ++i) if (quad == QGaussSimplex(i)) - { - QGaussSimplex tri(i); - - if (dim == 2) - return {ReferenceCells::Triangle, - dealii::hp::QCollection(tri, tri, tri)}; - else - return {ReferenceCells::Tetrahedron, - dealii::hp::QCollection(tri, tri, tri, tri)}; - } + return {ReferenceCells::get_simplex(), + dealii::hp::QCollection( + QGaussSimplex(i))}; for (unsigned int i = 1; i <= 5; ++i) if (quad == QWitherdenVincentSimplex(i)) - { - QWitherdenVincentSimplex tri(i); - - if (dim == 2) - return {ReferenceCells::Triangle, - dealii::hp::QCollection(tri, tri, tri)}; - else - return {ReferenceCells::Tetrahedron, - dealii::hp::QCollection(tri, tri, tri, tri)}; - } + return {ReferenceCells::get_simplex(), + dealii::hp::QCollection( + QWitherdenVincentSimplex(i))}; } if (dim == 3) @@ -92,6 +83,13 @@ namespace internal dealii::hp::QCollection(quad, tri, tri, tri, tri)}; } + // note: handle hypercubes last since normally this function is not + // called for hypercubes + for (unsigned int i = 1; i <= 5; ++i) + if (quad == QGauss(i)) + return {ReferenceCells::get_hypercube(), + dealii::hp::QCollection(QGauss(i))}; + if (do_assert) AssertThrow(false, ExcNotImplemented()); @@ -100,19 +98,41 @@ namespace internal + /** + * Return face quadrature rules. In contrast to + * get_face_quadrature_collection(), it does not return the quadrature + * face for each face but returns one of each type. The first entry + * of the returned pair might contain a quadrature rule defined on lines and + * quadrilaterals, while the second entry might contain a quadrature rule + * defined on a triangle. + */ template inline std::pair, Quadrature> get_unique_face_quadratures(const Quadrature &quad) { if (dim == 2 || dim == 3) - for (unsigned int i = 1; i <= 3; ++i) - if (quad == QGaussSimplex(i)) - { - if (dim == 2) - return {QGaussSimplex(i), Quadrature()}; - else - return {Quadrature(), QGaussSimplex(i)}; - } + { + for (unsigned int i = 1; i <= 4; ++i) + if (quad == QGaussSimplex(i)) + { + if (dim == 2) + return {QGaussSimplex(i), // line! + Quadrature()}; + else + return {Quadrature(), QGaussSimplex(i)}; + } + + for (unsigned int i = 1; i <= 5; ++i) + if (quad == QWitherdenVincentSimplex(i)) + { + if (dim == 2) + return {QWitherdenVincentSimplex(i), // line! + Quadrature()}; + else + return {Quadrature(), + QWitherdenVincentSimplex(i)}; + } + } if (dim == 3) for (unsigned int i = 1; i <= 3; ++i) @@ -124,9 +144,15 @@ namespace internal if (quad == QGaussPyramid(i)) return {QGauss(i), QGaussSimplex(i)}; + // note: handle hypercubes last since normally this function is not + // called for hypercubes + for (unsigned int i = 1; i <= 5; ++i) + if (quad == QGauss(i)) + return {QGauss(i), Quadrature()}; + AssertThrow(false, ExcNotImplemented()); - return {QGauss(1), QGauss(1)}; + return {Quadrature(), Quadrature()}; } } // end of namespace MatrixFreeFunctions } // end of namespace internal