Fix projection onto NC internal faces

[hughcars]

This PR fixes an error where if a boundary condition is set on an internal boundary with nonconformal refinement, face dofs would not be correctly accounted for on the boundary face. The fix consists of using bdr_faces from the closure method.

Can be seen by setting the field to 0 in ex1p, then project 1 onto the internal boundary (not setting others) using the DividingPlaneMesh mesh used in test_ncmesh.cpp.

Before the fix:

After the fix:

PR	Author	Editor	Reviewers	Assignment	Approval	Merge
#4259	@hughcars	@tzanio	@psocratis + @dylan-copeland	5/11/24	⌛due 5/25/24	⌛due 6/1/24

PR Checklist

• Code builds.
• Code passes make style.
• Update CHANGELOG:
  • Is this a new feature users need to be aware of? New or updated example or miniapp?
  • Does it make sense to create a new section in the CHANGELOG to group with other related features?
• Update INSTALL:
  • Had a new optional library been added? If so, what range of versions of this library are required? (Make sure the external library is compatible with our BSD license, e.g. it is not licensed under GPL!)
  • Have the version ranges for any required or optional libraries changed?
  • Does make or cmake have a new target?
  • Did the requirements or the installation process change? (rare)
• Update continuous integration server configurations if necessary (e.g. with new version requirements for each of MFEM's dependencies)
  • .github
  • .appveyor.yml
• Update .gitignore:
  • Check if make distclean; git status shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
  • Add new patterns (just for the new files above) and re-run the above test.
• New examples:
  • All sample runs at the top of the example source file work.
  • Update examples/makefile:
    • Add the example code to the appropriate SEQ_EXAMPLES and PAR_EXAMPLES variables.
    • Add any files generated by it to the clean target.
    • Add the example binary and any files generated by it to the top-level .gitignore file.
  • Update examples/CMakeLists.txt:
    • Add the example code to the ALL_EXE_SRCS variable.
    • Make sure THIS_TEST_OPTIONS is set correctly for the new example.
  • List the new example in doc/CodeDocumentation.dox.
  • If new examples directory (e.g.examples/pumi), list it in doc/CodeDocumentation.conf.in
  • Companion pull request for documentation in mfem/web repo:
    • Update or add example-specific documentation, see e.g. the src/examples.md.
    • Add the description, labels and screenshots in src/examples.md and src/img.
    • In examples.md, list the example under the appropriate categories, add new categories if necessary.
    • Add a short description of the example in the "Extensive Examples" section of features.md.
• New miniapps:
  • All sample runs at the top of the miniapp source file work.
  • Update top-level makefile and makefile in corresponding miniapp directory.
  • Add the miniapp binary and any files generated by it to the top-level .gitignore file.
  • Update CMake build system:
    • Update the CMakeLists.txt file in the miniapps directory, if the new miniapp is in a new directory.
    • Add/update the CMakeLists.txt file in the new miniapp directory.
    • Consider adding a new test for the new miniapp.
  • List the new miniapp in doc/CodeDocumentation.dox
  • If new miniapps directory (e.g.miniapps/nurbs), add it to MINIAPP_SUBDIRS in the makefile.
  • If new miniapps directory (e.g.miniapps/nurbs), list it in doc/CodeDocumentation.conf.in
  • Companion pull request for documentation in mfem/web repo:
    • Update or add miniapp-specific documentation, see e.g. the src/meshing.md and src/electromagnetics.md files.
    • Add the description, labels and screenshots in src/examples.md and src/img.
    • The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
    • Add a short description of the miniapp in the "Extensive Examples" section of features.md.
• New capability:
  • All new public, protected, and private classes, methods, data members, and functions have full Doxygen-style documentation in source comments. Documentation should include descriptions of member data, function arguments and return values, template parameters, and prerequisites for calling new functions.
  • Pointer arguments and return values must specify whether ownership is being transferred or lent with the call.
  • Any new functions should include descriptions of their intended use e.g. for internal use only, user-facing, etc., along with references to example code whenever possible/appropriate.
  • Consider adding new sample runs in existing examples to highlight the new capability.
  • Consider saving cool simulation pictures with the new capability in the Confluence gallery (LLNL only) or submitting them, via pull request, to the gallery section of the mfem/web repo.
  • If this is a major new feature, consider mentioning it in the short summary inside README (rare).
  • List major new classes in doc/CodeDocumentation.dox (rare).
• Update this checklist, if the new pull request affects it.
• Run make unittest to make sure all unit tests pass.
• Run the tests in tests/scripts.
• (LLNL only) After merging:
  • Update internal tests to include the new features.

[tzanio]

This PR is now under review (see the table in the PR description). To help with the review process, please do not force push to the branch.

[hughcars]

Testing this on ex1p can be done with

diff --git a/examples/ex1p.cpp b/examples/ex1p.cpp
index c3240d621..6d28bd392 100644
--- a/examples/ex1p.cpp
+++ b/examples/ex1p.cpp
@@ -66,6 +66,8 @@
 using namespace std;
 using namespace mfem;
 
+Mesh DividingPlaneMesh(bool tet_mesh, bool split);
+
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI and HYPRE.
@@ -128,7 +130,30 @@ int main(int argc, char *argv[])
    // 4. Read the (serial) mesh from the given mesh file on all processors.  We
    //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
    //    and volume meshes with the same code.
-   Mesh mesh(mesh_file, 1, 1);
+   Mesh mesh = DividingPlaneMesh(false, true); // set first arg to true to have tets
+   mesh.EnsureNCMesh(true);
+   auto OneSidedNCRefine = [](Mesh &mesh)
+   {
+      // Pick one element attached to the new boundary attribute and refine.
+      const auto interface_attr = mesh.bdr_attributes.Max();
+      Array<int> el_to_ref;
+      for (int nbe = 0; nbe < mesh.GetNBE(); nbe++)
+      {
+         if (mesh.GetBdrAttribute(nbe) == interface_attr)
+         {
+            int f, o, e1, e2;
+            mesh.GetBdrElementFace(nbe, &f, &o);
+            mesh.GetFaceElements(f, &e1, &e2);
+            el_to_ref.Append(e1);
+         }
+      }
+      mesh.GeneralRefinement(el_to_ref);
+      return;
+   };
+
+
+
    int dim = mesh.Dimension();
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In
@@ -136,13 +161,15 @@ int main(int argc, char *argv[])
    //    'ref_levels' to be the largest number that gives a final mesh with no
    //    more than 10,000 elements.
    {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
+      int ref_levels = 2;
       for (int l = 0; l < ref_levels; l++)
       {
          mesh.UniformRefinement();
       }
    }
+   OneSidedNCRefine(mesh);
 
    // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
@@ -150,7 +177,7 @@ int main(int argc, char *argv[])
    ParMesh pmesh(MPI_COMM_WORLD, mesh);
    mesh.Clear();
    {
-      int par_ref_levels = 2;
+      int par_ref_levels = 0;
       for (int l = 0; l < par_ref_levels; l++)
       {
          pmesh.UniformRefinement();
@@ -193,10 +220,11 @@ int main(int argc, char *argv[])
    //    by marking all the boundary attributes from the mesh as essential
    //    (Dirichlet) and converting them to a list of true dofs.
    Array<int> ess_tdof_list;
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
    if (pmesh.bdr_attributes.Size())
    {
-      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-      ess_bdr = 1;
+      ess_bdr = 0;
+      ess_bdr.Last() = 1;
       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
 
@@ -214,6 +242,8 @@ int main(int argc, char *argv[])
    ParGridFunction x(&fespace);
    x = 0.0;
 
+   x.ProjectBdrCoefficient(one, ess_bdr);
+
    // 11. Set up the parallel bilinear form a(.,.) on the finite element space
    //     corresponding to the Laplacian operator -Delta, by adding the
    //     Diffusion domain integrator.
@@ -310,3 +340,45 @@ int main(int argc, char *argv[])
 
    return 0;
 }
+
+
+
+
+Mesh DividingPlaneMesh(bool tet_mesh, bool split)
+{
+   auto mesh = Mesh("../../data/ref-cube.mesh");
+   {
+      Array<Refinement> refs;
+      refs.Append(Refinement(0, Refinement::X));
+      mesh.GeneralRefinement(refs);
+   }
+   delete mesh.ncmesh;
+   mesh.ncmesh = nullptr;
+   mesh.FinalizeTopology();
+   mesh.Finalize(true, true);
+
+   mesh.SetAttribute(0, 1);
+   mesh.SetAttribute(1, split ? 2 : 1);
+
+   // Introduce internal boundary elements
+   const int new_attribute = mesh.bdr_attributes.Max() + 1;
+   for (int f = 0; f < mesh.GetNumFaces(); ++f)
+   {
+      int e1, e2;
+      mesh.GetFaceElements(f, &e1, &e2);
+      if (e1 >= 0 && e2 >= 0 && mesh.GetAttribute(e1) != mesh.GetAttribute(e2))
+      {
+         // This is the internal face between attributes.
+         auto *new_elem = mesh.GetFace(f)->Duplicate(&mesh);
+         new_elem->SetAttribute(new_attribute);
+         mesh.AddBdrElement(new_elem);
+      }
+   }
+   if (tet_mesh)
+   {
+      mesh = Mesh::MakeSimplicial(mesh);
+   }
+   mesh.FinalizeTopology();
+   mesh.Finalize(true, true);
+   return mesh;
+}
\ No newline at end of file

then run with mpirun -np 2 ex1p -o 2 which should show the impact on internal boundary projection. Changing the first argument to DividingPlaneMesh to true will give tets and a similar result.

[dylan-copeland]

Extra space after the comma. Does coef need to be nonzero for this test to make sense? If so, this could be in a short comment.

[dylan-copeland]

Should there be a parallel test? It seems the issue is fixed in parallel, considering the ex1p test you provided. Do you think the serial test is sufficient, since the changes are only in gridfunc.cpp?