From: Peter Munch Date: Tue, 22 Dec 2020 21:05:56 +0000 (+0100) Subject: Generalize ShapeInfo::reinit X-Git-Tag: v9.3.0-rc1~467^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11372%2Fhead;p=dealii.git Generalize ShapeInfo::reinit --- diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 23d07de2bb..a69022a7aa 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -420,6 +420,12 @@ namespace internal */ unsigned int n_q_points_face; + /** + * Stores the number of quadrature points of a face in @p dim dimensions + * for simplex, wedge and pyramid reference cells. + */ + std::vector n_q_points_faces; + /** * Stores the number of DoFs per face in @p dim dimensions. */ diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 1ec957421e..4e77bb7e33 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -314,64 +314,97 @@ namespace internal q] = grad[d]; } - try - { - const dealii::ReferenceCell reference_cell = - dealii::ReferenceCells::get_simplex(); + { + const auto reference_cell_type = fe.reference_cell(); + + const auto temp = get_face_quadrature_collection(quad, false); + const auto &quad_face = temp.second; + + if (reference_cell_type != temp.first) + { + // TODO: this might happen if the quadrature rule and the + // the FE do not match + this->n_q_points_face = 0; + } + else + { + 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 auto quad_face = get_face_quadrature(quad); - this->n_q_points_face = quad_face.size(); + unsigned int n_q_points_face_max = 0; - const unsigned int n_face_orientations = dim == 2 ? 2 : 6; - const unsigned int n_faces = - dealii::internal::ReferenceCell::get_cell(reference_cell) - .n_faces(); + for (unsigned int i = 0; i < quad_face.size(); ++i) + n_q_points_face_max = + std::max(n_q_points_face_max, quad_face[i].size()); - const auto projected_quad_face = - QProjector::project_to_all_faces(reference_cell, - quad_face); + unsigned int n_max_vertices = 0; - shape_values_face.reinit( - {n_faces, n_face_orientations, n_dofs * n_q_points_face}); + for (unsigned int face_no = 0; face_no < quad_face.size(); + ++face_no) + n_max_vertices = std::max(n_max_vertices, + internal::ReferenceCell::get_face( + reference_cell_type, face_no) + .n_vertices()); - shape_gradients_face.reinit( - {n_faces, n_face_orientations, dim, n_dofs * n_q_points_face}); + const auto projected_quad_face = + QProjector::project_to_all_faces(reference_cell_type, + quad_face); - for (unsigned int f = 0; f < n_faces; ++f) - for (unsigned int o = 0; o < n_face_orientations; ++o) + const unsigned int n_max_face_orientations = + dim == 2 ? 2 : (2 * n_max_vertices); + + shape_values_face.reinit({quad_face.size(), + n_max_face_orientations, + n_dofs * n_q_points_face_max}); + + shape_gradients_face.reinit({quad_face.size(), + n_max_face_orientations, + dim, + n_dofs * n_q_points_face_max}); + + for (unsigned int f = 0; f < quad_face.size(); ++f) { - const auto offset = - QProjector::DataSetDescriptor::face( - reference_cell, - f, - (o ^ 1) & 1, // face_orientation - (o >> 1) & 1, // face_flip - (o >> 2) & 1, // face_rotation - n_q_points_face); - - for (unsigned int i = 0; i < n_dofs; ++i) - for (unsigned int q = 0; q < n_q_points_face; ++q) - { - const auto point = - projected_quad_face.point(q + offset); - - shape_values_face(f, o, i * n_q_points_face + q) = - fe.shape_value(i, point); - - const auto grad = fe.shape_grad(i, point); - - for (int d = 0; d < dim; ++d) - shape_gradients_face( - f, o, d, i * n_q_points_face + q) = grad[d]; - } + const unsigned int n_face_orientations = + dim == 2 ? 2 : + (2 * internal::ReferenceCell::get_face( + reference_cell_type, f) + .n_vertices()); + + const unsigned int n_q_points_face = quad_face[f].size(); + + for (unsigned int o = 0; o < n_face_orientations; ++o) + { + const auto offset = + QProjector::DataSetDescriptor::face( + reference_cell_type, + f, + (o ^ 1) & 1, // face_orientation + (o >> 1) & 1, // face_flip + (o >> 2) & 1, // face_rotation + quad_face); + + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const auto point = + projected_quad_face.point(q + offset); + + shape_values_face(f, o, i * n_q_points_face + q) = + fe.shape_value(i, point); + + const auto grad = fe.shape_grad(i, point); + + for (int d = 0; d < dim; ++d) + shape_gradients_face( + f, o, d, i * n_q_points_face + q) = grad[d]; + } + } } - } - catch (...) - { - // TODO: this might happen if the quadrature rule and the - // the FE do not match - this->n_q_points_face = 0; - } + } + } // TODO: also fill shape_hessians, inverse_shape_values, // shape_data_on_face, quadrature_data_on_face, diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h index b46f3c27f5..495c6a6e6d 100644 --- a/include/deal.II/matrix_free/util.h +++ b/include/deal.II/matrix_free/util.h @@ -22,6 +22,10 @@ #include +#include + +#include + #include DEAL_II_NAMESPACE_OPEN @@ -45,6 +49,70 @@ namespace internal return Quadrature(); } + template + inline std::pair> + get_face_quadrature_collection(const Quadrature &quad, + const bool do_assert = true) + { + if (dim == 2 || dim == 3) + { + for (unsigned int i = 1; i <= 4; ++i) + if (quad == Simplex::QGauss(i)) + { + Simplex::QGauss 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)}; + } + + for (unsigned int i = 1; i <= 5; ++i) + if (quad == Simplex::QWitherdenVincent(i)) + { + Simplex::QWitherdenVincent 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)}; + } + } + + if (dim == 3) + for (unsigned int i = 1; i <= 3; ++i) + if (quad == Simplex::QGaussWedge(i)) + { + QGauss quad(i); + Simplex::QGauss tri(i); + + return { + ReferenceCells::Wedge, + dealii::hp::QCollection(tri, tri, quad, quad, quad)}; + } + + if (dim == 3) + for (unsigned int i = 1; i <= 2; ++i) + if (quad == Simplex::QGaussPyramid(i)) + { + QGauss quad(i); + Simplex::QGauss tri(i); + + return { + ReferenceCells::Pyramid, + dealii::hp::QCollection(quad, tri, tri, tri, tri)}; + } + + if (do_assert) + AssertThrow(false, ExcNotImplemented()); + + return {ReferenceCells::Invalid, dealii::hp::QCollection()}; + } + } // end of namespace MatrixFreeFunctions } // end of namespace internal