diff --git a/src/smith/numerics/functional/domain.cpp b/src/smith/numerics/functional/domain.cpp index a80a97b54..2c5a18a57 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); }; } diff --git a/src/smith/numerics/functional/element_restriction.cpp b/src/smith/numerics/functional/element_restriction.cpp index 1b1fee9e8..db34ef2f0 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 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; } else { // In 3D, I think it's consistently CCW. So we directly use the orientation from GetFaceInformation diff --git a/src/smith/numerics/functional/geometric_factors.cpp b/src/smith/numerics/functional/geometric_factors.cpp index ae8e5cf6d..05d971ff6 100644 --- a/src/smith/numerics/functional/geometric_factors.cpp +++ b/src/smith/numerics/functional/geometric_factors.cpp @@ -73,14 +73,52 @@ 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); 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::InteriorBoundaryElements && 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. + // Note that faces on periodic boundaries in mfem mesh are considered interior faces with boundary attributes. + if (domain.type_ == Domain::Type::InteriorBoundaryElements && 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 9bc329d78..8a09f75d0 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 a58dd0dce..e7a7aea30 100644 --- a/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/smith/numerics/functional/tests/dg_ghost_face_index.cpp @@ -76,7 +76,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(std::format("One size = {}, The other side = {}, Jump = {}", smith::format::streamed(u_1), + SLIC_INFO(std::format("One side = {}, The other side = {}, Jump = {}", smith::format::streamed(u_1), smith::format::streamed(u_2), smith::format::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 000000000..f7b96e58d --- /dev/null +++ b/src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp @@ -0,0 +1,156 @@ + +// 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; + +// 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(std::string meshfile) +{ + using test_space = L2; + using trial_space = L2; + + 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"); + } + + 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::ofInteriorBoundaryElements(*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))); + + // 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; + 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_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>(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[]) +{ + ::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 1d4972fd9..68e0a32e9 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;