From df7dd7f428cc8f147e1a81ed090eb90f22600e2a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 12 Dec 2020 15:22:33 +0100 Subject: [PATCH] Enable FEFaceValues for wedges and pyramids --- source/base/qprojector.cc | 266 +++++++++++++----- source/fe/mapping_fe.cc | 154 ++++++++++ tests/simplex/poisson_01.cc | 35 ++- ...ilinos=true.with_simplex_support=on.output | 12 +- ...ilinos=true.with_simplex_support=on.output | 12 +- 5 files changed, 379 insertions(+), 100 deletions(-) diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 5515b54652..b0cd657d78 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include @@ -642,77 +643,88 @@ QProjector<3>::project_to_all_faces( const ReferenceCell::Type reference_cell_type, const hp::QCollection<2> &quadrature) { - if (reference_cell_type == ReferenceCell::Type::Tet) - { - // 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::array, 3>, double>, 4> faces = { - {{{{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*/}}}; + const auto support_points_tri = + [](const auto &face, const auto &orientation) -> std::vector> { + // determine support point of the current line with the correct + // orientation + switch (orientation) + { + case 1: + return {{face.first[0], face.first[1], face.first[2]}}; + case 3: + return {{face.first[1], face.first[0], face.first[2]}}; + case 5: + return {{face.first[2], face.first[0], face.first[1]}}; + case 0: + return {{face.first[0], face.first[2], face.first[1]}}; + case 2: + return {{face.first[1], face.first[2], face.first[0]}}; + case 4: + return {{face.first[2], face.first[1], face.first[0]}}; + default: + Assert(false, ExcNotImplemented()); + return {{}}; + } + }; - // linear polynomial to map the reference quadrature points correctly - // on faces - const Simplex::ScalarPolynomial<2> poly(1); + const auto support_points_quad = + [](const auto &face, const auto &orientation) -> std::vector> { + switch (orientation) + { + case 1: + return {{face.first[0], face.first[1], face.first[2], face.first[3]}}; + case 3: + return {{face.first[1], face.first[3], face.first[0], face.first[2]}}; + case 5: + return {{face.first[3], face.first[2], face.first[1], face.first[0]}}; + case 7: + return {{face.first[2], face.first[0], face.first[3], face.first[1]}}; + case 0: + return {{face.first[0], face.first[2], face.first[1], face.first[3]}}; + case 2: + return {{face.first[2], face.first[3], face.first[0], face.first[1]}}; + case 4: + return {{face.first[3], face.first[1], face.first[2], face.first[0]}}; + case 6: + return {{face.first[1], face.first[0], face.first[3], face.first[2]}}; + default: + Assert(false, ExcNotImplemented()); + return {{}}; + } + }; - // new (projected) quadrature points and weights - std::vector> points; - std::vector weights; + const auto process = [&](const auto faces) { + // new (projected) quadrature points and weights + std::vector> points; + std::vector weights; + + const Simplex::ScalarPolynomial<2> poly_tri(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) + { + // linear polynomial to map the reference quadrature points correctly + // on faces + const unsigned int n_shape_functions = faces[face_no].first.size(); + + const auto &poly = + n_shape_functions == 3 ? + static_cast &>(poly_tri) : + static_cast &>(poly_quad); - // loop over all faces (triangles) ... - for (unsigned int face_no = 0; face_no < faces.size(); ++face_no) // ... and over all possible orientations - for (unsigned int orientation = 0; orientation < 6; ++orientation) + for (unsigned int orientation = 0; + orientation < (n_shape_functions * 2); + ++orientation) { const auto &face = faces[face_no]; - std::array, 3> support_points; - - // determine support point of the current line with the correct - // orientation - switch (orientation) - { - case 1: - support_points = { - {face.first[0], face.first[1], face.first[2]}}; - break; - case 3: - support_points = { - {face.first[1], face.first[0], face.first[2]}}; - break; - case 5: - support_points = { - {face.first[2], face.first[0], face.first[1]}}; - break; - case 0: - support_points = { - {face.first[0], face.first[2], face.first[1]}}; - break; - case 2: - support_points = { - {face.first[1], face.first[2], face.first[0]}}; - break; - case 4: - support_points = { - {face.first[2], face.first[1], face.first[0]}}; - break; - default: - Assert(false, ExcNotImplemented()); - } - + const auto support_points = + n_shape_functions == 3 ? support_points_tri(face, orientation) : + support_points_quad(face, orientation); // the quadrature rule to be projected ... const auto &sub_quadrature_points = @@ -726,7 +738,7 @@ QProjector<3>::project_to_all_faces( Point<3> mapped_point; // map reference quadrature point - for (unsigned int i = 0; i < 3; ++i) + for (unsigned int i = 0; i < n_shape_functions; ++i) mapped_point += support_points[i] * poly.compute_value(i, sub_quadrature_points[j]); @@ -735,18 +747,16 @@ QProjector<3>::project_to_all_faces( // scale quadrature weight const double scaling = [&]() { - const auto &supp_pts = support_points; - - const unsigned int n_shape_functions = 3; - const unsigned int dim_ = 2; - const unsigned int spacedim = 3; + const auto & supp_pts = support_points; + const unsigned int dim_ = 2; + const unsigned int spacedim = 3; double result[spacedim][dim_]; std::vector> shape_derivatives( n_shape_functions); - for (unsigned int i = 0; i < 3; ++i) + for (unsigned int i = 0; i < n_shape_functions; ++i) shape_derivatives[i] = poly.compute_1st_derivative(i, sub_quadrature_points[j]); @@ -782,9 +792,91 @@ QProjector<3>::project_to_all_faces( weights.push_back(sub_quadrature_weights[j] * scaling); } } + } - // construct new quadrature rule - return {points, weights}; + // construct new quadrature rule + return Quadrature<3>(points, weights); + }; + + if (reference_cell_type == ReferenceCell::Type::Tet) + { + // 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>> faces = { + {{{{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(faces); + } + else if (reference_cell_type == ReferenceCell::Type::Wedge) + { + const std::vector>, double>> faces = { + {{{{Point<3>(0.0, 1.0, 0.0), + Point<3>(1.0, 0.0, 0.0), + Point<3>(0.0, 0.0, 0.0)}}, + 0.5}, + {{{Point<3>(1.0, 0.0, 1.0), + Point<3>(0.0, 1.0, 1.0), + Point<3>(0.0, 0.0, 1.0)}}, + 0.5}, + {{{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}, + {{{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}}}; + + return process(faces); + } + else if (reference_cell_type == ReferenceCell::Type::Pyramid) + { + const std::vector>, double>> faces = { + {{{{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(faces); } @@ -1331,18 +1423,36 @@ QProjector::DataSetDescriptor::face( const hp::QCollection &quadrature) { if (reference_cell_type == ReferenceCell::Type::Tri || - reference_cell_type == ReferenceCell::Type::Tet) + reference_cell_type == ReferenceCell::Type::Tet || + reference_cell_type == ReferenceCell::Type::Wedge || + reference_cell_type == ReferenceCell::Type::Pyramid) { unsigned int offset = 0; + static const unsigned int X = numbers::invalid_unsigned_int; + static const std::array scale_tri = {{2, 2, 2, X, X}}; + static const std::array scale_tet = {{6, 6, 6, 6, X}}; + static const std::array scale_wedge = {{6, 6, 8, 8, 8}}; + static const std::array scale_pyramid = { + {8, 6, 6, 6, 6}}; + + const auto &scale = + (reference_cell_type == ReferenceCell::Type::Tri) ? + scale_tri : + ((reference_cell_type == ReferenceCell::Type::Tet) ? + scale_tet : + ((reference_cell_type == ReferenceCell::Type::Wedge) ? + scale_wedge : + scale_pyramid)); + if (quadrature.size() == 1) - offset = quadrature[0].size() * face_no; + offset = scale[0] * quadrature[0].size() * face_no; else for (unsigned int i = 0; i < face_no; ++i) - offset += quadrature[i].size(); + offset += scale[i] * quadrature[i].size(); if (dim == 2) - return {2 * offset + + return {offset + face_orientation * quadrature[quadrature.size() == 1 ? 0 : face_no].size()}; else if (dim == 3) @@ -1350,7 +1460,7 @@ QProjector::DataSetDescriptor::face( const unsigned int orientation = (face_flip * 2 + face_rotation) * 2 + face_orientation; - return {6 * offset + + return {offset + orientation * quadrature[quadrature.size() == 1 ? 0 : face_no].size()}; } diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index a5c428ea4c..169f944bf5 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -246,6 +246,160 @@ MappingFE::InternalData::initialize_face( for (unsigned int i = 0; i < n_original_q_points; i++) unit_tangentials[7].emplace_back(t1); } + else if (this->fe.reference_cell_type() == ReferenceCell::Type::Wedge) + { + Tensor<1, dim> t1; + constexpr int d0 = 0; + constexpr int d1 = 1 % dim; + constexpr int d2 = 2 % dim; + + // face 0 + t1[d0] = 0; + t1[d1] = 1; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[0].emplace_back(t1); + + // face 0 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[5].emplace_back(t1); + + // face 1 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[1].emplace_back(t1); + + // face 1 + t1[d0] = 0; + t1[d1] = 1; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[6].emplace_back(t1); + + // face 2 + t1[d0] = -1 / std::sqrt(2.0); + t1[d1] = +1 / std::sqrt(2.0); + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[2].emplace_back(t1); + + // face 2 + t1[d0] = 0; + t1[d1] = 0; + t1[d2] = 1; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[7].emplace_back(t1); + + // face 3 + t1[d0] = +0; + t1[d1] = +0; + t1[d2] = +1; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[3].emplace_back(t1); + + // face 3 + t1[d0] = +0; + t1[d1] = +1; + t1[d2] = +0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[8].emplace_back(t1); + + // face 4 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[4].emplace_back(t1); + + // face 4 + t1[d0] = 0; + t1[d1] = 0; + t1[d2] = 1; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[9].emplace_back(t1); + } + else if (this->fe.reference_cell_type() == ReferenceCell::Type::Pyramid) + { + Tensor<1, dim> t1; + constexpr int d0 = 0; + constexpr int d1 = 1 % dim; + constexpr int d2 = 2 % dim; + + // face 0 + t1[d0] = 0; + t1[d1] = 1; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[0].emplace_back(t1); + + // face 0 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[5].emplace_back(t1); + + // face 1 + t1[d0] = 1.0 / sqrt(2.0); + t1[d1] = 0; + t1[d2] = 1.0 / sqrt(2.0); + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[1].emplace_back(t1); + + // face 1 + t1[d0] = 0; + t1[d1] = 1; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[6].emplace_back(t1); + + // face 2 + t1[d0] = 1.0 / sqrt(2.0); + t1[d1] = 0; + t1[d2] = -1.0 / sqrt(2.0); + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[2].emplace_back(t1); + + // face 2 + t1[d0] = 0; + t1[d1] = 1; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[7].emplace_back(t1); + + // face 3 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[3].emplace_back(t1); + + // face 3 + t1[d0] = 0; + t1[d1] = 1.0 / sqrt(2.0); + t1[d2] = 1.0 / sqrt(2.0); + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[8].emplace_back(t1); + + // face 4 + t1[d0] = 1; + t1[d1] = 0; + t1[d2] = 0; + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[4].emplace_back(t1); + + // face 4 + t1[d0] = 0; + t1[d1] = +1.0 / sqrt(2.0); + t1[d2] = -1.0 / sqrt(2.0); + for (unsigned int i = 0; i < n_original_q_points; i++) + unit_tangentials[9].emplace_back(t1); + } else { Assert(false, ExcNotImplemented()); diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc index 44e9f05800..56a2020e41 100644 --- a/tests/simplex/poisson_01.cc +++ b/tests/simplex/poisson_01.cc @@ -38,6 +38,8 @@ #include #include +#include + #include #include #include @@ -96,7 +98,7 @@ void test(const Triangulation &tria, const FiniteElement &fe, const Quadrature & quad, - const Quadrature & face_quad, + const hp::QCollection & face_quad, const Mapping & mapping, const double r_boundary, const bool do_use_fe_face_values = true) @@ -230,7 +232,7 @@ test(const Triangulation &tria, if (face->at_boundary() && (face->boundary_id() == 1)) { fe_face_values->reinit(cell, face); - for (unsigned int q = 0; q < face_quad.size(); ++q) + for (const auto q : fe_face_values->quadrature_point_indices()) for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_rhs(i) += (1.0 * // 1.0 @@ -289,7 +291,9 @@ test(const Triangulation &tria, data_out.build_patches(mapping); - std::ofstream output("result.vtk"); + static unsigned int counter = 0; + + std::ofstream output("result" + std::to_string(counter++) + ".vtk"); data_out.write_vtk(output); } } @@ -372,7 +376,8 @@ test_tet(const MPI_Comm &comm, const Parameters ¶ms) Simplex::QGauss quad(params.degree + 1); - Simplex::QGauss face_quad(params.degree + 1); + hp::QCollection face_quad{ + Simplex::QGauss(params.degree + 1)}; Simplex::FE_P fe_mapping(1); MappingFE mapping(fe_mapping); @@ -415,12 +420,12 @@ test_hex(const MPI_Comm &comm, const Parameters ¶ms) QGauss quad(params.degree + 1); - QGauss quad_face(params.degree + 1); + hp::QCollection face_quad{QGauss(params.degree + 1)}; MappingQ mapping(1); // 4) Perform test (independent of mesh type) - test(tria, fe, quad, quad_face, mapping, params.p2[0]); + test(tria, fe, quad, face_quad, mapping, params.p2[0]); } template @@ -501,13 +506,18 @@ test_wedge(const MPI_Comm &comm, const Parameters ¶ms) Simplex::QGaussWedge quad(params.degree + 1); - Quadrature face_quad; // not needed + hp::QCollection face_quad{ + Simplex::QGauss(params.degree + 1), + Simplex::QGauss(params.degree + 1), + QGauss(params.degree + 1), + QGauss(params.degree + 1), + QGauss(params.degree + 1)}; Simplex::FE_WedgeP fe_mapping(1); MappingFE mapping(fe_mapping); // 4) Perform test (independent of mesh type) - test(*tria, fe, quad, face_quad, mapping, params.p2[0], false); + test(*tria, fe, quad, face_quad, mapping, params.p2[0], true); } template @@ -588,13 +598,18 @@ test_pyramid(const MPI_Comm &comm, const Parameters ¶ms) Simplex::QGaussPyramid quad(params.degree + 1); - Quadrature face_quad; // not needed + hp::QCollection face_quad{ + QGauss(params.degree + 1), + Simplex::QGauss(params.degree + 1), + Simplex::QGauss(params.degree + 1), + Simplex::QGauss(params.degree + 1), + Simplex::QGauss(params.degree + 1)}; Simplex::FE_PyramidP fe_mapping(1); MappingFE mapping(fe_mapping); // 4) Perform test (independent of mesh type) - test(*tria, fe, quad, face_quad, mapping, params.p2[0], false); + test(*tria, fe, quad, face_quad, mapping, params.p2[0], true); } int diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output index f239fef1ac..1df99dcc92 100644 --- a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output +++ b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output @@ -25,13 +25,13 @@ DEAL:: with 134 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on WEDGE mesh: DEAL:: on parallel::fullydistributed::Triangulation -DEAL:cg::Starting value 0.0132550 -DEAL:cg::Convergence step 186 value 8.72381e-13 -DEAL:: with 186 CG iterations needed to obtain convergence +DEAL:cg::Starting value 0.0561953 +DEAL:cg::Convergence step 197 value 8.61398e-13 +DEAL:: with 197 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on PYRAMID mesh: DEAL:: on parallel::fullydistributed::Triangulation -DEAL:cg::Starting value 0.0230669 -DEAL:cg::Convergence step 79 value 7.70260e-13 -DEAL:: with 79 CG iterations needed to obtain convergence +DEAL:cg::Starting value 0.101163 +DEAL:cg::Convergence step 83 value 8.81456e-13 +DEAL:: with 83 CG iterations needed to obtain convergence DEAL:: diff --git a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output index 01f8dc29c4..6c0ab9a4a1 100644 --- a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output +++ b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output @@ -25,13 +25,13 @@ DEAL:: with 134 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on WEDGE mesh: DEAL:: on parallel::fullydistributed::Triangulation -DEAL:cg::Starting value 0.0132550 -DEAL:cg::Convergence step 186 value 8.72414e-13 -DEAL:: with 186 CG iterations needed to obtain convergence +DEAL:cg::Starting value 0.0561953 +DEAL:cg::Convergence step 197 value 8.60684e-13 +DEAL:: with 197 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on PYRAMID mesh: DEAL:: on parallel::fullydistributed::Triangulation -DEAL:cg::Starting value 0.0230669 -DEAL:cg::Convergence step 79 value 7.71392e-13 -DEAL:: with 79 CG iterations needed to obtain convergence +DEAL:cg::Starting value 0.101163 +DEAL:cg::Convergence step 83 value 9.01991e-13 +DEAL:: with 83 CG iterations needed to obtain convergence DEAL:: -- 2.39.5