improved some comments

.gitignore

@@ -229,6 +229,7 @@ miniapps/meshing/reflected.mesh
 miniapps/meshing/optimized*
 miniapps/meshing/perturbed*
 miniapps/meshing/polar-nc.mesh
+miniapps/meshing/ext-mesh-mapping
 
 miniapps/mtop/parheat
 miniapps/mtop/ParHeat*

mesh/ncmesh.hpp

@@ -153,6 +153,11 @@ public:
    int GetNFaces() const { return NFaces; }
    virtual int GetNGhostElements() const { return 0; }
 
+   /** NCMesh can change the vertex ordering after refinement, coarsening, or on
+       creation of the NCMesh object.  After update operation the Vertex ID Map
+       contains the remapping information. */
+   const Array<int> &GetVertexIDMap() {return vertex_nodeId;}
+
    /** Perform the given batch of refinements. Please note that in the presence
        of anisotropic splits additional refinements may be necessary to keep
        the mesh consistent. However, the function always performs at least the
@@ -282,7 +287,6 @@ public:
       }
    }
 
-
    // coarse/fine transforms
 
    /** Remember the current layer of leaf elements before the mesh is refined.

miniapps/meshing/CMakeLists.txt

@@ -108,6 +108,10 @@ if (MFEM_USE_MPI)
     MAIN pminimal-surface.cpp
     LIBRARIES mfem)
 
+  add_mfem_miniapp(ext-mesh-mapping
+    MAIN ext-mesh-mapping.cpp
+    LIBRARIES mfem)
+
   # Add parallel tests.
   if (MFEM_ENABLE_TESTING)
     add_test(NAME pmesh-optimizer_np=${MFEM_MPI_NP}

miniapps/meshing/ext-mesh-mapping.cpp

@@ -0,0 +1,317 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//             ---------------------------------------------------
+//                       External Mesh Mapping Miniapp
+//             ---------------------------------------------------
+//
+// This miniapp starts with a serial non-conforming mesh in a dummy format and
+// demonstrates how to build up a corresponding non-conforming MFEM Mesh and
+// then decompose it into a parallel ParMesh.  As part of this process we will
+// demonstrate how to obtain and compose the vertex ID mappings from the
+// various steps that can shuffle the vertices.  In the end this will let us
+// map between the vertex ID numbers in the external dummy mesh and the 
+// parallel non-conforming mesh constructed in MFEM.  This is the sort of
+// map you will need to maintain if you intend to add MFEM to an existing 
+// simulation code and share data between the MFEM mesh and the  and the 
+// field data in  and need them to exist and share data with other kinds of meshes in that
+// code.  If you are starting a new MFEM code, it makes much more sense to do everything with
+// MFEM meshes.
+//
+// Compile with: make ext-mesh-mapping
+//
+// Sample runs:   mpirun -np 4 ext-mesh-mapping
+
+#include "mfem.hpp"
+#include "ext-mesh-mapping.hpp"
+
+using namespace mfem;
+
+Mesh *build_mfem_mesh(DummyMesh *dmesh);
+void create_pmesh_to_mesh_emaps(Array<int> &partition, ParMesh *pmesh,
+                                Array<int> &emap);
+void create_pmesh_to_mesh_vmaps(ParMesh *pmesh, Mesh *mesh, Array<int> &emap,
+                                Array<int> &vmap);
+
+void print_dmesh_verts(DummyMesh *dmesh);
+void print_mesh_verts(Mesh *mesh, const Array<int> &vmap = Array<int>());
+void print_pmesh_verts(ParMesh *pmesh, const Array<int> &vmap = Array<int>());
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // This miniapp is designed to work on the test mesh with exactly 4 ranks
+   if (Mpi::WorldSize() != 4)
+   {
+      if (Mpi::Root())
+      {
+         std::cout << "Miniapp must run with exactly 4 MPI ranks." << std::endl;
+      }
+      exit(0);
+   }
+
+   // Initilize our dummy mesh and display the coordinates of the vertices
+   DummyMesh *dmesh = new DummyMesh();
+   print_dmesh_verts(dmesh);
+
+   // Now build the mfem Mesh object with all of the elements in it on all
+   // processors and capture the vertex id mapping that occurs when
+   // non-conforming meshes are finalized.
+   Mesh * mesh = build_mfem_mesh(dmesh);
+   Array<int> mesh_to_dmesh_vmap;
+   const Array<int> vmap = mesh->ncmesh->GetVertexIDMap();
+   mesh_to_dmesh_vmap = vmap;
+
+   // Print the vertices reordered using the vmap
+   print_mesh_verts(mesh, mesh_to_dmesh_vmap);
+
+   // Now enable parallel, given the following partition of the elements.
+   // Note that we only have local vertex ids in the pmesh object.
+   // The partition is the MPI Rank that each element lives on.
+   Array<int> partition({3,3,1,2,0});
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh, partition.GetData());
+
+   // Compute the mappings between the local element ids on each processor
+   // in pmesh and the original global IDs in the mesh.  This can be
+   // computed from the element partition array.
+   Array<int> pmesh_to_mesh_emap;
+   create_pmesh_to_mesh_emaps(partition, pmesh, pmesh_to_mesh_emap);
+
+   // Now compute the mappings between the local vertex ids on each
+   // processor in pmesh and the global ids in mesh.
+   Array<int> pmesh_to_mesh_vmap;
+   create_pmesh_to_mesh_vmaps(pmesh, mesh, pmesh_to_mesh_emap, pmesh_to_mesh_vmap);
+
+   // Now compose the maps to create a final map between the local vertex 
+   // numbering in the pmesh object on this processor and the original
+   // dmesh vertex numbering.
+   Array<int> final_vmap(pmesh->GetNV());
+   for (int local_vi = 0; local_vi < pmesh->GetNV(); ++local_vi)
+   {
+      final_vmap[local_vi] = mesh_to_dmesh_vmap[pmesh_to_mesh_vmap[local_vi]];
+   }
+
+   // Print the vertices with global IDs from the final_vmap
+   print_pmesh_verts(pmesh, final_vmap); 
+}
+
+// Build up an MFEM Mesh from the data in the Dummy Mesh.  Since we are
+// building the vertex and element lists in the same order as we found them
+// in the dmesh, the element id and vertex id mappings will be the 
+// identity maps.
+Mesh *build_mfem_mesh(DummyMesh *dmesh)
+{
+   // Initilize the the dimension and memory for the mesh
+   Mesh *mesh = new Mesh(2,    // The dimension of the mesh
+                         dmesh->num_vertices,
+                         dmesh->num_elements,
+                         dmesh->num_belements,
+                         2    // The dimension of the space the mesh lives in
+                        );
+
+
+   //Add the vertices into the mesh in the same order as they are in the dmesh
+   for (int vid = 0; vid < dmesh->num_vertices; ++vid)
+   {
+      mesh->AddVertex(dmesh->V[vid].x, dmesh->V[vid].y);
+   }
+
+   //Add the elements into the mesh in the same order as they are in the dmesh
+   //The dmesh is purely quads, but if it had other types we would use
+   //the other mesh.Add* methods here as well.
+   for (int eid = 0; eid < dmesh->num_elements; ++eid)
+   {
+      std::array<int,4> vid = dmesh->E[eid].vertex_ids;
+
+      //The order of vertices in the MFEM quad elements is different from the
+      //order of the vertices in out DummyMesh format so we must permute them
+      //here.  See the ref-*.mesh files in mfem/data to establish to ordering
+      //of the vertices in the MFEM elements.
+      //  Dummy    MFEM
+      //  2--3     3--2
+      //  |  |     |  |
+      //  0--1     0--1
+      mesh->AddQuad(vid[0], vid[1], vid[3], vid[2]);
+   }
+
+   //Add the boundary elements into the mesh in the same order they are in the dmesh
+   for (int bid = 0; bid < dmesh->num_belements; ++bid)
+   {
+      std::array<int,2> vid = dmesh->B[bid].vertex_ids;
+      mesh->AddBdrSegment(vid[0], vid[1]);
+   }
+
+   //Finally add vertex parents to mark element 0 for anisotropic refinement
+   for (int vpi = 0; vpi < dmesh->num_vparents; ++vpi)
+   {
+      mesh->AddVertexParents(std::get<0>(dmesh->VP[vpi]),std::get<1>(dmesh->VP[vpi]),
+                             std::get<2>(dmesh->VP[vpi]));
+   }
+
+   //This will make the mesh usable
+   mesh->FinalizeMesh();
+
+   return mesh;
+}
+
+
+// We can use the partition array to compute the mapping between the local
+// element ids on each processor in pmesh and the global element ids in mesh.
+void create_pmesh_to_mesh_emaps(Array<int> &partition, ParMesh *pmesh,
+                                Array<int> &emap)
+{
+   int my_rank = Mpi::WorldRank();
+   emap.SetSize(pmesh->GetNE());
+
+   int local_eid = 0;
+   std::vector<int> indices;
+   auto it = std::find(partition.begin(), partition.end(), my_rank);
+   while (it != partition.end())
+   {
+      emap[local_eid] = std::distance(partition.begin(), it);
+      local_eid++;
+      it++;
+      it = std::find(it, partition.end(), my_rank);
+   }
+}
+
+
+// We can use the local elements with local vertex ids defined in pmesh and the global
+// elements defined with global vertex ids in mesh to define a mapping from the
+// local vertex ids on each processor to their global vertex id numbers.
+void create_pmesh_to_mesh_vmaps(ParMesh *pmesh, Mesh *mesh, Array<int> &emap,
+                                Array<int> &vmap)
+{
+   vmap.SetSize(pmesh->GetNV());
+   for (int local_eid = 0; local_eid < pmesh->GetNE(); ++local_eid)
+   {
+      int global_eid = emap[local_eid];
+      Array<int> local_elem_verts, global_elem_verts;
+      pmesh->GetElement(local_eid)->GetVertices(local_elem_verts);
+      mesh->GetElement(global_eid)->GetVertices(global_elem_verts);
+      for (int vi = 0; vi < local_elem_verts.Size(); ++vi)
+      {
+         vmap[local_elem_verts[vi]] = global_elem_verts[vi];
+      }
+   }
+}
+
+
+void print_dmesh_verts(DummyMesh *dmesh)
+{
+   if (Mpi::Root())
+   {
+      std::cout << "6-----7-----8" << std::endl;
+      std::cout << "|     |     |" << std::endl;
+      std::cout << "|  3  |  4  |" << std::endl;
+      std::cout << "|     |     |" << std::endl;
+      std::cout << "3-----4-----5" << std::endl;
+      std::cout << "|  1  |     |" << std::endl;
+      std::cout << "9----10  2  |" << std::endl;
+      std::cout << "|  0  |     |" << std::endl;
+      std::cout << "0-----1-----2" << std::endl;
+
+      std::cout << "DummyMesh vertices:  "  << std::endl;
+      for (int vid = 0; vid < dmesh->num_vertices; ++vid)
+      {
+         std::cout << vid << ":  " << dmesh->V[vid].x << ", " << dmesh->V[vid].y <<
+                   std::endl;
+      }
+      std::cout << std::endl;
+   }
+}
+
+void print_mesh_verts(Mesh *mesh, const Array<int> &vmap)
+{
+   if (Mpi::Root())
+   {
+      std::cout << std::endl << "MFEM Mesh Vertices:  "  << std::endl;
+      if (vmap.Size() > 0)
+      {
+         std::cout << "(Remapped vertex ids)" << std::endl;
+      }
+      for (int vid = 0; vid < mesh->GetNV(); ++vid)
+      {
+         int id = vmap.Size() > 0 ? vmap[vid] : vid;
+         double *vertex = mesh->GetVertex(id);
+         std::cout << vid << ":  " <<  vertex[0] << ", " << vertex[1] << std::endl;
+      }
+      std::cout << std::endl;
+   }
+}
+
+
+void print_pmesh_verts(ParMesh *pmesh, const Array<int> &vmap)
+{
+   int my_rank = Mpi::WorldRank();
+   int num_rank = Mpi::WorldSize();
+
+   Array<int> num_verts(num_rank);
+   int my_num_verts = pmesh->GetNV();
+   MPI_Allgather(&my_num_verts, 1, MPI_INTEGER,
+                 num_verts.GetData(), 1, MPI_INTEGER, MPI_COMM_WORLD);
+
+   int max_num_verts = *std::max_element(num_verts.begin(), num_verts.end());
+   Array<int> id_data(max_num_verts);
+   Array<double> x_data(max_num_verts);
+   Array<double> y_data(max_num_verts);
+
+   //Send the data to rank 0
+   for (int vid = 0; vid < pmesh->GetNV(); ++vid)
+   {
+      int id = vmap.Size() > 0 ? vmap[vid] : vid;
+      double *vertex = pmesh->GetVertex(vid);
+      id_data[vid] = id;
+      x_data[vid] =  vertex[0];
+      y_data[vid] =  vertex[1];
+   }
+   if (my_rank != 0)
+   {
+      MPI_Send(id_data.GetData(), pmesh->GetNV(), MPI_INTEGER, 0, 0, MPI_COMM_WORLD);
+      MPI_Send(x_data.GetData(), pmesh->GetNV(), MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
+      MPI_Send(y_data.GetData(), pmesh->GetNV(), MPI_DOUBLE, 0, 2, MPI_COMM_WORLD);
+   }
+
+
+   if (my_rank == 0)
+   {
+      std::cout << "MFEM ParMesh Vertices on each processor:  "  << std::endl;
+      if (vmap.Size() > 0)
+      {
+         std::cout << "(Remapped vertex ids)" << std::endl;
+      }
+
+      for (int p = 0; p < num_rank; ++p)
+      {
+         if (p != 0)
+         {
+            MPI_Status status;
+            MPI_Recv(id_data.GetData(), num_verts[p], MPI_INTEGER, p, 0, MPI_COMM_WORLD,
+                     &status);
+            MPI_Recv(x_data.GetData(), num_verts[p], MPI_DOUBLE, p, 1, MPI_COMM_WORLD,
+                     &status);
+            MPI_Recv(y_data.GetData(), num_verts[p], MPI_DOUBLE, p, 2, MPI_COMM_WORLD,
+                     &status);
+         }
+
+         for (int vid = 0; vid < num_verts[p]; ++vid)
+         {
+            std::cout << "rank (" << p << ") id (" << id_data[vid] << "):  "
+                      << x_data[vid]  << ", " << y_data[vid] << std::endl;
+         }
+      }
+      std::cout << std::endl;
+   }
+
+   MPI_Barrier(MPI_COMM_WORLD);
+}

miniapps/meshing/ext-mesh-mapping.hpp

@@ -0,0 +1,117 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// MFEM External Mesh Mapping Miniapp
+
+#include <array>
+
+// A very simple example of an external mesh class with a minimal 
+// non-conforming mesh in it.
+// The element and vertex ids are as follows:
+// 6-----7-----8
+// |     |     |
+// |  3  |  4  |
+// |     |     |
+// 3-----4-----5
+// |  1  |     |
+// 9----10  2  |
+// |  0  |     |
+// 0-----1-----2
+class DummyMesh
+{
+   class Vertex
+   {
+   public:
+      Vertex(): x(0.0), y(0.0) {}
+
+      void Set(double x_, double y_)
+      {
+         x = x_;
+         y = y_;
+      }
+
+      double x;
+      double y;
+   };
+
+   class Element
+   {
+   public:
+      Element() {}
+
+      void Set(int id0, int id1, int id2, int id3)
+      {
+         vertex_ids[0] = id0; vertex_ids[1] = id1;
+         vertex_ids[2] = id2; vertex_ids[3] = id3;
+      }
+
+      const int num_vertices = 4;
+      std::array<int,4> vertex_ids;
+   };
+
+   class BElement
+   {
+   public:
+      BElement() {}
+
+      void Set(int id0, int id1)
+      {
+         vertex_ids[0] = id0; vertex_ids[1] = id1;
+      }
+
+      const int num_vertices = 2;
+      std::array<int,2> vertex_ids;
+   };
+
+public:
+   DummyMesh() : num_vertices(11), num_elements(5), num_belements(9),
+      num_vparents(2)
+   {
+      // Vertices
+      V[6].Set(0.0, 2.0);  V[7].Set(1.0, 2.0);  V[8].Set(2.0, 2.0);
+      V[3].Set(0.0, 1.0);  V[4].Set(1.0, 1.0);  V[5].Set(2.0, 1.0);
+      V[0].Set(0.0, 0.0);  V[1].Set(1.0, 0.0);  V[2].Set(2.0, 0.0);
+
+      V[9].Set(0.0, 0.5);  V[10].Set(1.0, 0.5);
+
+      // Elements in this dummy format have their vertex indices listed
+      // lexographic order rather than going around the element as is normal
+      // in MFEM.
+      E[0].Set(0,1,9,10);
+      E[1].Set(9,10,3,4);
+      E[2].Set(1,2,4,5);
+      E[3].Set(3,4,6,7);
+      E[4].Set(4,5,7,8);
+
+      // Boundary Elements
+      // Bottom
+      B[0].Set(0,1);       B[1].Set(1,2);
+      // Top
+      B[2].Set(6,7);       B[3].Set(7,8);
+      // Left
+      B[4].Set(0,9);       B[5].Set(9,3);       B[6].Set(3,6);
+      // Right
+      B[7].Set(2,5);       B[8].Set(5,8);
+
+      // Set the vertex parents
+      VP[0] = std::make_tuple(9,0,3);
+      VP[1] = std::make_tuple(10,1,4);
+   }
+
+   const int num_vertices;
+   const int num_elements;
+   const int num_belements;
+   const int num_vparents;
+   Vertex V[11];                    // Mesh vertices
+   Element E[5];                    // Mesh elements
+   BElement B[9];                   // Mesh boundary elements
+   std::tuple<int,int,int> VP[2];   // Vertex parents (vid, par1id, par2id)
+};

miniapps/meshing/makefile

@@ -27,7 +27,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_MINIAPPS = mobius-strip klein-bottle toroid trimmer twist mesh-explorer\
 	shaper extruder mesh-optimizer minimal-surface polar-nc reflector
-PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting
+PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting ext-mesh-mapping
+
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -65,7 +66,7 @@ MESH_FILES = amr-quad-q2.mesh blade.mesh cube.mesh icf.mesh jagged.mesh\
  square01.mesh stretched2D.mesh
 $(MESH_FILES): %: $(SRC)%
 	ln -sf $(<) .
-mesh-optimizer pmesh-optimizer pmesh-fitting: | $(MESH_FILES)
+mesh-optimizer pmesh-optimizer pmesh-fitting ext-mesh-mapping: | $(MESH_FILES)
 .PHONY: copy-data
 copy-data: | $(MESH_FILES)
 endif
@@ -96,10 +97,13 @@ pminimal-surface-test-par: pminimal-surface
 reflector-test-seq: reflector
 	@$(call mfem-test-file,$<,, Meshing miniapp,reflected.mesh)
 
+
 # Testing: Specific execution options
 mesh-explorer-test-seq:
 	@true
-shaper-test-seq:
+shaper-test-seq:	
+	@true
+ext-mesh-mapping-test-par: 
 	@true
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
@@ -115,6 +119,7 @@ clean-build:
 	rm -f mesh-explorer shaper extruder trimmer reflector
 	rm -f mesh-optimizer pmesh-optimizer pmesh-fitting polar-nc
 	rm -f minimal-surface pminimal-surface
+	rm -f ext-mesh-mapping
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec: