Li/periodic geometric factor

src/smith/numerics/functional/domain.cpp

@@ -375,6 +375,80 @@ Domain Domain::ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::ve
   return domain_of_boundary_elems<3>(mesh, func);
 }
 
+/// @brief constructs a domain from some subset of the interior face elements in a mesh
+template <int d>
+Domain domain_of_interior_faces(const mesh_t& mesh, std::function<bool(std::vector<tensor<double, d>>, int)> predicate)
+{
+  assert(mesh.SpaceDimension() == d);
+  Domain output{mesh, d - 1, Domain::Type::InteriorFaces};
+
+  mfem::Array<int> 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<int> vertex_ids;
+    mesh.GetFaceVertices(f, vertex_ids);
+
+    auto x = gather<d>(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<bool(std::vector<vec2>, int)> func)
+{
+  return domain_of_interior_faces<2>(mesh, func);
+}
+
+Domain Domain::ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec3>, 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
  *

src/smith/numerics/functional/domain.hpp

@@ -167,6 +167,16 @@ struct Domain {
   /// @overload
   static Domain ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec3>, 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<bool(std::vector<vec2>, int)> func);
+
+  /// @overload
+  static Domain ofInteriorFaces(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
+
   /// @brief get elements by geometry type
   const std::vector<int>& get(mfem::Geometry::Type geom) const
   {

src/smith/numerics/functional/geometric_factors.cpp

@@ -73,15 +73,40 @@ 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);
 
+  // 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);
 

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

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

src/smith/numerics/functional/tests/dg_periodic_bc_index.cpp

@@ -0,0 +1,103 @@
+
+// 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;
+
+template <int dim, int p>
+void L2_periodic_index_test(mfem::Element::Type element_type)
+{
+  using test_space = L2<p, dim>;
+  using trial_space = L2<p, dim>;
+
+  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<mfem::Vector> translations = {x_translation, y_translation, z_translation};
+  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::ofInteriorFaces(*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)));
+
+        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();
+}

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;