Li/periodic geometric factor

src/smith/numerics/functional/domain.cpp

@@ -473,7 +473,7 @@ mfem::Array<int> Domain::dof_list(const smith::fes_t* fes) const
     GetDofs = [&](int i, mfem::Array<int>& vdofs) { return fes->GetElementDofs(i, vdofs); };
   }
 
-  if (type_ == Type::BoundaryElements) {
+  if (type_ == Type::BoundaryElements || type_ == Type::InteriorBoundaryElements) {
     GetDofs = [&](int i, mfem::Array<int>& vdofs) { return fes->GetFaceDofs(i, vdofs); };
   }
 

src/smith/numerics/functional/element_restriction.cpp

@@ -502,15 +502,9 @@ axom::Array<DoF, 2, smith::detail::host_memory_space> 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<DoF, 2, smith::detail::host_memory_space> 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<DoF, 2, smith::detail::host_memory_space> 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

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<const mfem::ParGridFunction*>(domain.mesh_.GetNodes());
   mfem::ParFiniteElementSpace* fes = nodes->ParFESpace();
-  // const mfem::GridFunction* nodes = domain.mesh_.GetNodes();
-  // const mfem::FiniteElementSpace* fes = nodes->FESpace();
 
   const std::vector<int>& 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);

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}

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);

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 <fstream>
+#include <iostream>
+
+#include "mfem.hpp"
+
+#include <gtest/gtest.h>
+
+#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 <int dim, int p>
+void L2_periodic_index_test(std::string meshfile)
+{
+  using test_space = L2<p, dim>;
+  using trial_space = L2<p, dim>;
+
+  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<mfem::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<int> periodicMap = initial_mesh.CreatePeriodicVertexMapping(translations, tol);
+
+  auto mesh = mesh::refineAndDistribute(mfem::Mesh::MakePeriodic(initial_mesh, periodicMap));
+
+  auto [test_fespace, test_fec] = smith::generateParFiniteElementSpace<test_space>(mesh.get());
+  auto [trial_fespace, trial_fec] = smith::generateParFiniteElementSpace<trial_space>(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<test_space(trial_space)> residual(test_fespace.get(), {trial_fespace.get()});
+
+  Domain periodic_faces = Domain::ofInteriorBoundaryElements(*mesh, by_attr<dim>({1, 2, 3, 4, 5, 6}));
+
+  // Define the integral of jumps over all interior faces
+  residual.AddInteriorFaceIntegral(
+      Dimension<dim - 1>{}, DependsOn<0>{},
+      [=](double /*t*/, auto X, auto velocity) {
+        // compute the surface normal
+        auto dX_dxi = get<DERIVATIVE>(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();
+}

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;