From acc4598ac89094b291934a0fed9474f504f20987 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 25 Jan 2023 21:44:21 -0500 Subject: [PATCH] Un-hardcode some orientation info in FE_SimplexP{_Bubbles}. --- source/fe/fe_simplex_p.cc | 71 ++++++------------------------- source/fe/fe_simplex_p_bubbles.cc | 62 ++++++++++++++++++--------- 2 files changed, 57 insertions(+), 76 deletions(-) diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc index 18ab94784d..e5998c5daf 100644 --- a/source/fe/fe_simplex_p.cc +++ b/source/fe/fe_simplex_p.cc @@ -91,71 +91,28 @@ namespace return unit_points; } + const auto reference_cell = ReferenceCells::get_simplex(); // Piecewise constants are a special case: use a support point at the // centroid and only the centroid if (degree == 0) { - unit_points.emplace_back( - ReferenceCells::get_simplex().template barycenter()); + unit_points.emplace_back(reference_cell.template barycenter()); return unit_points; } - if (dim == 1) - { - // We don't really have dim = 1 support for simplex elements yet, but - // its convenient for populating the face array - Assert(degree <= 2, ExcNotImplemented()); - if (degree >= 1) - { - unit_points.emplace_back(0.0); - unit_points.emplace_back(1.0); + Assert(degree <= 2, ExcNotImplemented()); + for (const auto &vertex_no : reference_cell.vertex_indices()) + unit_points.emplace_back(reference_cell.template vertex(vertex_no)); - if (degree == 2) - unit_points.emplace_back(0.5); - } - } - else if (dim == 2) - { - Assert(degree <= 2, ExcNotImplemented()); - if (degree >= 1) - { - unit_points.emplace_back(0.0, 0.0); - unit_points.emplace_back(1.0, 0.0); - unit_points.emplace_back(0.0, 1.0); - - if (degree == 2) - { - unit_points.emplace_back(0.5, 0.0); - unit_points.emplace_back(0.5, 0.5); - unit_points.emplace_back(0.0, 0.5); - } - } - } - else if (dim == 3) - { - Assert(degree <= 2, ExcNotImplemented()); - if (degree >= 1) - { - unit_points.emplace_back(0.0, 0.0, 0.0); - unit_points.emplace_back(1.0, 0.0, 0.0); - unit_points.emplace_back(0.0, 1.0, 0.0); - unit_points.emplace_back(0.0, 0.0, 1.0); - - if (degree == 2) - { - unit_points.emplace_back(0.5, 0.0, 0.0); - unit_points.emplace_back(0.5, 0.5, 0.0); - unit_points.emplace_back(0.0, 0.5, 0.0); - unit_points.emplace_back(0.0, 0.0, 0.5); - unit_points.emplace_back(0.5, 0.0, 0.5); - unit_points.emplace_back(0.0, 0.5, 0.5); - } - } - } - else - { - Assert(false, ExcNotImplemented()); - } + if (degree == 2) + for (const auto &line_no : reference_cell.line_indices()) + { + const auto v0 = reference_cell.template vertex( + reference_cell.line_to_cell_vertices(line_no, 0)); + const auto v1 = reference_cell.template vertex( + reference_cell.line_to_cell_vertices(line_no, 1)); + unit_points.emplace_back((v0 + v1) / 2.0); + } return unit_points; } diff --git a/source/fe/fe_simplex_p_bubbles.cc b/source/fe/fe_simplex_p_bubbles.cc index 260e5a3847..d52d4ed11e 100644 --- a/source/fe/fe_simplex_p_bubbles.cc +++ b/source/fe/fe_simplex_p_bubbles.cc @@ -67,8 +67,8 @@ namespace FE_P_BubblesImplementation FE_SimplexP fe_p(degree); std::vector> points = fe_p.get_unit_support_points(); - const Point centroid = - fe_p.reference_cell().template barycenter(); + const auto reference_cell = fe_p.reference_cell(); + const Point centroid = reference_cell.template barycenter(); switch (dim) { @@ -85,11 +85,27 @@ namespace FE_P_BubblesImplementation { if (degree == 2) { - const double q13 = 1.0 / 3.0; - points.emplace_back(q13, q13, 0.0); - points.emplace_back(q13, 0.0, q13); - points.emplace_back(0.0, q13, q13); - points.emplace_back(q13, q13, q13); + for (const auto &face_no : reference_cell.face_indices()) + { + Point midpoint; + for (const auto face_vertex_no : + reference_cell.face_reference_cell(0).vertex_indices()) + { + const auto vertex_no = + reference_cell.face_to_cell_vertices( + face_no, + face_vertex_no, + ReferenceCell::default_combined_face_orientation()); + + midpoint += + reference_cell.template vertex(vertex_no); + } + + midpoint /= + reference_cell.face_reference_cell(0).n_vertices(); + points.push_back(midpoint); + } + points.push_back(centroid); } return points; @@ -117,8 +133,8 @@ namespace FE_P_BubblesImplementation BarycentricPolynomials get_basis(const unsigned int degree) { - const Point centroid = - ReferenceCells::get_simplex().template barycenter(); + const auto reference_cell = ReferenceCells::get_simplex(); + const Point centroid = reference_cell.template barycenter(); auto M = [](const unsigned int d) { return BarycentricPolynomial::monomial(d); @@ -140,8 +156,8 @@ namespace FE_P_BubblesImplementation // in 2D and 3D we add a centroid bubble function auto c_bubble = BarycentricPolynomial() + 1; - for (unsigned int d = 0; d < dim + 1; ++d) - c_bubble = c_bubble * M(d); + for (const auto &vertex : reference_cell.vertex_indices()) + c_bubble = c_bubble * M(vertex); c_bubble = c_bubble / c_bubble.value(centroid); std::vector> bubble_functions; @@ -154,14 +170,22 @@ namespace FE_P_BubblesImplementation // need 'face bubble' functions in addition to the centroid. // Furthermore we need to subtract them off from the other // functions so that we end up with an interpolatory basis - auto b0 = 27 * M(0) * M(1) * M(2); - bubble_functions.push_back(b0 - b0.value(centroid) * c_bubble); - auto b1 = 27 * M(0) * M(1) * M(3); - bubble_functions.push_back(b1 - b1.value(centroid) * c_bubble); - auto b2 = 27 * M(0) * M(2) * M(3); - bubble_functions.push_back(b2 - b2.value(centroid) * c_bubble); - auto b3 = 27 * M(1) * M(2) * M(3); - bubble_functions.push_back(b3 - b3.value(centroid) * c_bubble); + for (const auto &face_no : reference_cell.face_indices()) + { + std::vector vertices; + for (const auto face_vertex_no : + reference_cell.face_reference_cell(0).vertex_indices()) + vertices.push_back(reference_cell.face_to_cell_vertices( + face_no, + face_vertex_no, + ReferenceCell::default_combined_face_orientation())); + + Assert(vertices.size() == 3, ExcInternalError()); + auto b = + 27.0 * M(vertices[0]) * M(vertices[1]) * M(vertices[2]); + bubble_functions.push_back(b - + b.value(centroid) * c_bubble); + } bubble_functions.push_back(c_bubble); } -- 2.39.5