From 632bb5ed4f53a80ba299605dd0febba80c3092ff Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 20 Mar 2026 18:28:03 -0700 Subject: [PATCH 1/6] fix geometric factor of periodic mesh --- src/smith/numerics/functional/domain.cpp | 74 ++++++++++++ src/smith/numerics/functional/domain.hpp | 10 ++ .../numerics/functional/geometric_factors.cpp | 27 ++++- .../numerics/functional/tests/CMakeLists.txt | 1 + .../functional/tests/dg_ghost_face_index.cpp | 2 +- .../functional/tests/dg_periodic_bc_index.cpp | 109 ++++++++++++++++++ src/smith/numerics/functional/typedefs.hpp | 4 - 7 files changed, 220 insertions(+), 7 deletions(-) create mode 100644 src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp diff --git a/src/smith/numerics/functional/domain.cpp b/src/smith/numerics/functional/domain.cpp index 9acbe708c4..148ac34c5d 100644 --- a/src/smith/numerics/functional/domain.cpp +++ b/src/smith/numerics/functional/domain.cpp @@ -375,6 +375,80 @@ Domain Domain::ofBoundaryElements(const mesh_t& mesh, std::function(mesh, func); } +/// @brief constructs a domain from some subset of the interior face elements in a mesh +template +Domain domain_of_interior_faces(const mesh_t& mesh, std::function>, int)> predicate) +{ + assert(mesh.SpaceDimension() == d); + Domain output{mesh, d - 1, Domain::Type::InteriorFaces}; + + mfem::Array face_id_to_bdr_id = mesh.GetFaceToBdrElMap(); + mfem::Vector vertices; + mesh.GetVertices(vertices); + + int edge_id = 0; + int tri_id = 0; + int quad_id = 0; + + for (int f = 0; f < mesh.GetNumFaces(); f++) { + // discard faces with the wrong type + if (!mesh.GetFaceInformation(f).IsInterior()) continue; + + auto geom = mesh.GetFaceGeometry(f); + + mfem::Array vertex_ids; + mesh.GetFaceVertices(f, vertex_ids); + + auto x = gather(vertices, vertex_ids); + + int bdr_id = face_id_to_bdr_id[f]; + int attr = (bdr_id >= 0) ? mesh.GetBdrAttribute(bdr_id) : -1; + + bool add = predicate(x, attr); + + switch (geom) { + case mfem::Geometry::SEGMENT: + if (add) { + output.edge_ids_.push_back(edge_id); + output.mfem_edge_ids_.push_back(f); + } + edge_id++; + break; + case mfem::Geometry::TRIANGLE: + if (add) { + output.tri_ids_.push_back(tri_id); + output.mfem_tri_ids_.push_back(f); + } + tri_id++; + break; + case mfem::Geometry::SQUARE: + if (add) { + output.quad_ids_.push_back(quad_id); + output.mfem_quad_ids_.push_back(f); + } + quad_id++; + break; + default: + SLIC_ERROR("unsupported element type"); + break; + } + } + + output.insert_shared_interior_face_list(); + + return output; +} + +Domain Domain::ofInteriorFaces(const mesh_t& mesh, std::function, int)> func) +{ + return domain_of_interior_faces<2>(mesh, func); +} + +Domain Domain::ofInteriorFaces(const mesh_t& mesh, std::function, int)> func) +{ + return domain_of_interior_faces<3>(mesh, func); +} + /** * @brief Get local dofs that are part of a domain, but are owned by a neighboring MPI rank * diff --git a/src/smith/numerics/functional/domain.hpp b/src/smith/numerics/functional/domain.hpp index a7fe5747e0..b96cbe2482 100644 --- a/src/smith/numerics/functional/domain.hpp +++ b/src/smith/numerics/functional/domain.hpp @@ -167,6 +167,16 @@ struct Domain { /// @overload static Domain ofBoundaryElements(const mesh_t& mesh, std::function, int)> func); + /** + * @brief create a domain from some subset of the interior faces (spatial dim == geometry dim + 1) in an mfem::Mesh + * @param mesh the entire mesh + * @param func predicate function for determining which interior faces will be included in this domain + */ + static Domain ofInteriorFaces(const mesh_t& mesh, std::function, int)> func); + + /// @overload + static Domain ofInteriorFaces(const mesh_t& mesh, std::function, int)> func); + /// @brief get elements by geometry type const std::vector& get(mfem::Geometry::Type geom) const { diff --git a/src/smith/numerics/functional/geometric_factors.cpp b/src/smith/numerics/functional/geometric_factors.cpp index ae8e5cf6d9..e97fd59365 100644 --- a/src/smith/numerics/functional/geometric_factors.cpp +++ b/src/smith/numerics/functional/geometric_factors.cpp @@ -73,8 +73,6 @@ GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry:: { const mfem::ParGridFunction* nodes = static_cast(domain.mesh_.GetNodes()); mfem::ParFiniteElementSpace* fes = nodes->ParFESpace(); - // const mfem::GridFunction* nodes = domain.mesh_.GetNodes(); - // const mfem::FiniteElementSpace* fes = nodes->FESpace(); const std::vector& element_ids = domain.get_mfem_ids(geom); @@ -82,6 +80,31 @@ GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry:: mfem::Vector X_e(int(restriction.ESize())); restriction.Gather(*nodes, X_e); + // For periodic meshes, the mesh nodes are in DG space, which will double the nodes_per_elem for interior faces. + // Therefore, we must discard half of the entries to recover H1 coordinates and correctly compute geometric factor. + // Note that faces on periodic boundaries in mfem mesh are considered interior faces with boundary attributes. + if (domain.type_ == Domain::Type::InteriorFaces && fes->IsDGSpace()) { + const uint64_t new_nodes_per_elem = restriction.nodes_per_elem / 2; + mfem::Vector X_h1(X_e.Size() / 2); + + int H1_id = 0; + for (uint64_t i = 0; i < restriction.num_elements; i++) { + for (uint64_t c = 0; c < restriction.components; c++) { + for (uint64_t j = 0; j < restriction.nodes_per_elem; j++) { + // ignore the coordinate dofs on the other side of the face + // coordinate on both side of the face should have the same values + if (j >= new_nodes_per_elem) { continue; } + + uint64_t E_id = (i * restriction.components + c) * restriction.nodes_per_elem + j; + X_h1[H1_id++] = X_e[int(E_id)]; + } + } + } + + X_e.SetSize(X_h1.Size()); + X_e = X_h1; + } + // assumes all elements are the same order int p = fes->GetElementOrder(0); diff --git a/src/smith/numerics/functional/tests/CMakeLists.txt b/src/smith/numerics/functional/tests/CMakeLists.txt index 9bc329d78d..8a09f75d04 100644 --- a/src/smith/numerics/functional/tests/CMakeLists.txt +++ b/src/smith/numerics/functional/tests/CMakeLists.txt @@ -46,6 +46,7 @@ set(functional_parallel_test_sources functional_comparisons.cpp functional_comparison_L2.cpp dg_ghost_face_index.cpp + dg_periodic_bc_index.cpp ) smith_add_tests(SOURCES ${functional_parallel_test_sources} diff --git a/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp b/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp index d65878dddb..1c3622bccd 100644 --- a/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp @@ -75,7 +75,7 @@ void L2_index_test(std::string meshfile) // note: the orientation convention is such that the normal // computed as above will point from from side 1->2 auto [u_1, u_2] = velocity; - SLIC_INFO(axom::fmt::format("One size = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), + SLIC_INFO(axom::fmt::format("One side = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), axom::fmt::streamed(u_2), axom::fmt::streamed(u_1 - u_2))); auto a = dot(u_2 - u_1, n); diff --git a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp new file mode 100644 index 0000000000..acb7e3293d --- /dev/null +++ b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp @@ -0,0 +1,109 @@ + +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include + +#include "mfem.hpp" + +#include + +#include "smith/infrastructure/application_manager.hpp" +#include "smith/smith_config.hpp" +#include "smith/mesh_utils/mesh_utils_base.hpp" +#include "smith/numerics/stdfunction_operator.hpp" +#include "smith/numerics/functional/functional.hpp" +#include "smith/numerics/functional/tensor.hpp" + +using namespace smith; +using namespace smith::profiling; + +template +void L2_periodic_index_test(mfem::Element::Type element_type) +{ + using test_space = L2; + using trial_space = L2; + + auto initial_mesh = mfem::Mesh(mfem::Mesh::MakeCartesian3D(4, 4, 4, element_type, 1.0, 1.0, 1.0)); + + mfem::Vector x_translation({1.0, 0.0, 0.0}); + mfem::Vector y_translation({0.0, 1.0, 0.0}); + mfem::Vector z_translation({0.0, 0.0, 1.0}); + std::vector translations = {x_translation, y_translation, z_translation}; + double tol = 1e-6; + + std::vector periodicMap = initial_mesh.CreatePeriodicVertexMapping(translations, tol); + + auto mesh = mesh::refineAndDistribute(mfem::Mesh::MakePeriodic(initial_mesh, periodicMap)); + + auto [test_fespace, test_fec] = smith::generateParFiniteElementSpace(mesh.get()); + auto [trial_fespace, trial_fec] = smith::generateParFiniteElementSpace(mesh.get()); + + // Initialize the ParGridFunction by dof coordinates + mfem::ParGridFunction U_gf(trial_fespace.get()); + mfem::VectorFunctionCoefficient vcoef(dim, [](const mfem::Vector& X, mfem::Vector& F) { + int d = X.Size(); + F.SetSize(d); + for (int i = 0; i < d; ++i) { + F(i) = X(i); + } + }); + U_gf.ProjectCoefficient(vcoef); + + mfem::Vector U(trial_fespace->TrueVSize()); + U_gf.GetTrueDofs(U); + + // Construct the new functional object using the specified test and trial spaces + Functional residual(test_fespace.get(), {trial_fespace.get()}); + + Domain periodic_faces = Domain::ofInteriorFaces(*mesh, by_attr({1, 2, 3, 4, 5, 6})); + + // Define the integral of jumps over all interior faces + residual.AddInteriorFaceIntegral( + Dimension{}, DependsOn<0>{}, + [=](double /*t*/, auto X, auto velocity) { + // compute the surface normal + auto dX_dxi = get(X); + auto n = normalize(cross(dX_dxi)); + + // extract the velocity values from each side of the interface + // note: the orientation convention is such that the normal + // computed as above will point from from side 1->2 + auto [u_1, u_2] = velocity; + SLIC_INFO(axom::fmt::format("One side = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), + axom::fmt::streamed(u_2), axom::fmt::streamed(u_1 - u_2))); + + auto a = dot(u_1 - u_2, n) - 1.0; + + auto f_1 = u_1 * a; + auto f_2 = u_2 * a; + return smith::tuple{f_1, f_2}; + }, + periodic_faces); + + double t = 0.0; + + auto value = residual(t, U); + EXPECT_NEAR(0., value.Norml2(), 1.e-12); +} + +TEST(periodic_index, L2_test_tets_linear) +{ + L2_periodic_index_test<3, 1>(mfem::Element::Type::TETRAHEDRON); +} + +TEST(periodic_index, L2_test_hex_linear) +{ + L2_periodic_index_test<3, 1>(mfem::Element::Type::HEXAHEDRON); +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/numerics/functional/typedefs.hpp b/src/smith/numerics/functional/typedefs.hpp index 1d4972fd9f..68e0a32e96 100644 --- a/src/smith/numerics/functional/typedefs.hpp +++ b/src/smith/numerics/functional/typedefs.hpp @@ -8,10 +8,6 @@ namespace smith { -// sam: this is a kludge-- it looks like the DG spaces need to use some interface defined ONLY on -// mfem::ParMesh / mfeme::ParFiniteElementSpace, but I'm not ready to pull the trigger on a big -// interface change like that, so these typedefs mark the parts that would need to eventually change - /// @cond using mesh_t = mfem::ParMesh; using fes_t = mfem::ParFiniteElementSpace; From d49971dde7a95e79432c8565044b8118a0df8182 Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 20 Mar 2026 18:33:35 -0700 Subject: [PATCH 2/6] fix style --- src/smith/numerics/functional/geometric_factors.cpp | 4 +++- .../numerics/functional/tests/dg_periodic_bc_index.cpp | 10 ++-------- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/smith/numerics/functional/geometric_factors.cpp b/src/smith/numerics/functional/geometric_factors.cpp index e97fd59365..46dd10b79b 100644 --- a/src/smith/numerics/functional/geometric_factors.cpp +++ b/src/smith/numerics/functional/geometric_factors.cpp @@ -93,7 +93,9 @@ GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry:: for (uint64_t j = 0; j < restriction.nodes_per_elem; j++) { // ignore the coordinate dofs on the other side of the face // coordinate on both side of the face should have the same values - if (j >= new_nodes_per_elem) { continue; } + if (j >= new_nodes_per_elem) { + continue; + } uint64_t E_id = (i * restriction.components + c) * restriction.nodes_per_elem + j; X_h1[H1_id++] = X_e[int(E_id)]; diff --git a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp index acb7e3293d..8a4e8819d0 100644 --- a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp +++ b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp @@ -91,15 +91,9 @@ void L2_periodic_index_test(mfem::Element::Type element_type) EXPECT_NEAR(0., value.Norml2(), 1.e-12); } -TEST(periodic_index, L2_test_tets_linear) -{ - L2_periodic_index_test<3, 1>(mfem::Element::Type::TETRAHEDRON); -} +TEST(periodic_index, L2_test_tets_linear) { L2_periodic_index_test<3, 1>(mfem::Element::Type::TETRAHEDRON); } -TEST(periodic_index, L2_test_hex_linear) -{ - L2_periodic_index_test<3, 1>(mfem::Element::Type::HEXAHEDRON); -} +TEST(periodic_index, L2_test_hex_linear) { L2_periodic_index_test<3, 1>(mfem::Element::Type::HEXAHEDRON); } int main(int argc, char* argv[]) { From 79ff7ffba8419697485b80ce3b6712853f9480c1 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 25 Mar 2026 15:54:58 -0700 Subject: [PATCH 3/6] correct local data size --- .../numerics/functional/geometric_factors.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/smith/numerics/functional/geometric_factors.cpp b/src/smith/numerics/functional/geometric_factors.cpp index 46dd10b79b..5283170eb4 100644 --- a/src/smith/numerics/functional/geometric_factors.cpp +++ b/src/smith/numerics/functional/geometric_factors.cpp @@ -78,7 +78,20 @@ GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry:: auto restriction = smith::ElementRestriction(fes, geom, element_ids); mfem::Vector X_e(int(restriction.ESize())); - restriction.Gather(*nodes, X_e); + + // Since nodes are in DG space for periodic meshes, we need a correctly sized local data vector for element + // restriction operator. Unlike the preparation operation in functional (appendFaceNbrData), we don't need + // to fill the entries in [ nodes->Size(), nodes->Size() + nodes->FaceNbrData().Size() ) range. This is because + // we only need one set of coordinates to compute geometric factor in H1 interpolation and we will only use + // the set local to the processor. + if (domain.type_ == Domain::Type::InteriorFaces && fes->IsDGSpace()) { + mfem::Vector nodes_L(nodes->Size() + nodes->FaceNbrData().Size()); + nodes_L.SetVector(*nodes, 0); + + restriction.Gather(nodes_L, X_e); + } else { + restriction.Gather(*nodes, X_e); + } // For periodic meshes, the mesh nodes are in DG space, which will double the nodes_per_elem for interior faces. // Therefore, we must discard half of the entries to recover H1 coordinates and correctly compute geometric factor. From 827a6a8e96f60cac19e45c66090768ca3f7d4367 Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 21 Apr 2026 09:34:49 -0700 Subject: [PATCH 4/6] fix 2d mesh edge orientation --- .../functional/element_restriction.cpp | 41 +++++++++------- .../functional/tests/dg_periodic_bc_index.cpp | 49 +++++++++++++++---- 2 files changed, 61 insertions(+), 29 deletions(-) diff --git a/src/smith/numerics/functional/element_restriction.cpp b/src/smith/numerics/functional/element_restriction.cpp index 1b1fee9e86..db15820548 100644 --- a/src/smith/numerics/functional/element_restriction.cpp +++ b/src/smith/numerics/functional/element_restriction.cpp @@ -502,15 +502,9 @@ axom::Array GetFaceDofs(const smith::f if (mesh->Dimension() == 2) { mesh->GetElementEdges(elem, elem_side_ids, orientations); - // mfem returns {-1, 1} for edge orientations, - // but {0, 1, ... , n} for face orientations. - // Here, we renumber the edge orientations to - // {0 (no permutation), 1 (reversed)} so the values can be - // consistently used as indices into a permutation table - for (auto& o : orientations) { - o = (o == -1) ? 1 : 0; - } - + // mfem returns {-1, 1} for edge orientations, but {0, 1, ... , n} for face orientations. + // Therefore, for 2D meshes, we discard orientations obtained by GetElementEdges and + // use orientations from GetFaceInformation } else { mesh->GetElementFaces(elem, elem_side_ids, orientations); } @@ -528,13 +522,23 @@ axom::Array GetFaceDofs(const smith::f mfem::Geometry::Type elem_geom = mesh->GetElementGeometry(elem); - // mfem uses different conventions for boundary element orientations in 2D and 3D. - // In 2D, mfem's official edge orientations on the boundary will always be a mix of - // CW and CCW, so we have to discard mfem's orientation information in order - // to get a consistent winding. + // mfem uses different conventions for interface element (edge/face) orientations in 2D and 3D. + // In 2D, mfem's official edge orientations will always be a mix of CW and CCW, so we have to discard + // mfem's orientation information in order to get a consistent winding. // - // In 3D, mfem does use a consistently CCW winding for boundary faces (I think). - int orientation = (mesh->Dimension() == 2 && info.IsBoundary()) ? 0 : orientations[i]; + // In 3D, mfem does use a consistently CCW winding for faces. + int orientation = -1; + if (mesh->Dimension() == 2) { + if (elem == info.element[0].index) { + orientation = info.element[0].orientation; + } else if (elem == info.element[1].index) { + orientation = info.element[1].orientation; + } + } else if (mesh->Dimension() == 3) { + orientation = orientations[i]; + } else { + SLIC_ERROR("face orientation is not determined."); + } // 4. extract only the dofs that correspond to side `i` for (auto k : face_perm(orientation)) { @@ -583,10 +587,9 @@ axom::Array GetFaceDofs(const smith::f int ghost_orientation; if (mesh->Dimension() == 2) { // In 2D, orientations[i] sometimes is inverted as compared to the ones from GetFaceInformation - // In this case, we stay with the orientations constructed previously by GetElementEdges, which is strictly - // CCW So we set ghost orientation to be exactly the opposite of the original orientation (orientations[i]) - // side note: even if you use info.element[1].orientation, looks like it's still fine. - // The only difference is the order those face dof indices gets registered in std::vector face_dofs + // In this case, we previous discard orientations obtained by GetElementEdges and use the ones from + // GetFaceInformation, which is strictly CCW. Then, we set ghost orientation to be exactly the opposite + // of the original orientation ghost_orientation = (orientation == 0) ? 1 : 0; } else { // In 3D, I think it's consistently CCW. So we directly use the orientation from GetFaceInformation diff --git a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp index 00820bc131..f2f8cd21b2 100644 --- a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp +++ b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp @@ -22,20 +22,40 @@ using namespace smith; using namespace smith::profiling; +// This test is only for translational periodic meshes +// At xmax, the mfem mesh nodal coordinates are (xmax, y) | (xmin, y). Therefore, u1 - u2 = xmax - xmin = L +// and n = {1, 0}. Consequently, dot(u1 - u2, n) = L. The same logic applies to other periodic boundaries. template -void L2_periodic_index_test(mfem::Element::Type element_type) +void L2_periodic_index_test(std::string meshfile) { using test_space = L2; using trial_space = L2; - auto initial_mesh = mfem::Mesh(mfem::Mesh::MakeCartesian3D(4, 4, 4, element_type, 1.0, 1.0, 1.0)); + auto initial_mesh = buildMeshFromFile(meshfile); + initial_mesh.UniformRefinement(); + initial_mesh.UniformRefinement(); + + // 2D patch mesh has side length of 1.0 and 3D patch mesh has side length of 2.0 + std::vector translations; + if (dim == 2) { + mfem::Vector x_translation({1.0, 0.0}); + mfem::Vector y_translation({0.0, 1.0}); + + translations.push_back(x_translation); + translations.push_back(y_translation); + } else if (dim == 3) { + mfem::Vector x_translation({2.0, 0.0, 0.0}); + mfem::Vector y_translation({0.0, 2.0, 0.0}); + mfem::Vector z_translation({0.0, 0.0, 2.0}); + + translations.push_back(x_translation); + translations.push_back(y_translation); + translations.push_back(z_translation); + } else { + SLIC_ERROR_ROOT("This test is only for 2D and 3D meshes"); + } - mfem::Vector x_translation({1.0, 0.0, 0.0}); - mfem::Vector y_translation({0.0, 1.0, 0.0}); - mfem::Vector z_translation({0.0, 0.0, 1.0}); - std::vector translations = {x_translation, y_translation, z_translation}; double tol = 1e-6; - std::vector periodicMap = initial_mesh.CreatePeriodicVertexMapping(translations, tol); auto mesh = mesh::refineAndDistribute(mfem::Mesh::MakePeriodic(initial_mesh, periodicMap)); @@ -77,7 +97,8 @@ void L2_periodic_index_test(mfem::Element::Type element_type) SLIC_INFO(axom::fmt::format("One side = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), axom::fmt::streamed(u_2), axom::fmt::streamed(u_1 - u_2))); - auto a = dot(u_1 - u_2, n) - 1.0; + // The (dim - 1.0) is because 2D patch mesh has side length of 1.0 and 3D patch mesh has side length of 2.0 + auto a = dot(u_1 - u_2, n) - (dim - 1.0); auto f_1 = u_1 * a; auto f_2 = u_2 * a; @@ -91,9 +112,17 @@ void L2_periodic_index_test(mfem::Element::Type element_type) EXPECT_NEAR(0., value.Norml2(), 1.e-12); } -TEST(periodic_index, L2_test_tets_linear) { L2_periodic_index_test<3, 1>(mfem::Element::Type::TETRAHEDRON); } +TEST(periodic_index, L2_test_tris_linear) { L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); } +TEST(periodic_index, L2_test_tris_quadratic) { L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); } + +TEST(periodic_index, L2_test_quads_linear) { L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); } +TEST(periodic_index, L2_test_quads_quadratic) { L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); } + +TEST(periodic_index, L2_test_tets_linear) { L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } +TEST(periodic_index, L2_test_tets_quadratic) { L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } -TEST(periodic_index, L2_test_hex_linear) { L2_periodic_index_test<3, 1>(mfem::Element::Type::HEXAHEDRON); } +TEST(periodic_index, L2_test_hex_linear) { L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } +TEST(periodic_index, L2_test_hex_quadratic) { L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } int main(int argc, char* argv[]) { From a037f9f414f5c408e54b0cd2e658a2e9f736d64c Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 21 Apr 2026 10:42:02 -0700 Subject: [PATCH 5/6] fix style --- .../functional/element_restriction.cpp | 2 +- .../functional/tests/dg_periodic_bc_index.cpp | 40 +++++++++++++++---- 2 files changed, 33 insertions(+), 9 deletions(-) diff --git a/src/smith/numerics/functional/element_restriction.cpp b/src/smith/numerics/functional/element_restriction.cpp index db15820548..db34ef2f07 100644 --- a/src/smith/numerics/functional/element_restriction.cpp +++ b/src/smith/numerics/functional/element_restriction.cpp @@ -587,7 +587,7 @@ axom::Array GetFaceDofs(const smith::f int ghost_orientation; if (mesh->Dimension() == 2) { // In 2D, orientations[i] sometimes is inverted as compared to the ones from GetFaceInformation - // In this case, we previous discard orientations obtained by GetElementEdges and use the ones from + // In this case, we previously discarded orientations obtained by GetElementEdges and use the ones from // GetFaceInformation, which is strictly CCW. Then, we set ghost orientation to be exactly the opposite // of the original orientation ghost_orientation = (orientation == 0) ? 1 : 0; diff --git a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp index f2f8cd21b2..f7b96e58d2 100644 --- a/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp +++ b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp @@ -112,17 +112,41 @@ void L2_periodic_index_test(std::string meshfile) EXPECT_NEAR(0., value.Norml2(), 1.e-12); } -TEST(periodic_index, L2_test_tris_linear) { L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); } -TEST(periodic_index, L2_test_tris_quadratic) { L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); } +TEST(periodic_index, L2_test_tris_linear) +{ + L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); +} +TEST(periodic_index, L2_test_tris_quadratic) +{ + L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_tris.mesh"); +} -TEST(periodic_index, L2_test_quads_linear) { L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); } -TEST(periodic_index, L2_test_quads_quadratic) { L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); } +TEST(periodic_index, L2_test_quads_linear) +{ + L2_periodic_index_test<2, 1>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); +} +TEST(periodic_index, L2_test_quads_quadratic) +{ + L2_periodic_index_test<2, 2>(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"); +} -TEST(periodic_index, L2_test_tets_linear) { L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } -TEST(periodic_index, L2_test_tets_quadratic) { L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } +TEST(periodic_index, L2_test_tets_linear) +{ + L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); +} +TEST(periodic_index, L2_test_tets_quadratic) +{ + L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"); +} -TEST(periodic_index, L2_test_hex_linear) { L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } -TEST(periodic_index, L2_test_hex_quadratic) { L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } +TEST(periodic_index, L2_test_hex_linear) +{ + L2_periodic_index_test<3, 1>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); +} +TEST(periodic_index, L2_test_hex_quadratic) +{ + L2_periodic_index_test<3, 2>(SMITH_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); +} int main(int argc, char* argv[]) { From 0d1f4a662cddc39af7ed417d3299a7cf0ab13a41 Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 1 May 2026 09:20:55 -0700 Subject: [PATCH 6/6] make sure we can impose Dirichlet BC on internal BCs (eg. periodic BC) --- src/smith/numerics/functional/domain.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/numerics/functional/domain.cpp b/src/smith/numerics/functional/domain.cpp index a80a97b543..2c5a18a575 100644 --- a/src/smith/numerics/functional/domain.cpp +++ b/src/smith/numerics/functional/domain.cpp @@ -473,7 +473,7 @@ mfem::Array Domain::dof_list(const smith::fes_t* fes) const GetDofs = [&](int i, mfem::Array& vdofs) { return fes->GetElementDofs(i, vdofs); }; } - if (type_ == Type::BoundaryElements) { + if (type_ == Type::BoundaryElements || type_ == Type::InteriorBoundaryElements) { GetDofs = [&](int i, mfem::Array& vdofs) { return fes->GetFaceDofs(i, vdofs); }; }