From 2deff43ada77021f1068da3139059437a0b3aa7f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 8 Feb 2023 18:05:09 -0700 Subject: [PATCH] Do not pass information that is not necessary. We are passing areas of faces to lambda functions, but they do not actually use it. So avoid doing that. --- source/base/qprojector.cc | 310 ++++++++++++++++++-------------------- 1 file changed, 145 insertions(+), 165 deletions(-) diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index fac848e931..d42eedf21a 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -644,213 +644,193 @@ QProjector<3>::project_to_all_faces(const ReferenceCell & reference_cell, const hp::QCollection<2> &quadrature) { const auto support_points_tri = - [](const std::pair>, double> &face, - const unsigned char orientation) -> std::vector> { + [](const std::vector> &face, + const unsigned char orientation) -> std::vector> { const boost::container::small_vector, 8> temp = ReferenceCells::Triangle.permute_by_combined_orientation>( - face.first, orientation); + face, orientation); return std::vector>(temp.begin(), temp.end()); }; const auto support_points_quad = - [](const std::pair>, double> &face, - const unsigned char orientation) -> std::vector> { + [](const std::vector> &face, + const unsigned char orientation) -> std::vector> { const boost::container::small_vector, 8> temp = ReferenceCells::Quadrilateral.permute_by_combined_orientation>( - face.first, orientation); + face, orientation); return std::vector>(temp.begin(), temp.end()); }; - const auto process = - [&](const std::vector>, double>> &faces) { - // new (projected) quadrature points and weights - std::vector> points; - std::vector weights; + const auto process = [&](const std::vector>> &faces) { + // new (projected) quadrature points and weights + std::vector> points; + std::vector weights; - const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1); - const TensorProductPolynomials<2> poly_quad( - Polynomials::generate_complete_Lagrange_basis( - {Point<1>(0.0), Point<1>(1.0)})); + const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1); + const TensorProductPolynomials<2> poly_quad( + Polynomials::generate_complete_Lagrange_basis( + {Point<1>(0.0), Point<1>(1.0)})); - // loop over all faces (triangles) ... - for (unsigned int face_no = 0; face_no < faces.size(); ++face_no) - { - // We will use linear polynomials to map the reference quadrature - // points correctly to on faces. There are as many linear shape - // functions as there are vertices in the face. - const unsigned int n_linear_shape_functions = - faces[face_no].first.size(); - - const auto &poly = - (n_linear_shape_functions == 3 ? - static_cast &>(poly_tri) : - static_cast &>(poly_quad)); - - // ... and over all possible orientations - for (unsigned char orientation = 0; - orientation < reference_cell.n_face_orientations(face_no); - ++orientation) - { - const auto &face = faces[face_no]; + // loop over all faces (triangles) ... + for (unsigned int face_no = 0; face_no < faces.size(); ++face_no) + { + // We will use linear polynomials to map the reference quadrature + // points correctly to on faces. There are as many linear shape + // functions as there are vertices in the face. + const unsigned int n_linear_shape_functions = faces[face_no].size(); - const auto support_points = - n_linear_shape_functions == 3 ? - support_points_tri(face, orientation) : - support_points_quad(face, orientation); + const auto &poly = + (n_linear_shape_functions == 3 ? + static_cast &>(poly_tri) : + static_cast &>(poly_quad)); - // the quadrature rule to be projected ... - const auto &sub_quadrature_points = - quadrature[quadrature.size() == 1 ? 0 : face_no].get_points(); - const auto &sub_quadrature_weights = - quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights(); + // ... and over all possible orientations + for (unsigned char orientation = 0; + orientation < reference_cell.n_face_orientations(face_no); + ++orientation) + { + const auto &face = faces[face_no]; - // loop over all quadrature points - for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j) - { - Point<3> mapped_point; + const auto support_points = + n_linear_shape_functions == 3 ? + support_points_tri(face, orientation) : + support_points_quad(face, orientation); - // map reference quadrature point - for (unsigned int i = 0; i < n_linear_shape_functions; ++i) - mapped_point += - support_points[i] * - poly.compute_value(i, sub_quadrature_points[j]); + // the quadrature rule to be projected ... + const auto &sub_quadrature_points = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_points(); + const auto &sub_quadrature_weights = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights(); - points.push_back(mapped_point); + // loop over all quadrature points + for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j) + { + Point<3> mapped_point; - // scale quadrature weight - const double scaling = [&]() { - const auto & supp_pts = support_points; - const unsigned int dim_ = 2; - const unsigned int spacedim = 3; + // map reference quadrature point + for (unsigned int i = 0; i < n_linear_shape_functions; ++i) + mapped_point += + support_points[i] * + poly.compute_value(i, sub_quadrature_points[j]); - double result[spacedim][dim_]; + points.push_back(mapped_point); - std::vector> shape_derivatives( - n_linear_shape_functions); + // scale quadrature weight + const double scaling = [&]() { + const auto & supp_pts = support_points; + const unsigned int dim_ = 2; + const unsigned int spacedim = 3; - for (unsigned int i = 0; i < n_linear_shape_functions; ++i) - shape_derivatives[i] = - poly.compute_1st_derivative(i, - sub_quadrature_points[j]); + double result[spacedim][dim_]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim_; ++j) - result[i][j] = shape_derivatives[0][j] * supp_pts[0][i]; - for (unsigned int k = 1; k < n_linear_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim_; ++j) - result[i][j] += - shape_derivatives[k][j] * supp_pts[k][i]; + std::vector> shape_derivatives( + n_linear_shape_functions); - DerivativeForm<1, dim_, spacedim> contravariant; + for (unsigned int i = 0; i < n_linear_shape_functions; ++i) + shape_derivatives[i] = + poly.compute_1st_derivative(i, sub_quadrature_points[j]); + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim_; ++j) + result[i][j] = shape_derivatives[0][j] * supp_pts[0][i]; + for (unsigned int k = 1; k < n_linear_shape_functions; ++k) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < dim_; ++j) - contravariant[i][j] = result[i][j]; + result[i][j] += + shape_derivatives[k][j] * supp_pts[k][i]; + DerivativeForm<1, dim_, spacedim> contravariant; - Tensor<1, spacedim> DX_t[dim_]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim_; ++j) - DX_t[j][i] = contravariant[i][j]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim_; ++j) + contravariant[i][j] = result[i][j]; - Tensor<2, dim_> G; - for (unsigned int i = 0; i < dim_; ++i) - for (unsigned int j = 0; j < dim_; ++j) - G[i][j] = DX_t[i] * DX_t[j]; - return std::sqrt(determinant(G)); - }(); + Tensor<1, spacedim> DX_t[dim_]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim_; ++j) + DX_t[j][i] = contravariant[i][j]; - weights.push_back(sub_quadrature_weights[j] * scaling); - } - } - } + Tensor<2, dim_> G; + for (unsigned int i = 0; i < dim_; ++i) + for (unsigned int j = 0; j < dim_; ++j) + G[i][j] = DX_t[i] * DX_t[j]; - // construct new quadrature rule - return Quadrature<3>(points, weights); - }; + return std::sqrt(determinant(G)); + }(); + + weights.push_back(sub_quadrature_weights[j] * scaling); + } + } + } + + // construct new quadrature rule + return Quadrature<3>(points, weights); + }; if (reference_cell == ReferenceCells::Tetrahedron) { - // reference faces (defined by its support points and its area) - // note: the area is later not used as a scaling factor but recomputed - const std::vector>, double>> - face_vertex_locations_and_area = { - {{{{Point<3>(0.0, 0.0, 0.0), - Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 1.0, 0.0)}}, - 0.5}, - {{{Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 0.0, 0.0), - Point<3>(0.0, 0.0, 1.0)}}, - 0.5}, - {{{Point<3>(0.0, 0.0, 0.0), - Point<3>(0.0, 1.0, 0.0), - Point<3>(0.0, 0.0, 1.0)}}, - 0.5}, - {{{Point<3>(0.0, 1.0, 0.0), - Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 0.0, 1.0)}}, - 0.5 * sqrt(3.0) /*equilateral triangle*/}}}; - - return process(face_vertex_locations_and_area); + const std::vector>> face_vertex_locations = { + {{Point<3>(0.0, 0.0, 0.0), + Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 1.0, 0.0)}}, + {{Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 1.0)}}, + {{Point<3>(0.0, 0.0, 0.0), + Point<3>(0.0, 1.0, 0.0), + Point<3>(0.0, 0.0, 1.0)}}, + {{Point<3>(0.0, 1.0, 0.0), + Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 1.0)}}}; + + return process(face_vertex_locations); } else if (reference_cell == ReferenceCells::Wedge) { - const std::vector>, double>> - face_vertex_locations_and_area = {{{{{Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 0.0, 0.0), - Point<3>(0.0, 1.0, 0.0)}}, - 0.5}, - {{{Point<3>(0.0, 0.0, 1.0), - Point<3>(1.0, 0.0, 1.0), - Point<3>(0.0, 1.0, 1.0)}}, - 0.5}, - {{{Point<3>(0.0, 0.0, 0.0), - Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 0.0, 1.0), - Point<3>(1.0, 0.0, 1.0)}}, - 1.0}, - {{{Point<3>(1.0, 0.0, 0.0), - Point<3>(0.0, 1.0, 0.0), - Point<3>(1.0, 0.0, 1.0), - Point<3>(0.0, 1.0, 1.0)}}, - std::sqrt(2.0)}, - {{{Point<3>(0.0, 1.0, 0.0), - Point<3>(0.0, 0.0, 0.0), - Point<3>(0.0, 1.0, 1.0), - Point<3>(0.0, 0.0, 1.0)}}, - 1.0}}}; - - return process(face_vertex_locations_and_area); + const std::vector>> face_vertex_locations = { + {{Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 0.0), + Point<3>(0.0, 1.0, 0.0)}}, + {{Point<3>(0.0, 0.0, 1.0), + Point<3>(1.0, 0.0, 1.0), + Point<3>(0.0, 1.0, 1.0)}}, + {{Point<3>(0.0, 0.0, 0.0), + Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 1.0), + Point<3>(1.0, 0.0, 1.0)}}, + {{Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 1.0, 0.0), + Point<3>(1.0, 0.0, 1.0), + Point<3>(0.0, 1.0, 1.0)}}, + {{Point<3>(0.0, 1.0, 0.0), + Point<3>(0.0, 0.0, 0.0), + Point<3>(0.0, 1.0, 1.0), + Point<3>(0.0, 0.0, 1.0)}}}; + + return process(face_vertex_locations); } else if (reference_cell == ReferenceCells::Pyramid) { - const std::vector>, double>> - face_vertex_locations_and_area = {{{{{Point<3>(-1.0, -1.0, 0.0), - Point<3>(+1.0, -1.0, 0.0), - Point<3>(-1.0, +1.0, 0.0), - Point<3>(+1.0, +1.0, 0.0)}}, - 4.0}, - {{{Point<3>(-1.0, -1.0, 0.0), - Point<3>(-1.0, +1.0, 0.0), - Point<3>(+0.0, +0.0, 1.0)}}, - std::sqrt(2.0)}, - {{{Point<3>(+1.0, +1.0, 0.0), - Point<3>(+1.0, -1.0, 0.0), - Point<3>(+0.0, +0.0, 1.0)}}, - std::sqrt(2.0)}, - {{{Point<3>(+1.0, -1.0, 0.0), - Point<3>(-1.0, -1.0, 0.0), - Point<3>(+0.0, +0.0, 1.0)}}, - std::sqrt(2.0)}, - {{{Point<3>(-1.0, +1.0, 0.0), - Point<3>(+1.0, +1.0, 0.0), - Point<3>(+0.0, +0.0, 1.0)}}, - std::sqrt(2.0)}}}; - - return process(face_vertex_locations_and_area); + const std::vector>> face_vertex_locations = { + {{Point<3>(-1.0, -1.0, 0.0), + Point<3>(+1.0, -1.0, 0.0), + Point<3>(-1.0, +1.0, 0.0), + Point<3>(+1.0, +1.0, 0.0)}}, + {{Point<3>(-1.0, -1.0, 0.0), + Point<3>(-1.0, +1.0, 0.0), + Point<3>(+0.0, +0.0, 1.0)}}, + {{Point<3>(+1.0, +1.0, 0.0), + Point<3>(+1.0, -1.0, 0.0), + Point<3>(+0.0, +0.0, 1.0)}}, + {{Point<3>(+1.0, -1.0, 0.0), + Point<3>(-1.0, -1.0, 0.0), + Point<3>(+0.0, +0.0, 1.0)}}, + {{Point<3>(-1.0, +1.0, 0.0), + Point<3>(+1.0, +1.0, 0.0), + Point<3>(+0.0, +0.0, 1.0)}}}; + + return process(face_vertex_locations); } -- 2.39.5