Merge remote-tracking branch 'origin/master' into entity-sets-dev

# Conflicts:
#	mesh/pncmesh.hpp

examples/ex1_es.cpp

@@ -0,0 +1,266 @@
+//                                MFEM Example 1
+//
+// Compile with: make ex1
+//
+// Sample runs:  ex1 -m ../data/square-disc.mesh
+//               ex1 -m ../data/star.mesh
+//               ex1 -m ../data/escher.mesh
+//               ex1 -m ../data/fichera.mesh
+//               ex1 -m ../data/square-disc-p2.vtk -o 2
+//               ex1 -m ../data/square-disc-p3.mesh -o 3
+//               ex1 -m ../data/square-disc-nurbs.mesh -o -1
+//               ex1 -m ../data/disc-nurbs.mesh -o -1
+//               ex1 -m ../data/pipe-nurbs.mesh -o -1
+//               ex1 -m ../data/star-surf.mesh
+//               ex1 -m ../data/square-disc-surf.mesh
+//               ex1 -m ../data/inline-segment.mesh
+//               ex1 -m ../data/amr-quad.mesh
+//               ex1 -m ../data/amr-hex.mesh
+//               ex1 -m ../data/fichera-amr.mesh
+//               ex1 -m ../data/mobius-strip.mesh
+//               ex1 -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "./star-set.mesh";
+   int order = 1;
+   int rs = -1;
+   int ra = 0;
+   int bt = EntitySets::INVALID;
+   const char *bs = "Origin";
+   bool static_cond = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&rs, "-rs", "--refine-serial",
+                  "Number of serial refinement levels");
+   args.AddOption(&ra, "-ra", "--refine-adaptive",
+                  "Number of adaptive refinement levels");
+   args.AddOption(&bt, "-bt", "--bc-entity-type",
+                  "");
+   args.AddOption(&bs, "-bs", "--bc-entity-set-name",
+                  "");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    largest number that gives a final mesh with no more than 50,000
+   //    elements.
+   {
+      int ref_levels = ( rs >= 0 ) ? rs :
+                       (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+
+   if ( ra > 0 )
+   {
+      cout << "calling EnsureNCMesh" << endl;
+      mesh->EnsureNCMesh();
+      cout << "back from EnsureNCMesh" << endl;
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   cout << "Calling RandomRefinement " << ra << " times." << endl;
+   for (int l = 0; l < ra; l++)
+   {
+      mesh->RandomRefinement(0.2);
+   }
+   cout << "Done with refinement" << endl;
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   if ( mesh->ncmesh )
+   {
+      mesh->ncmesh->PrintStats(cout);
+
+      ofstream ofsV("vp.out");
+      ofstream ofsE("ce.out");
+      mesh->ncmesh->PrintVertexParents(ofsV);
+      mesh->ncmesh->PrintCoarseElements(ofsE);
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: "
+        << fespace->GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined 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;
+   if ( bt == EntitySets::INVALID )
+   {
+      if (mesh->bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(mesh->bdr_attributes.Max());
+         ess_bdr = 1;
+         fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+   }
+   else
+   {
+      fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                    ess_tdof_list);
+   }
+
+   cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm *b = new LinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 7. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 8. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 9. Assemble the bilinear form and the corresponding linear system,
+   //    applying any necessary transformations such as: eliminating boundary
+   //    conditions, applying conforming constraints for non-conforming AMR,
+   //    static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   SparseMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   cout << "Size of linear system: " << A.Height() << endl;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   //     solve the system A X = B with PCG.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+#else
+   // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(A);
+   umf_solver.Mult(B, X);
+#endif
+
+   // 11. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 12. Save the refined mesh and the solution. This output can be viewed later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+   ofstream mesh_ofs("refined.mesh");
+   mesh_ofs.precision(8);
+   mesh->Print(mesh_ofs);
+   ofstream sol_ofs("sol.gf");
+   sol_ofs.precision(8);
+   x.Save(sol_ofs);
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 14. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   if (order > 0) { delete fec; }
+   delete mesh;
+
+   return 0;
+}

examples/ex1p_es.cpp

@@ -0,0 +1,388 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+//               The following are examples of using EntitySets to define
+//               homogeneous Dirichlet boundary condition.  These examples
+//               require a modified mesh file and a specialized version of
+//               example 1 called "ex1p_es".
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 0 -bs Origin
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 1 -bs Axes
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh
+//                                    -bt 1 -bs "Negative Axes"
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh
+//                                    -bt 2 -bs "Interior Corner"
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh
+//                                    -bt 2 -bs "Exterior Corner"
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh
+//                                    -bt 3 -bs "Interior Corner"
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh
+//                                    -bt 3 -bs "Exterior Corner"
+//               mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 3 -bs "Steps"
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "./star-set.mesh";
+   int order = 1;
+   int rs = -1;
+   int rp = 2;
+   int ra = 0;
+   int bt = EntitySets::INVALID;
+   const char *bs = "Origin";
+   bool static_cond = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&rs, "-rs", "--refine-serial",
+                  "Number of serial refinement levels");
+   args.AddOption(&rp, "-rp", "--refine-parallel",
+                  "Number of parallel refinement levels");
+   args.AddOption(&ra, "-ra", "--refine-adaptive",
+                  "Number of adaptive refinement levels");
+   args.AddOption(&bt, "-bt", "--bc-entity-type",
+                  "");
+   args.AddOption(&bs, "-bs", "--bc-entity-set-name",
+                  "");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. 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 = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels = ( rs >= 0 ) ? rs :
+                       (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
+         mesh->UniformRefinement();
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+      if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     At this point we have a serial mesh containing an EntitySets
+     object which stores the current node/edge/face/element indices
+     for each entity in each set.  This data is duplicated on each MPI
+     rank.
+    */
+   if ( ra > 0 )
+   {
+      cout << "calling EnsureNCMesh" << endl;
+      mesh->EnsureNCMesh();
+      cout << "back from EnsureNCMesh" << endl;
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     We now have an NCEntitySets object which stores the node indices
+     describing each enity in each node/edge/face set and the element
+     indices for the elements in each element set.  This data is
+     duplicated on each MPI rank.
+    */
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   cout << "creating ParMesh from serial mesh" << endl;
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   cout << "done creating ParMesh from serial mesh" << endl;
+   delete mesh;
+   if ( pmesh->pent_sets )
+   {
+      cout << "pmesh->pent_sets is non NULL" << endl;
+      pmesh->pent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "pmesh->pent_sets is NULL" << endl;
+   }
+   /*
+     We now have a ParEntitySets object which marshals the data stored
+     in EntitySets objects.  The data has now been pruned so that each
+     rank only contains indices of local entities.
+
+     The NCEntitySets object remains unchanged...
+
+     If we have an NC mesh a different path is taken and the
+     EntitySets are ignored.
+
+     1) ParNCMesh is created from NCMesh
+       a) Creates a ParNCEntitySets object from ncmesh (every rank contains
+          information to find every entity)
+     2) ParNCMesh is pruned which involves renumbering elements and vertices
+     3) ParMesh is initialized from ParNCMesh
+     4) ParNCMesh::OnMeshUpdated is called
+     5) Mesh::GenerateNCFaceInfo is called
+    */
+   {
+      int par_ref_levels = rp;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
+         pmesh->UniformRefinement();
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+      if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
+   }
+   /*
+     RandomRefinement will end up calling
+     ParMesh::NonconformingRefinement which will create a new ParMesh
+     object using the ParNCMesh object and then call
+     ParMesh::OnMeshUpdated on this new mesh.
+    */
+
+   for (int l = 0; l < ra; l++)
+   {
+      pmesh->RandomRefinement(0.2);
+   }
+   if ( ra > 0 )
+   {
+      if ( pmesh->pent_sets )
+      {
+         cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
+         pmesh->pent_sets->PrintSetInfo(cout);
+      }
+      else
+      {
+         cout << "pmesh->pent_sets is NULL post random refinement" << endl;
+      }
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (pmesh->GetNodes())
+   {
+      fec = pmesh->GetNodes()->OwnFEC();
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    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;
+   if ( bt == EntitySets::INVALID )
+   {
+      if (pmesh->bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+         ess_bdr = 1;
+         fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+   }
+   else
+   {
+      fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                    ess_tdof_list);
+   }
+   for (int i=0; i<num_procs; i++)
+   {
+      if (myid == i)
+      {
+         cout << "Number of Dirichlet dofs on proc " << i << ": "
+              << ess_tdof_list.Size() << endl;
+      }
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm *b = new ParLinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   ParGridFunction x(fespace);
+   x = 0.0;
+
+   // 10. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+   }
+
+   // 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
+   //     preconditioner from hypre.
+   HypreSolver *amg = new HypreBoomerAMG(A);
+   HyprePCG *pcg = new HyprePCG(A);
+   pcg->SetTol(1e-12);
+   pcg->SetMaxIter(200);
+   pcg->SetPrintLevel(2);
+   pcg->SetPreconditioner(*amg);
+   pcg->Mult(B, X);
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 14. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 16. Free the used memory.
+   delete pcg;
+   delete amg;
+   delete a;
+   delete b;
+   delete fespace;
+   if (order > 0) { delete fec; }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/ex3p_es.cpp

@@ -0,0 +1,411 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p
+//
+// Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh
+//               mpirun -np 4 ex3p -m ../data/fichera.mesh
+//               mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
+//               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
+//
+// Description:  This example code solves a simple electromagnetic diffusion
+//               problem corresponding to the second order definite Maxwell
+//               equation curl curl E + E = f with boundary condition
+//               E x n = <given tangential field>. Here, we use a given exact
+//               solution E and compute the corresponding r.h.s. f.
+//               We discretize with Nedelec finite elements in 2D or 3D.
+//
+//               The example demonstrates the use of H(curl) finite element
+//               spaces with the curl-curl and the (vector finite element) mass
+//               bilinear form, as well as the computation of discretization
+//               error when the exact solution is known. Static condensation is
+//               also illustrated.
+//
+//               We recommend viewing examples 1-2 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+//void E_exact(const Vector &, Vector &);
+//void f_exact(const Vector &, Vector &);
+//double freq = 1.0, kappa;
+void f_const(const Vector &, Vector &);
+
+int dim;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   int rs = -1;
+   int rp = 2;
+   int ra = 0;
+   int bt = EntitySets::INVALID;
+   const char *bs = "Origin";
+   bool static_cond = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   /*
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   */
+   args.AddOption(&rs, "-rs", "--refine-serial",
+                  "Number of serial refinement levels");
+   args.AddOption(&rp, "-rp", "--refine-parallel",
+                  "Number of parallel refinement levels");
+   args.AddOption(&ra, "-ra", "--refine-adaptive",
+                  "Number of adaptive refinement levels");
+   args.AddOption(&bt, "-bt", "--bc-entity-type",
+                  "");
+   args.AddOption(&bs, "-bs", "--bc-entity-set-name",
+                  "");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   // kappa = freq * M_PI;
+
+   // 3. 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 = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels = ( rs >= 0 ) ? rs :
+                       (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
+         mesh->UniformRefinement();
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+      if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     At this point we have a serial mesh containing an EntitySets
+     object which stores the current node/edge/face/element indices
+     for each entity in each set.  This data is duplicated on each MPI
+     rank.
+    */
+   if ( ra > 0 )
+   {
+      cout << "calling EnsureNCMesh" << endl;
+      mesh->EnsureNCMesh();
+      cout << "back from EnsureNCMesh" << endl;
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     We now have an NCEntitySets object which stores the node indices
+     describing each enity in each node/edge/face set and the element
+     indices for the elements in each element set.  This data is
+     duplicated on each MPI rank.
+    */
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   if ( pmesh->pent_sets )
+   {
+      cout << "pmesh->pent_sets is non NULL" << endl;
+      pmesh->pent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "pmesh->pent_sets is NULL" << endl;
+   }
+   {
+      int par_ref_levels = rp;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
+         pmesh->UniformRefinement();
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+      if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
+   }
+   pmesh->ReorientTetMesh();
+   pmesh->ent_sets->PrintSetInfo(cout);
+
+   for (int l = 0; l < ra; l++)
+   {
+      pmesh->RandomRefinement(0.2);
+   }
+   if ( ra > 0 )
+   {
+      if ( pmesh->pent_sets )
+      {
+         cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
+         pmesh->pent_sets->PrintSetInfo(cout);
+      }
+      else
+      {
+         cout << "pmesh->pent_sets is NULL post random refinement" << endl;
+      }
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    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;
+   if ( bt == EntitySets::INVALID )
+   {
+      if (pmesh->bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+         ess_bdr = 1;
+         fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+   }
+   else
+   {
+      fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                    ess_tdof_list);
+   }
+   if (myid == 0)
+   {
+      cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
+   }
+
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_const);
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x by projecting the exact
+   //    solution. Note that only values from the boundary edges will be used
+   //    when eliminating the non-homogeneous boundary condition to modify the
+   //    r.h.s. vector b.
+   ParGridFunction x(fespace);
+   // VectorFunctionCoefficient E(sdim, E_exact);
+   // x.ProjectCoefficient(E);
+   x = 0.0;
+
+   // 10. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *sigma = new ConstantCoefficient(1.0);
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+   }
+
+   // 12. Define and apply a parallel PCG solver for AX=B with the AMS
+   //     preconditioner from hypre.
+   ParFiniteElementSpace *prec_fespace =
+      (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
+   HypreSolver *ams = new HypreAMS(A, prec_fespace);
+   HyprePCG *pcg = new HyprePCG(A);
+   pcg->SetTol(1e-12);
+   pcg->SetMaxIter(500);
+   pcg->SetPrintLevel(2);
+   pcg->SetPreconditioner(*ams);
+   pcg->Mult(B, X);
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   /*
+   // 14. Compute and print the L^2 norm of the error.
+   {
+      double err = x.ComputeL2Error(E);
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+   */
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete pcg;
+   delete ams;
+   delete a;
+   delete sigma;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+/*
+void E_exact(const Vector &x, Vector &E)
+{
+   if (dim == 3)
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(2));
+      E(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { E(2) = 0.0; }
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
+      f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}
+*/
+void f_const(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = 1.0;
+      f(1) = 1.0;
+      f(2) = 1.0;
+   }
+   else
+   {
+      f(0) = 1.0;
+      f(1) = 1.0;
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

examples/ex4p_es.cpp

@@ -0,0 +1,438 @@
+//                       MFEM Example 4 - Parallel Version
+//
+// Compile with: make ex4p
+//
+// Sample runs:  mpirun -np 4 ex4p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex4p -m ../data/star.mesh
+//               mpirun -np 4 ex4p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex4p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex4p -m ../data/escher.mesh -o 2 -sc
+//               mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex4p -m ../data/fichera-q3.mesh -o 2 -sc
+//               mpirun -np 4 ex4p -m ../data/square-disc-nurbs.mesh -o 3
+//               mpirun -np 4 ex4p -m ../data/beam-hex-nurbs.mesh -o 3
+//               mpirun -np 4 ex4p -m ../data/periodic-square.mesh -no-bc
+//               mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
+//               mpirun -np 4 ex4p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
+//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
+//
+// Description:  This example code solves a simple 2D/3D H(div) diffusion
+//               problem corresponding to the second order definite equation
+//               -grad(alpha div F) + beta F = f with boundary condition F dot n
+//               = <given normal field>. Here, we use a given exact solution F
+//               and compute the corresponding r.h.s. f.  We discretize with
+//               Raviart-Thomas finite elements.
+//
+//               The example demonstrates the use of H(div) finite element
+//               spaces with the grad-div and H(div) vector finite element mass
+//               bilinear form, as well as the computation of discretization
+//               error when the exact solution is known. Bilinear form
+//               hybridization and static condensation are also illustrated.
+//
+//               We recommend viewing examples 1-3 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, F, and r.h.s., f. See below for implementation.
+//void F_exact(const Vector &, Vector &);
+//void f_exact(const Vector &, Vector &);
+//double freq = 1.0, kappa;
+void f_const(const Vector &, Vector &);
+
+int dim;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int rs = -1;
+   int rp = 2;
+   int ra = 0;
+   int bt = EntitySets::INVALID;
+   const char *bs = "Origin";
+   bool set_bc = true;
+   bool static_cond = false;
+   bool hybridization = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&set_bc, "-bc", "--impose-bc", "-no-bc", "--dont-impose-bc",
+                  "Impose or not essential boundary conditions.");
+   args.AddOption(&rs, "-rs", "--refine-serial",
+                  "Number of serial refinement levels");
+   args.AddOption(&rp, "-rp", "--refine-parallel",
+                  "Number of parallel refinement levels");
+   args.AddOption(&ra, "-ra", "--refine-adaptive",
+                  "Number of adaptive refinement levels");
+   args.AddOption(&bt, "-bt", "--bc-entity-type",
+                  "");
+   args.AddOption(&bs, "-bs", "--bc-entity-set-name",
+                  "");
+   // args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+   //                 " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&hybridization, "-hb", "--hybridization", "-no-hb",
+                  "--no-hybridization", "Enable hybridization.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   // kappa = freq * M_PI;
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume, as well as periodic meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels = ( rs >= 0 ) ? rs :
+                       (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
+         mesh->UniformRefinement();
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+      if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     At this point we have a serial mesh containing an EntitySets
+     object which stores the current node/edge/face/element indices
+     for each entity in each set.  This data is duplicated on each MPI
+     rank.
+    */
+   if ( ra > 0 )
+   {
+      cout << "calling EnsureNCMesh" << endl;
+      mesh->EnsureNCMesh();
+      cout << "back from EnsureNCMesh" << endl;
+   }
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+   /*
+     We now have an NCEntitySets object which stores the node indices
+     describing each enity in each node/edge/face set and the element
+     indices for the elements in each element set.  This data is
+     duplicated on each MPI rank.
+    */
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them (this is needed in the ADS solver below).
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   if ( pmesh->pent_sets )
+   {
+      cout << "pmesh->pent_sets is non NULL" << endl;
+      pmesh->pent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "pmesh->pent_sets is NULL" << endl;
+   }
+   {
+      int par_ref_levels = rp;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   for (int l = 0; l < ra; l++)
+   {
+      pmesh->RandomRefinement(0.2);
+   }
+   if ( ra > 0 )
+   {
+      if ( pmesh->pent_sets )
+      {
+         cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
+         pmesh->pent_sets->PrintSetInfo(cout);
+      }
+      else
+      {
+         cout << "pmesh->pent_sets is NULL post random refinement" << endl;
+      }
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Raviart-Thomas finite elements of the specified order.
+   FiniteElementCollection *fec = new RT_FECollection(order-1, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    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;
+   if ( bt == EntitySets::INVALID )
+   {
+      if (pmesh->bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+         ess_bdr = set_bc ? 1 : 0;
+         fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+   }
+   else
+   {
+      fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                    ess_tdof_list);
+   }
+   if (myid == 0)
+   {
+      cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_const);
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x by projecting the exact
+   //    solution. Note that only values from the boundary faces will be used
+   //    when eliminating the non-homogeneous boundary condition to modify the
+   //    r.h.s. vector b.
+   ParGridFunction x(fespace);
+   // VectorFunctionCoefficient F(sdim, F_exact);
+   // x.ProjectCoefficient(F);
+   x = 0.0;
+
+   // 10. Set up the parallel bilinear form corresponding to the H(div)
+   //     diffusion operator grad alpha div + beta I, by adding the div-div and
+   //     the mass domain integrators.
+   Coefficient *alpha = new ConstantCoefficient(1.0);
+   Coefficient *beta  = new ConstantCoefficient(1.0);
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation,
+   //     hybridization, etc.
+   FiniteElementCollection *hfec = NULL;
+   ParFiniteElementSpace *hfes = NULL;
+   if (static_cond)
+   {
+      a->EnableStaticCondensation();
+   }
+   else if (hybridization)
+   {
+      hfec = new DG_Interface_FECollection(order-1, dim);
+      hfes = new ParFiniteElementSpace(pmesh, hfec);
+      a->EnableHybridization(hfes, new NormalTraceJumpIntegrator(),
+                             ess_tdof_list);
+   }
+   a->Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   HYPRE_Int glob_size = A.GetGlobalNumRows();
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << glob_size << endl;
+   }
+
+   // 12. Define and apply a parallel PCG solver for A X = B with the 2D AMS or
+   //     the 3D ADS preconditioners from hypre. If using hybridization, the
+   //     system is preconditioned with hypre's BoomerAMG.
+   HypreSolver *prec = NULL;
+   CGSolver *pcg = new CGSolver(A.GetComm());
+   pcg->SetOperator(A);
+   pcg->SetRelTol(1e-12);
+   pcg->SetMaxIter(500);
+   pcg->SetPrintLevel(1);
+   if (hybridization) { prec = new HypreBoomerAMG(A); }
+   else
+   {
+      ParFiniteElementSpace *prec_fespace =
+         (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
+      if (dim == 2)   { prec = new HypreAMS(A, prec_fespace); }
+      else            { prec = new HypreADS(A, prec_fespace); }
+   }
+   pcg->SetPreconditioner(*prec);
+   pcg->Mult(B, X);
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+   /*
+   // 14. Compute and print the L^2 norm of the error.
+   {
+      double err = x.ComputeL2Error(F);
+      if (myid == 0)
+      {
+         cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+   */
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete pcg;
+   delete prec;
+   delete hfes;
+   delete hfec;
+   delete a;
+   delete alpha;
+   delete beta;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+/*
+// The exact solution (for non-surface meshes)
+void F_exact(const Vector &p, Vector &F)
+{
+   int dim = p.Size();
+
+   double x = p(0);
+   double y = p(1);
+   // double z = (dim == 3) ? p(2) : 0.0;
+
+   F(0) = cos(kappa*x)*sin(kappa*y);
+   F(1) = cos(kappa*y)*sin(kappa*x);
+   if (dim == 3)
+   {
+      F(2) = 0.0;
+   }
+}
+
+// The right hand side
+void f_exact(const Vector &p, Vector &f)
+{
+   int dim = p.Size();
+
+   double x = p(0);
+   double y = p(1);
+   // double z = (dim == 3) ? p(2) : 0.0;
+
+   double temp = 1 + 2*kappa*kappa;
+
+   f(0) = temp*cos(kappa*x)*sin(kappa*y);
+   f(1) = temp*cos(kappa*y)*sin(kappa*x);
+   if (dim == 3)
+   {
+      f(2) = 0;
+   }
+}
+*/
+void f_const(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = 1.0;
+      f(1) = 1.0;
+      f(2) = 1.0;
+   }
+   else
+   {
+      f(0) = 1.0;
+      f(1) = 1.0;
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

examples/ex6p_es.cpp

@@ -0,0 +1,325 @@
+//                       MFEM Example 6 - Parallel Version
+//
+// Compile with: make ex6p
+//
+// Sample runs:  mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 1
+//               mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/square-disc-nurbs.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/star.mesh -o 3
+//               mpirun -np 4 ex6p -m ../data/escher.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/fichera.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/disc-nurbs.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/ball-nurbs.mesh
+//               mpirun -np 4 ex6p -m ../data/pipe-nurbs.mesh
+//               mpirun -np 4 ex6p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/square-disc-surf.mesh -o 2
+//               mpirun -np 4 ex6p -m ../data/amr-quad.mesh
+//
+// Description:  This is a version of Example 1 with a simple adaptive mesh
+//               refinement loop. The problem being solved is again the Laplace
+//               equation -Delta u = 1 with homogeneous Dirichlet boundary
+//               conditions. The problem is solved on a sequence of meshes which
+//               are locally refined in a conforming (triangles, tetrahedrons)
+//               or non-conforming (quadrilaterals, hexahedra) manner according
+//               to a simple ZZ error estimator.
+//
+//               The example demonstrates MFEM's capability to work with both
+//               conforming and nonconforming refinements, in 2D and 3D, on
+//               linear, curved and surface meshes. Interpolation of functions
+//               from coarse to fine meshes, as well as persistent GLVis
+//               visualization are also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static int max_dofs = 100000;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "./star-set.mesh";
+   int order = 1;
+   int bt = EntitySets::INVALID;
+   const char *bs = "";
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&max_dofs, "-md", "--max-dofs",
+                  "Maximum number of degrees of freedom.");
+   args.AddOption(&bt, "-bt", "--bc-entity-type",
+                  "");
+   args.AddOption(&bs, "-bs", "--bc-entity-set-name",
+                  "");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. 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 = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution.
+   //    Also project a NURBS mesh to a piecewise-quadratic curved mesh. Make
+   //    sure that the mesh is non-conforming.
+   if (mesh->NURBSext)
+   {
+      mesh->UniformRefinement();
+      mesh->SetCurvature(2);
+   }
+   mesh->EnsureNCMesh();
+   if ( mesh->ent_sets )
+   {
+      cout << "mesh->ent_sets is non NULL" << endl;
+      mesh->ent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "mesh->ent_sets is NULL" << endl;
+   }
+
+   // 5. Define a parallel mesh by partitioning the serial mesh.
+   //    Once the parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   if ( pmesh.pent_sets )
+   {
+      cout << "pmesh->pent_sets is non NULL" << endl;
+      pmesh.pent_sets->PrintSetInfo(cout);
+   }
+   else
+   {
+      cout << "pmesh->pent_sets is NULL" << endl;
+   }
+
+   // 6. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+   Array<int> ess_tdof_list;
+   if ( bt == EntitySets::INVALID )
+   {
+      if (pmesh.bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+         ess_bdr = 1;
+         fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+   }
+   else
+   {
+      fespace.GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                   ess_tdof_list);
+   }
+
+   // 7. As in Example 1p, we set up bilinear and linear forms corresponding to
+   //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
+   //    problem yet, this will be done in the main loop.
+   ParBilinearForm a(&fespace);
+   ParLinearForm b(&fespace);
+
+   ConstantCoefficient one(1.0);
+
+   BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
+   a.AddDomainIntegrator(integ);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+
+   // 8. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   ParGridFunction x(&fespace);
+   x = 0;
+
+   // 9. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream sout;
+   if (visualization)
+   {
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (myid == 0)
+         {
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+            cout << "GLVis visualization disabled.\n";
+         }
+         visualization = false;
+      }
+
+      sout.precision(8);
+   }
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     with L2 projection in the smoothing step to better handle hanging
+   //     nodes and parallel partitioning. We need to supply a space for the
+   //     discontinuous flux (L2) and a space for the smoothed flux (H(div) is
+   //     used here).
+   L2_FECollection flux_fec(order, dim);
+   ParFiniteElementSpace flux_fes(&pmesh, &flux_fec, sdim);
+   RT_FECollection smooth_flux_fec(order-1, dim);
+   ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
+   // Another possible option for the smoothed flux space:
+   // H1_FECollection smooth_flux_fec(order, dim);
+   // ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec, dim);
+   L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fes, smooth_flux_fes);
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.7);
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   // const int max_dofs = 100000;
+   for (int it = 0; ; it++)
+   {
+      HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
+      if (myid == 0)
+      {
+         cout << "\nAMR iteration " << it << endl;
+         cout << "Number of unknowns: " << global_dofs << endl;
+      }
+
+      // 13. Assemble the stiffness matrix and the right-hand side. Note that
+      //     MFEM doesn't care at this point that the mesh is nonconforming
+      //     and parallel. The FE space is considered 'cut' along hanging
+      //     edges/faces, and also across processor boundaries.
+      a.Assemble();
+      b.Assemble();
+
+      // 14. Create the parallel linear system: eliminate boundary conditions,
+      //     constrain hanging nodes and nodes across processor boundaries.
+      //     The system will be solved for true (unconstrained/unique) DOFs only.
+      // Array<int> ess_tdof_list;
+      if ( bt == EntitySets::INVALID )
+      {
+         if (pmesh.bdr_attributes.Size())
+         {
+            Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+            ess_bdr = 1;
+            fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+         }
+      }
+      else
+      {
+         fespace.GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
+                                      ess_tdof_list);
+      }
+
+      HypreParMatrix A;
+      Vector B, X;
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 15. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
+      //     preconditioner from hypre.
+      HypreBoomerAMG amg;
+      amg.SetPrintLevel(0);
+      CGSolver pcg(A.GetComm());
+      pcg.SetPreconditioner(amg);
+      pcg.SetOperator(A);
+      pcg.SetRelTol(1e-6);
+      pcg.SetMaxIter(200);
+      pcg.SetPrintLevel(3); // print the first and the last iterations only
+      pcg.Mult(B, X);
+
+      // 16. Extract the parallel grid function corresponding to the finite element
+      //     approximation X. This is the local solution on each processor.
+      a.RecoverFEMSolution(X, b, x);
+
+      // 17. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout << "solution\n" << pmesh << x << flush;
+      }
+
+      if (global_dofs > max_dofs)
+      {
+         if (myid == 0)
+         {
+            cout << "Reached the maximum number of dofs. Stop." << endl;
+         }
+         break;
+      }
+
+      // 18. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+      refiner.Apply(pmesh);
+      if (refiner.Stop())
+      {
+         if (myid == 0)
+         {
+            cout << "Stopping criterion satisfied. Stop." << endl;
+         }
+         break;
+      }
+
+      // 19. Update the finite element space (recalculate the number of DOFs,
+      //     etc.) and create a grid function update matrix. Apply the matrix
+      //     to any GridFunctions over the space. In this case, the update
+      //     matrix is an interpolation matrix so the updated GridFunction will
+      //     still represent the same function as before refinement.
+      fespace.Update();
+      x.Update();
+
+      // 20. Load balance the mesh, and update the space and solution. Currently
+      //     available only for nonconforming meshes.
+      if (pmesh.Nonconforming())
+      {
+         pmesh.Rebalance();
+
+         // Update the space and the GridFunction. This time the update matrix
+         // redistributes the GridFunction among the processors.
+         fespace.Update();
+         x.Update();
+      }
+
+      // 21. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+
+   MPI_Finalize();
+   return 0;
+}

examples/fichera-mixed-set.mesh

@@ -0,0 +1,162 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+3
+
+elements
+14
+1 4 13 15 21 25
+1 4 12 13 15 21
+1 4 13 21 22 25
+1 4 15 24 21 25
+1 4 13 15 25 16
+1 5  0  1  4  3  9 10 13 12
+1 5  8  9 12 11 17 18 21 20
+1 5  2  3  6  5 11 12 15 14
+1 6  3  4  6 12 13 15
+1 6  4  7  6 13 16 15
+1 6 12 13 21  9 10 18
+1 6 13 22 21 10 19 18
+1 6 11 14 20 12 15 21
+1 6 15 21 24 14 20 23
+
+boundary
+30
+ 1 3  5  6  3  2
+ 2 2  3  6  4
+ 2 2  4  6  7
+ 3 3  3  4  1  0
+ 4 3 11 12  9  8
+ 5 3  2  3 12 11
+ 6 3  0  1 10  9
+ 7 2  9 10 18
+ 7 2 10 19 18
+ 8 3  8  9 18 17
+ 9 3  1  4 13 10
+10 3  4  7 16 13
+11 2 13 16 25
+11 2 13 25 22
+12 3 10 13 22 19
+13 3  7  6 15 16
+14 3  6  5 14 15
+15 3 15 14 23 24
+16 2 16 15 25
+16 2 15 24 25
+17 3  5  2 11 14
+18 3  3  0  9 12
+19 3 11  8 17 20
+20 2 11 20 14
+20 2 14 20 23
+21 3 17 18 21 20
+22 3 18 19 22 21
+23 2 21 22 25
+23 2 21 25 24
+24 3 20 21 24 23
+
+vertices
+26
+3
+0 -1 -1
+1 -1 -1
+-1 0 -1
+0 0 -1
+1 0 -1
+-1 1 -1
+0 1 -1
+1 1 -1
+-1 -1 0
+0 -1 0
+1 -1 0
+-1 0 0
+0 0 0
+1 0 0
+-1 1 0
+0 1 0
+1 1 0
+-1 -1 1
+0 -1 1
+1 -1 1
+-1 0 1
+0 0 1
+1 0 1
+-1 1 1
+0 1 1
+1 1 1
+
+MFEM sets v1.0
+
+vertex_sets
+1
+
+Origin
+1
+12
+
+edge_sets
+2
+
+Axes
+3
+12 13
+12 15
+12 21
+
+Negative Axes
+3
+12 9
+12 11
+12 3
+
+face_sets
+2
+
+Interior Corner
+3
+3 11 12  9 8
+3  2  3 12 11
+3  3  0  9 12
+
+Exterior Corner
+15
+2 13 16 25
+2 13 25 22
+2 16 15 25
+2 15 24 25
+2 21 22 25
+2 21 25 24
+3 10 13 22 19
+3 4 7 16 13
+3 1 4 13 10
+3 7 6 15 16
+3 6 5 14 15
+3 15 14 23 24
+3 20 21 24 23
+3 18 19 22 21
+3 17 18 21 20
+
+element_sets
+3
+
+Interior Corner
+3
+5 6 7
+
+Exterior Corner
+5
+0 1 2 3 4
+
+Steps
+3
+6 8 9

examples/fichera-set.mesh

@@ -0,0 +1,145 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+7
+1 5 0 1 4 3 9 10 13 12
+1 5 3 4 7 6 12 13 16 15
+1 5 2 3 6 5 11 12 15 14
+1 5 8 9 12 11 17 18 21 20
+1 5 9 10 13 12 18 19 22 21
+1 5 12 13 16 15 21 22 25 24
+1 5 11 12 15 14 20 21 24 23
+
+boundary
+24
+1 3 5 6 3 2
+2 3 6 7 4 3
+3 3 3 4 1 0
+4 3 11 12 9 8
+5 3 2 3 12 11
+6 3 0 1 10 9
+7 3 9 10 19 18
+8 3 8 9 18 17
+9 3 1 4 13 10
+10 3 4 7 16 13
+11 3 13 16 25 22
+12 3 10 13 22 19
+13 3 7 6 15 16
+14 3 6 5 14 15
+15 3 15 14 23 24
+16 3 16 15 24 25
+17 3 5 2 11 14
+18 3 3 0 9 12
+19 3 11 8 17 20
+20 3 14 11 20 23
+21 3 17 18 21 20
+22 3 18 19 22 21
+23 3 21 22 25 24
+24 3 20 21 24 23
+
+vertices
+26
+3
+0 -1 -1
+1 -1 -1
+-1 0 -1
+0 0 -1
+1 0 -1
+-1 1 -1
+0 1 -1
+1 1 -1
+-1 -1 0
+0 -1 0
+1 -1 0
+-1 0 0
+0 0 0
+1 0 0
+-1 1 0
+0 1 0
+1 1 0
+-1 -1 1
+0 -1 1
+1 -1 1
+-1 0 1
+0 0 1
+1 0 1
+-1 1 1
+0 1 1
+1 1 1
+
+MFEM sets v1.0
+
+vertex_sets
+1
+
+Origin
+1
+12
+
+edge_sets
+2
+
+Axes
+3
+12 13
+12 15
+12 21
+
+Negative Axes
+3
+12 9
+12 11
+12 3
+
+face_sets
+2
+
+Interior Corner
+3
+3 11 12  9 8
+3  2  3 12 11
+3  3  0  9 12
+
+Exterior Corner
+12
+3 13 16 25 22
+3 16 15 24 25
+3 21 22 25 24
+3 10 13 22 19
+3 4 7 16 13
+3 1 4 13 10
+3 7 6 15 16
+3 6 5 14 15
+3 15 14 23 24
+3 20 21 24 23
+3 18 19 22 21
+3 17 18 21 20
+
+element_sets
+3
+
+Interior Corner
+3
+0 2 3
+
+Exterior Corner
+1
+5
+
+Steps
+2
+1 3

examples/star-mixed-set.mesh

@@ -0,0 +1,158 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+30
+1 3 0 11 26 14
+1 3 0 14 27 17
+1 3 0 17 28 20
+1 3 0 20 29 23
+1 3 0 23 30 11
+1 2 11 1 26
+1 2 1 12 26
+1 3 26 12 3 13
+1 2 26 13 2
+1 2 14 26 2
+1 2 14 2 27
+1 2 2 15 27
+1 3 27 15 5 16
+1 2 27 16 4
+1 2 17 27 4
+1 2 17 4 28
+1 2 4 18 28
+1 3 28 18 7 19
+1 2 28 19 6
+1 2 20 28 6
+1 2 20 6 29
+1 2 6 21 29
+1 3 29 21 9 22
+1 2 29 22 8
+1 2 23 29 8
+1 2 23 8 30
+1 2 8 24 30
+1 3 30 24 10 25
+1 2 30 25 1
+1 2 11 30 1
+
+boundary
+20
+1 1 13 2
+1 1 12 3
+1 1 16 4
+1 1 15 5
+1 1 19 6
+1 1 18 7
+1 1 22 8
+1 1 21 9
+1 1 25 1
+1 1 24 10
+1 1 3 13
+1 1 1 12
+1 1 5 16
+1 1 2 15
+1 1 7 19
+1 1 4 18
+1 1 9 22
+1 1 6 21
+1 1 10 25
+1 1 8 24
+
+vertices
+31
+2
+0 0
+1 0
+0.309017 0.951057
+1.30902 0.951057
+-0.809017 0.587785
+-0.5 1.53884
+-0.809017 -0.587785
+-1.61803 0
+0.309017 -0.951057
+-0.5 -1.53884
+1.30902 -0.951057
+0.5 0
+1.15451 0.475529
+0.809019 0.951057
+0.154508 0.475529
+-0.0954915 1.24495
+-0.654508 1.06331
+-0.404508 0.293893
+-1.21352 0.293893
+-1.21352 -0.293892
+-0.404508 -0.293893
+-0.654508 -1.06331
+-0.0954915 -1.24495
+0.154508 -0.475529
+0.809019 -0.951057
+1.15451 -0.475529
+0.654509 0.475529
+-0.25 0.769421
+-0.809016 0
+-0.25 -0.76942
+0.654509 -0.475529
+
+MFEM sets v1.0
+
+vertex_sets
+3
+
+Origin
+1
+0
+
+Tent
+5
+1 2 4 6 8
+
+Gazebo
+5
+3 5 7 9 10
+
+edge_sets
+2
+
+Columbine
+5
+1 11
+2 14
+4 17
+6 20
+8 23
+
+Lily
+5
+0 11
+0 14
+0 17
+0 20
+0 23
+
+element_sets
+3
+
+Flying Squirrel
+3
+7 17 27
+
+Sea Lion
+4
+12 17 22 27
+
+Pinwheel
+5
+8 13 18 23 28

examples/star-set.mesh

@@ -0,0 +1,143 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+20
+1 3 0 11 26 14
+1 3 0 14 27 17
+1 3 0 17 28 20
+1 3 0 20 29 23
+1 3 0 23 30 11
+1 3 11 1 12 26
+1 3 26 12 3 13
+1 3 14 26 13 2
+1 3 14 2 15 27
+1 3 27 15 5 16
+1 3 17 27 16 4
+1 3 17 4 18 28
+1 3 28 18 7 19
+1 3 20 28 19 6
+1 3 20 6 21 29
+1 3 29 21 9 22
+1 3 23 29 22 8
+1 3 23 8 24 30
+1 3 30 24 10 25
+1 3 11 30 25 1
+
+boundary
+20
+1 1 13 2
+1 1 12 3
+1 1 16 4
+1 1 15 5
+1 1 19 6
+1 1 18 7
+1 1 22 8
+1 1 21 9
+1 1 25 1
+1 1 24 10
+1 1 3 13
+1 1 1 12
+1 1 5 16
+1 1 2 15
+1 1 7 19
+1 1 4 18
+1 1 9 22
+1 1 6 21
+1 1 10 25
+1 1 8 24
+
+vertices
+31
+2
+0 0
+1 0
+0.309017 0.951057
+1.30902 0.951057
+-0.809017 0.587785
+-0.5 1.53884
+-0.809017 -0.587785
+-1.61803 0
+0.309017 -0.951057
+-0.5 -1.53884
+1.30902 -0.951057
+0.5 0
+1.15451 0.475529
+0.809019 0.951057
+0.154508 0.475529
+-0.0954915 1.24495
+-0.654508 1.06331
+-0.404508 0.293893
+-1.21352 0.293893
+-1.21352 -0.293892
+-0.404508 -0.293893
+-0.654508 -1.06331
+-0.0954915 -1.24495
+0.154508 -0.475529
+0.809019 -0.951057
+1.15451 -0.475529
+0.654509 0.475529
+-0.25 0.769421
+-0.809016 0
+-0.25 -0.76942
+0.654509 -0.475529
+
+MFEM sets v1.0
+
+vertex_sets
+3
+
+Origin
+1
+0
+
+Tent
+5
+1 2 4 6 8
+
+Gazebo
+5
+3 5 7 9 10
+
+edge_sets
+2
+
+Columbine
+5
+1 11
+2 14
+4 17
+6 20
+8 23
+
+Lily
+5
+0 11
+0 14
+0 17
+0 20
+0 23
+
+element_sets
+2
+
+Flying Squirrel
+3
+6 12 18
+
+Sea Lion
+4
+9 12 15 18

fem/fespace.cpp

@@ -561,6 +561,155 @@ void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
    }
 }
 
+void FiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
+                                           int set_index,
+                                           Array<int> &ess_vdofs,
+                                           int component) const
+{
+   Array<int> vdofs, dofs;
+
+   ess_vdofs.SetSize(GetVSize());
+   ess_vdofs = 0;
+
+   MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
+               "entity set information");
+   if (!mesh->ent_sets->SetExists(type, set_index))
+   {
+      ostringstream oss; oss << "Entity set of type \""
+                             << EntitySets::GetTypeName(type)
+                             << "\" and index " << set_index
+                             << " was not found.";
+
+      MFEM_VERIFY(false, oss.str().c_str());
+   }
+
+   set<int>::iterator it;
+   for (it=(*mesh->ent_sets)(type, set_index).begin();
+        it!=(*mesh->ent_sets)(type, set_index).end(); it++)
+   {
+      int ent_index = *it;
+      cout << "collecting vdofs for entity " << ent_index << "->";
+      if (component < 0)
+      {
+         switch (type)
+         {
+            case EntitySets::VERTEX:
+               GetVertexVDofs(ent_index, vdofs);
+               break;
+            case EntitySets::EDGE:
+               GetEdgeVDofs(ent_index, vdofs);
+               break;
+            case EntitySets::FACE:
+               GetFaceVDofs(ent_index, vdofs);
+               break;
+            case EntitySets::ELEMENT:
+               GetElementVDofs(ent_index, vdofs);
+               break;
+            default:
+               mfem_error("GetEssentialVDofs: Invalid entity type");
+         }
+         vdofs.Print(cout);
+         mark_dofs(vdofs, ess_vdofs);
+      }
+      else
+      {
+         switch (type)
+         {
+            case EntitySets::VERTEX:
+               GetVertexDofs(ent_index, dofs);
+               break;
+            case EntitySets::EDGE:
+               GetEdgeDofs(ent_index, dofs);
+               break;
+            case EntitySets::FACE:
+               GetFaceDofs(ent_index, dofs);
+               break;
+            case EntitySets::ELEMENT:
+               GetElementDofs(ent_index, dofs);
+               break;
+            default:
+               mfem_error("GetEssentialDofs: Invalid entity type");
+         }
+         for (int d = 0; d < dofs.Size(); d++)
+         { dofs[d] = DofToVDof(dofs[d], component); }
+         mark_dofs(dofs, ess_vdofs);
+      }
+   }
+
+   if (mesh->ncmesh)
+   {
+      Array<int> es_verts, es_edges, es_faces;
+      mesh->ncmesh->GetEntitySetClosure(type, set_index,
+                                        es_verts, es_edges, es_faces);
+      cout << "returned from get closure" << endl;
+      for (int i = 0; i < es_verts.Size(); i++)
+      {
+         if (es_verts[i] < GetNV())
+         {
+            if (component < 0)
+            {
+               GetVertexVDofs(es_verts[i], vdofs);
+               mark_dofs(vdofs, ess_vdofs);
+            }
+            else
+            {
+               GetVertexDofs(es_verts[i], dofs);
+               for (int d = 0; d < dofs.Size(); d++)
+               { dofs[d] = DofToVDof(dofs[d], component); }
+               mark_dofs(dofs, ess_vdofs);
+            }
+         }
+      }
+      for (int i = 0; i < es_edges.Size(); i++)
+      {
+         if (es_edges[i] < GetMesh()->GetNEdges())
+         {
+            if (component < 0)
+            {
+               GetEdgeVDofs(es_edges[i], vdofs);
+               mark_dofs(vdofs, ess_vdofs);
+            }
+            else
+            {
+               GetEdgeDofs(es_edges[i], dofs);
+               for (int d = 0; d < dofs.Size(); d++)
+               { dofs[d] = DofToVDof(dofs[d], component); }
+               mark_dofs(dofs, ess_vdofs);
+            }
+         }
+      }
+      for (int i = 0; i < es_faces.Size(); i++)
+      {
+         if (es_faces[i] < GetMesh()->GetNFaces())
+         {
+            if (component < 0)
+            {
+               GetFaceVDofs(es_faces[i], vdofs);
+               mark_dofs(vdofs, ess_vdofs);
+            }
+            else
+            {
+               GetFaceDofs(es_faces[i], dofs);
+               for (int d = 0; d < dofs.Size(); d++)
+               { dofs[d] = DofToVDof(dofs[d], component); }
+               mark_dofs(dofs, ess_vdofs);
+            }
+         }
+      }
+   }
+}
+
+void FiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
+                                           const string & set_name,
+                                           Array<int> &ess_vdofs,
+                                           int component) const
+{
+   MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
+               "entity set information");
+   GetEssentialVDofs(type, mesh->ent_sets->GetSetIndex(type, set_name),
+                     ess_vdofs, component);
+}
+
 void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
                                               Array<int> &ess_tdof_list,
                                               int component)
@@ -579,6 +728,36 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
    MarkerToList(ess_tdofs, ess_tdof_list);
 }
 
+void FiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
+                                              int set_index,
+                                              Array<int> &ess_tdof_list,
+                                              int component)
+{
+   Array<int> ess_vdofs, ess_tdofs;
+   GetEssentialVDofs(type, set_index, ess_vdofs, component);
+   const SparseMatrix *R = GetConformingRestriction();
+   if (!R)
+   {
+      ess_tdofs.MakeRef(ess_vdofs);
+   }
+   else
+   {
+      R->BooleanMult(ess_vdofs, ess_tdofs);
+   }
+   MarkerToList(ess_tdofs, ess_tdof_list);
+}
+
+void FiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
+                                              const string & set_name,
+                                              Array<int> &ess_tdof_list,
+                                              int component)
+{
+   MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
+               "entity set information");
+   GetEssentialTrueDofs(type, mesh->ent_sets->GetSetIndex(type, set_name),
+                        ess_tdof_list, component);
+}
+
 void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
                                              int component)
 {

fem/fespace.hpp

@@ -778,6 +778,19 @@ public:
                                   Array<int> &ess_vdofs,
                                   int component = -1) const;
 
+   /** Mark degrees of freedom associated with the entity set with the
+       specified entity type and set index. */
+   virtual void GetEssentialVDofs(EntitySets::EntityType type, int set_index,
+                                  Array<int> &ess_vdofs,
+                                  int component = -1) const;
+
+   /** Mark degrees of freedom associated with the entity set with the
+       specified entity type and set index. */
+   virtual void GetEssentialVDofs(EntitySets::EntityType type,
+                                  const std::string & set_name,
+                                  Array<int> &ess_vdofs,
+                                  int component = -1) const;
+
    /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess.
        For spaces with 'vdim' > 1, the 'component' parameter can be used
@@ -786,6 +799,19 @@ public:
                                      Array<int> &ess_tdof_list,
                                      int component = -1);
 
+   /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       entity set specified by the given entity type and set index. */
+   virtual void GetEssentialTrueDofs(EntitySets::EntityType type, int set_index,
+                                     Array<int> &ess_tdof_list,
+                                     int component = -1);
+
+   /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       entity set specified by the given entity type and set name. */
+   virtual void GetEssentialTrueDofs(EntitySets::EntityType type,
+                                     const std::string & set_name,
+                                     Array<int> &ess_tdof_list,
+                                     int component = -1);
+
    /** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
        with 'vdim' > 1, the 'component' parameter can be used to restricts the
        marked tDOFs to the specified component. Equivalent to

fem/pfespace.cpp

@@ -23,6 +23,8 @@
 #include <limits>
 #include <list>
 
+using namespace std;
+
 namespace mfem
 {
 
@@ -1018,6 +1020,30 @@ void ParFiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
    }
 }
 
+void ParFiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
+                                              int set_index,
+                                              Array<int> &ess_dofs,
+                                              int component) const
+{
+   FiniteElementSpace::GetEssentialVDofs(type, set_index, ess_dofs, component);
+
+   if (Conforming())
+   {
+      // Make sure that processors without boundary elements mark
+      // their boundary dofs (if they have any).
+      Synchronize(ess_dofs);
+   }
+}
+
+void ParFiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
+                                              const string & set_name,
+                                              Array<int> &ess_vdofs,
+                                              int component) const
+{
+   GetEssentialVDofs(type, pmesh->ent_sets->GetSetIndex(type, set_name),
+                     ess_vdofs, component);
+}
+
 void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
                                                  &bdr_attr_is_ess,
                                                  Array<int> &ess_tdof_list,
@@ -1047,6 +1073,27 @@ void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
    MarkerToList(true_ess_dofs, ess_tdof_list);
 }
 
+void ParFiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
+                                                 int set_index,
+                                                 Array<int> &ess_tdof_list,
+                                                 int component)
+{
+   Array<int> ess_dofs, true_ess_dofs;
+
+   GetEssentialVDofs(type, set_index, ess_dofs, component);
+   GetRestrictionMatrix()->BooleanMult(ess_dofs, true_ess_dofs);
+   MarkerToList(true_ess_dofs, ess_tdof_list);
+}
+
+void ParFiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
+                                                 const string & set_name,
+                                                 Array<int> &ess_tdof_list,
+                                                 int component)
+{
+   GetEssentialTrueDofs(type, pmesh->ent_sets->GetSetIndex(type, set_name),
+                        ess_tdof_list, component);
+}
+
 int ParFiniteElementSpace::GetLocalTDofNumber(int ldof) const
 {
    if (Nonconforming())

fem/pfespace.hpp

@@ -355,12 +355,38 @@ public:
                                   Array<int> &ess_dofs,
                                   int component = -1) const;
 
+   /** Mark degrees of freedom associated with the entity set with the
+       specified entity type and set index. */
+   virtual void GetEssentialVDofs(EntitySets::EntityType type, int set_index,
+                                  Array<int> &ess_vdofs,
+                                  int component = -1) const;
+
+   /** Mark degrees of freedom associated with the entity set with the
+       specified entity type and set index. */
+   virtual void GetEssentialVDofs(EntitySets::EntityType type,
+                                  const std::string & set_name,
+                                  Array<int> &ess_vdofs,
+                                  int component = -1) const;
+
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
    virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
                                      Array<int> &ess_tdof_list,
                                      int component = -1);
 
+   /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       entity set specified by the given entity type and set index. */
+   virtual void GetEssentialTrueDofs(EntitySets::EntityType type, int set_index,
+                                     Array<int> &ess_tdof_list,
+                                     int component = -1);
+
+   /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       entity set specified by the given entity type and set name. */
+   virtual void GetEssentialTrueDofs(EntitySets::EntityType type,
+                                     const std::string & set_name,
+                                     Array<int> &ess_tdof_list,
+                                     int component = -1);
+
    /** If the given ldof is owned by the current processor, return its local
        tdof number, otherwise return -1 */
    int GetLocalTDofNumber(int ldof) const;

general/sets.cpp

@@ -69,6 +69,7 @@ void IntegerSet::Recreate(const int n, const int *p)
 
    me.Sort();
 
+   // Remove duplicate entries
    for (j = 0, i = 1; i < n; i++)
       if (me[i] != me[j])
       {

general/sets.hpp

@@ -36,7 +36,7 @@ public:
    IntegerSet(const int n, const int *p) { Recreate(n, p); }
 
    /// Return the size of the set.
-   int Size() { return me.Size(); }
+   int Size() const { return me.Size(); }
 
    /// Return a reference to the sorted array of all the set entries.
    operator Array<int>& () { return me; }
@@ -50,6 +50,8 @@ public:
    /// Return 1 if the sets are equal and 0 otherwise.
    int operator==(IntegerSet &s);
 
+   inline const int & operator[](int i) const { return me[i]; }
+
    /** @brief Create an integer set from C-array 'p' of 'n' integers.
        Overwrites any existing set data. */
    void Recreate(const int n, const int *p);
@@ -64,7 +66,7 @@ private:
 public:
 
    /// Return the number of integer sets in the list.
-   int Size() { return TheList.Size(); }
+   int Size() const { return TheList.Size(); }
 
    /// Return the value of the first element of the ith set.
    int PickElementInSet(int i) { return TheList[i]->PickElement(); }
@@ -84,6 +86,9 @@ public:
    /// Write the list of sets into table 't'.
    void AsTable(Table &t);
 
+   inline const IntegerSet & operator[](int i) const { return *TheList[i]; }
+   inline IntegerSet & operator[](int i) { return *TheList[i]; }
+
    ~ListOfIntegerSets();
 };
 

general/stable3d.cpp

@@ -61,7 +61,7 @@ inline void Sort3 (int &r, int &c, int &f)
    }
 }
 
-int STable3D::Push (int r, int c, int f)
+int STable3D::Push (int r, int c, int f, int t)
 {
    STable3DNode *node;
 
@@ -86,6 +86,7 @@ int STable3D::Push (int r, int c, int f)
 #endif
    node->Column = c;
    node->Floor  = f;
+   node->Tier   = t;
    node->Number = NElem;
    node->Prev   = Rows[r];
    Rows[r] = node;
@@ -109,9 +110,9 @@ int STable3D::operator() (int r, int c, int f) const
          }
    }
 
-   MFEM_ABORT("(r,c,f) = (" << r << "," << c << "," << f << ")");
+   // MFEM_ABORT("(r,c,f) = (" << r << "," << c << "," << f << ")");
 
-   return 0;
+   return -1;
 }
 
 int STable3D::Index (int r, int c, int f) const
@@ -152,13 +153,13 @@ int STable3D::Push4 (int r, int c, int f, int t)
    switch (i)
    {
       case 0:
-         return Push (c,f,t);
+         return Push (c,f,t,r);
       case 1:
-         return Push (r,f,t);
+         return Push (r,f,t,c);
       case 2:
-         return Push (r,c,t);
+         return Push (r,c,t,f);
       case 3:
-         return Push (r,c,f);
+         return Push (r,c,f,t);
    }
 
    return -1;
@@ -218,6 +219,7 @@ void STable3D::Print(std::ostream & out) const
          out << row
              << ' ' << node_p->Column
              << ' ' << node_p->Floor
+             << ' ' << node_p->Tier
              << ' ' << node_p->Number
              << endl;
          node_p = node_p->Prev;

general/stable3d.hpp

@@ -15,6 +15,8 @@
 #include "mem_alloc.hpp"
 #include "../general/globals.hpp"
 
+#include <iostream>
+
 namespace mfem
 {
 
@@ -22,7 +24,7 @@ class STable3DNode
 {
 public:
    STable3DNode *Prev;
-   int Column, Floor, Number;
+   int Column, Floor, Tier, Number;
 };
 
 /** @brief Symmetric 3D Table stored as an array of rows each of which has a
@@ -47,7 +49,7 @@ public:
 
    /** @brief Check to see if this entry is in the table and add it to the table
        if it is not there. Returns the number assigned to the table entry. */
-   int Push (int r, int c, int f);
+   int Push (int r, int c, int f, int t = -1);
 
    /// Return the number assigned to the table entry. Abort if it's not there.
    int operator() (int r, int c, int f) const;
@@ -66,13 +68,30 @@ public:
        not there. */
    int operator() (int r, int c, int f, int t) const;
 
+   /// Return the number of rows added to the table.
+   int NumberOfRows() const { return Size; }
+
    /// Return the number of elements added to the table.
-   int NumberOfElements() { return NElem; }
+   int NumberOfElements() const { return NElem; }
 
    /// Print out all of the table elements.
    void Print(std::ostream &out = mfem::out) const;
 
    ~STable3D ();
+
+   class RowIterator
+   {
+   private:
+      STable3DNode *n;
+   public:
+      RowIterator (const STable3D &t, int r) { n = t.Rows[r]; }
+      int operator!() { return (n != NULL); }
+      void operator++() { n = n->Prev; }
+      int Column() { return (n->Column); }
+      int Floor()  { return (n->Floor); }
+      int Tier()   { return (n->Tier); }
+      int Index()  { return (n->Number); }
+   };
 };
 
 }

mesh/entsets.hpp

@@ -0,0 +1,219 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_ENTITY_SETS
+#define MFEM_ENTITY_SETS
+
+#include "../config/config.hpp"
+#include "../general/table.hpp"
+#include "../general/stable3d.hpp"
+#include <limits>
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+namespace mfem
+{
+
+class Mesh;
+class NCMesh;
+class NCEntitySets;
+
+class EntitySets
+{
+   friend class Mesh;
+   friend class NCMesh;
+   friend class NCEntitySets;
+
+public:
+   enum EntityType {INVALID = -1, VERTEX = 0, EDGE = 1, FACE = 2, ELEMENT = 3};
+
+   static std::map<EntityType,std::string> EntityTypeNames;
+
+   EntitySets(Mesh & mesh);
+   EntitySets(const EntitySets & ent_sets);
+   EntitySets(Mesh & mesh, NCMesh &ncmesh);
+
+   virtual ~EntitySets();
+
+   static const std::string & GetTypeName(EntityType t);
+
+   bool SetExists(EntityType t, unsigned int s) const;
+   bool SetExists(EntityType t, const std::string & s) const;
+
+   void Load(std::istream &input);
+   void Print(std::ostream &output) const;
+   virtual void PrintSetInfo(std::ostream &output) const;
+
+   inline Mesh *GetMesh() const { return mesh_; }
+
+   unsigned int GetNumSets(EntityType t) const;
+
+   const std::string & GetSetName(EntityType t, unsigned int s) const;
+   unsigned int        GetNumEntities(EntityType t, unsigned int s) const;
+
+   int          GetSetIndex(EntityType t, const std::string & s) const;
+   unsigned int GetNumEntities(EntityType t, const std::string & s) const;
+
+   inline std::set<int> & operator()(EntityType t, unsigned int s)
+   { return sets_[t][s]; }
+   inline const std::set<int> & operator()(EntityType t, unsigned int s) const
+   { return sets_[t][s]; }
+
+   const Table * GetEdgeVertexTable() const { return edge_vertex_; }
+   const Table * GetFaceVertexTable() const { return face_vertex_; }
+   const Table * GetFaceEdgeTable()   const { return face_edge_; }
+
+   // void Prune(int nelems);
+
+protected:
+
+   void SetNumSets(EntityType t, unsigned int n)
+   { sets_[t].resize(n); set_names_[t].resize(n); }
+   void SetSetName(EntityType t, int s, const std::string & name)
+   { set_names_[t][s] = name; set_index_by_name_[t][name] = s; }
+
+   /// Make local copies of edge_vertex, face_vertex, and face_edge tables.
+   void CopyMeshTables();
+
+   /// Refine quadrilateral mesh.
+   virtual void QuadUniformRefinement();
+
+   /// Refine hexahedral mesh.
+   virtual void HexUniformRefinement();
+
+   /// Refine 2D mesh.
+   virtual void UniformRefinement2D();
+
+   /// Refine 3D mesh.
+   virtual void UniformRefinement3D();
+
+private:
+
+   static void skip_comment_lines(std::istream &is, const char comment_char)
+   {
+      while (1)
+      {
+         is >> std::ws;
+         if (is.peek() != comment_char) { break; }
+         is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+      }
+   }
+   // Check for, and remove, a trailing '\r'.
+   static void filter_dos(std::string &line)
+   {
+      if (!line.empty() && *line.rbegin() == '\r')
+      { line.resize(line.size()-1); }
+   }
+
+   static std::map<EntityType,std::string> init_type_names();
+
+   void LoadEntitySets(std::istream &input, EntityType t,
+                       const std::string & header);
+
+   void PrintEntitySets(std::ostream &output, EntityType t,
+                        const std::string & header) const;
+
+   void PrintEdgeSets(std::ostream &output) const;
+
+   void PrintFaceSets(std::ostream &output) const;
+
+   void PrintEntitySetInfo(std::ostream & output, EntityType t,
+                           const std::string & ent_name) const;
+
+   void CopyEntitySets(const EntitySets & ent_sets, EntityType t);
+   void BuildEntitySets(NCMesh &ncmesh, EntityType t);
+
+protected:
+
+   Mesh     * mesh_;
+   Table    * edge_vertex_;
+   Table    * face_vertex_;
+   Table    * face_edge_;
+
+   int NumOfVertices_;
+   int NumOfEdges_;
+   int NumOfElements_;
+
+   /** The node/edge/face/element indices needed by the finite element
+       space to look up DoFs. */
+   std::vector<std::vector<std::set<int> > > sets_;
+
+   /// Names of each entity set
+   std::vector<std::vector<std::string> >    set_names_;
+
+   /// Indices of each entity set indexed by set name
+   std::vector<std::map<std::string, int> >  set_index_by_name_;
+};
+
+class NCEntitySets
+{
+   friend class EntitySets;
+
+public:
+   NCEntitySets(const EntitySets & ent_sets, NCMesh &ncmesh);
+   NCEntitySets(const NCEntitySets & ncent_sets);
+
+   bool SetExists(EntitySets::EntityType t, unsigned int s) const;
+   bool SetExists(EntitySets::EntityType t, const std::string & s) const;
+
+   unsigned int GetNumSets(EntitySets::EntityType t) const;
+
+   static int GetEntitySize(EntitySets::EntityType t);
+
+   const std::string & GetSetName(EntitySets::EntityType t, int s) const;
+   unsigned int        GetNumEntities(EntitySets::EntityType t, int s) const;
+   void                GetEntityIndex(EntitySets::EntityType t, int s,
+                                      int i, Array<int> & inds) const;
+
+   int          GetSetIndex(EntitySets::EntityType t,
+                            const std::string & s) const;
+   unsigned int GetNumEntities(EntitySets::EntityType t,
+                               const std::string & s) const;
+   void         GetEntityIndex(EntitySets::EntityType t,
+                               const std::string & s, int i,
+                               Array<int> & inds) const;
+
+   inline std::vector<int> & operator()(EntitySets::EntityType t, int s)
+   { return sets_[t][s]; }
+   inline const std::vector<int> & operator()(EntitySets::EntityType t,
+                                              int s) const
+   { return sets_[t][s]; }
+   inline int & operator()(EntitySets::EntityType t, int s, int i)
+   { return sets_[t][s][i]; }
+   inline int operator()(EntitySets::EntityType t, int s, int i) const
+   { return sets_[t][s][i]; }
+
+private:
+   void CopyNCEntitySets(const NCEntitySets & ncent_sets,
+                         EntitySets::EntityType t);
+
+protected:
+
+   NCMesh * ncmesh_;
+
+   /// The nodes defining the node/edge/face/element sets
+   std::vector<std::vector<std::vector<int> > > sets_;
+
+   /// Names of each entity set
+   std::vector<std::vector<std::string> >      set_names_;
+
+   /// Indices of each entity set indexed by set name
+   std::vector<std::map<std::string, int> >    set_index_by_name_;
+
+   /// Number of indices per entity
+   static const int entity_size_[4];
+};
+
+} // namespace mfem
+
+#endif // MFEM_ENTITY_SETS

mesh/mesh.cpp

@@ -1177,13 +1177,15 @@ void Mesh::Init()
    own_nodes = 1;
    NURBSext = NULL;
    ncmesh = NULL;
+   ent_sets = NULL;
    last_operation = Mesh::NONE;
 }
 
 void Mesh::InitTables()
 {
    el_to_edge =
-      el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
+      el_to_face = el_to_el = bel_to_edge = face_edge =
+                                               face_vertex = edge_vertex = NULL;
 }
 
 void Mesh::SetEmpty()
@@ -1205,6 +1207,7 @@ void Mesh::DestroyTables()
    }
 
    delete face_edge;
+   delete face_vertex;
    delete edge_vertex;
 }
 
@@ -1212,6 +1215,8 @@ void Mesh::DestroyPointers()
 {
    if (own_nodes) { delete Nodes; }
 
+   delete ent_sets;
+
    delete ncmesh;
 
    delete NURBSext;
@@ -3346,6 +3351,12 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    // Copy the edge-to-vertex Table, edge_vertex
    edge_vertex = (mesh.edge_vertex) ? new Table(*mesh.edge_vertex) : NULL;
 
+   // Copy the face-to-vertex Table, edge_vertex
+   face_vertex = (mesh.face_vertex) ? new Table(*mesh.face_vertex) : NULL;
+
+   // Do not copy any of the coarse (c_*), fine (f_*) or fine/coarse (fc_*)
+   // data members.
+
    // Copy the attributes and bdr_attributes
    mesh.attributes.Copy(attributes);
    mesh.bdr_attributes.Copy(bdr_attributes);
@@ -3396,6 +3407,9 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
       Nodes = mesh.Nodes;
       own_nodes = 0;
    }
+
+   // Copy entity sets if present in the input mesh
+   ent_sets = (mesh.ent_sets) ? new EntitySets(*mesh.ent_sets) : NULL;
 }
 
 Mesh::Mesh(Mesh &&mesh) : Mesh()
@@ -5768,6 +5782,38 @@ Table *Mesh::GetEdgeVertexTable() const
    return edge_vertex;
 }
 
+Table *Mesh::GetFaceVertexTable() const
+{
+   if (face_vertex)
+   {
+      return face_vertex;
+   }
+
+   STable3D * faces_tbl = GetFacesTable();
+
+   int nfaces = faces_tbl->NumberOfElements();
+   face_vertex = new Table(nfaces, 4);
+   for (int i = 0; i < NumOfVertices; i++)
+   {
+      for (STable3D::RowIterator it(*faces_tbl, i); !it; ++it)
+      {
+         int j = it.Index();
+         face_vertex->Push(j, i);
+         face_vertex->Push(j, it.Column());
+         face_vertex->Push(j, it.Floor());
+         if ( it.Tier() > 0 )
+         {
+            face_vertex->Push(j, it.Tier());
+         }
+      }
+   }
+   face_vertex->Finalize();
+
+   delete faces_tbl;
+
+   return face_vertex;
+}
+
 Table *Mesh::GetVertexToElementTable()
 {
    int i, j, nv, *v;
@@ -6402,7 +6448,7 @@ void Mesh::GenerateNCFaceInfo()
    }
 }
 
-STable3D *Mesh::GetFacesTable()
+STable3D *Mesh::GetFacesTable() const
 {
    STable3D *faces_tbl = new STable3D(NumOfVertices);
    for (int i = 0; i < NumOfElements; i++)
@@ -7657,6 +7703,11 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
       NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
 
+   if ( ent_sets )
+   {
+      ent_sets->CopyMeshTables();
+   }
+
    int quad_counter = 0;
    for (int i = 0; i < NumOfElements; i++)
    {
@@ -7792,6 +7843,11 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
 
    if (update_nodes) { UpdateNodes(); }
 
+   if ( ent_sets )
+   {
+      ent_sets->UniformRefinement2D();
+   }
+
 #ifdef MFEM_DEBUG
    if (!Nodes || update_nodes)
    {
@@ -7822,6 +7878,11 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
       GetElementToFaceTable();
    }
 
+   if ( ent_sets )
+   {
+      ent_sets->CopyMeshTables();
+   }
+
    Array<int> f2qf_loc;
    Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
    f2qf.SetSize(0);
@@ -8148,7 +8209,6 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
                }
                AverageVertices(vv, 4, oface + f2qf[f[fi]]);
             }
-
             for (int ei = 0; ei < 9; ei++)
             {
                for (int k = 0; k < 2; k++)
@@ -8492,6 +8552,11 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    sequence++;
 
    if (update_nodes) { UpdateNodes(); }
+
+   if (ent_sets)
+   {
+      ent_sets->UniformRefinement3D();
+   }
 }
 
 void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
@@ -8961,6 +9026,8 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
 
    mfem::Swap(geom_factors, other.geom_factors);
 
+   mfem::Swap(ent_sets, other.ent_sets);
+
 #ifdef MFEM_USE_MEMALLOC
    TetMemory.Swap(other.TetMemory);
 #endif

mesh/mesh.hpp

@@ -20,6 +20,7 @@
 #include "vertex.hpp"
 #include "vtk.hpp"
 #include "ncmesh.hpp"
+#include "entsets.hpp"
 #include "../fem/eltrans.hpp"
 #include "../fem/coefficient.hpp"
 #include "../general/zstr.hpp"
@@ -54,9 +55,11 @@ class Mesh
 #ifdef MFEM_USE_MPI
    friend class ParMesh;
    friend class ParNCMesh;
+   friend class ParEntitySets;
 #endif
    friend class NCMesh;
    friend class NURBSExtension;
+   friend class EntitySets;
 
 #ifdef MFEM_USE_ADIOS2
    friend class adios2stream;
@@ -166,6 +169,7 @@ protected:
    Array<int> be_to_face;
    mutable Table *face_edge;
    mutable Table *edge_vertex;
+   mutable Table *face_vertex;
 
    IsoparametricTransformation Transformation, Transformation2;
    IsoparametricTransformation BdrTransformation;
@@ -216,6 +220,8 @@ public:
    Array<FaceGeometricFactors*>
    face_geom_factors; ///< Optional face geometric factors.
 
+   EntitySets *ent_sets;
+
    // Global parameter that can be used to control the removal of unused
    // vertices performed when reading a mesh in MFEM format. The default value
    // (true) is set in mesh_readers.cpp.
@@ -287,7 +293,7 @@ protected:
    void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert);
    void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert);
 
-   STable3D *GetFacesTable();
+   STable3D *GetFacesTable() const;
    STable3D *GetElementToFaceTable(int ret_ftbl = 0);
 
    /** Red refinement. Element with index i is refined. The default
@@ -1067,9 +1073,12 @@ public:
    /// Returns the face-to-edge Table (3D)
    Table *GetFaceEdgeTable() const;
 
-   /// Returns the edge-to-vertex Table (3D)
+   /// Returns the edge-to-vertex Table (2D or 3D)
    Table *GetEdgeVertexTable() const;
 
+   /// Returns the face-to-vertex Table (2d or 3D)
+   Table *GetFaceVertexTable() const;
+
    /// Return the indices and the orientations of all faces of element i.
    void GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const;
 

mesh/mesh_readers.cpp

@@ -100,6 +100,14 @@ void Mesh::ReadMFEMMesh(std::istream &input, int version, int &curved)
       curved = 1;
    }
 
+   ent_sets = new EntitySets(*this);
+   ent_sets->Load(input);
+   if ( ent_sets->GetNumSets(EntitySets::FACE) > 0 && faces.Size() == 0 )
+   {
+      GetElementToFaceTable();
+      GenerateFaces();
+   }
+
    // When visualizing solutions on non-conforming grids, PETSc
    // may dump additional vertices
    if (remove_unused_vertices) { RemoveUnusedVertices(); }

mesh/ncmesh.cpp

@@ -185,6 +185,10 @@ NCMesh::NCMesh(const Mesh *mesh)
       face->attribute = be->GetAttribute();
    }
 
+   // Store entity set information if present in the Mesh
+   ncent_sets = (mesh->ent_sets) ?
+                new NCEntitySets(*mesh->ent_sets, *this) : NULL;
+
    // copy top-level vertex coordinates (leave empty if the mesh is curved)
    if (!mesh->Nodes)
    {
@@ -216,6 +220,10 @@ NCMesh::NCMesh(const NCMesh &other)
    other.free_element_ids.Copy(free_element_ids);
    other.root_state.Copy(root_state);
    other.coordinates.Copy(coordinates);
+
+   // Copy the entity set information
+   ncent_sets = (other.ncent_sets) ? new NCEntitySets(*other.ncent_sets) : NULL;
+
    Update();
 }
 
@@ -254,8 +262,11 @@ NCMesh::~NCMesh()
          DeleteUnusedFaces(elemFaces);
       }
    }
+
    // NOTE: in release mode, we just throw away all faces and nodes at once
 #endif
+
+   delete ncent_sets;
 }
 
 NCMesh::Node::~Node()
@@ -2504,6 +2515,42 @@ void NCMesh::OnMeshUpdated(Mesh *mesh)
       if (face->index < 0) { face->index = NFaces + (nghosts++); }
    }
    MFEM_ASSERT(nghosts == NGhostFaces, "");
+
+   if (ncent_sets)
+   {
+      std::cout << "NCMesh::OnMeshUpdated ncent_sets is non NULL" << std::endl;
+      if (!mesh->ent_sets)
+      {
+         std::cout << "NCMesh::OnMeshUpdated creating ent_sets from NCMesh" << std::endl;
+         mesh->ent_sets = new EntitySets(*mesh, *this);
+         std::cout << "NCMesh::OnMeshUpdated done creating ent_sets from NCMesh" <<
+                   std::endl;
+      }
+   }
+
+   std::ostringstream ossN;
+   ossN << "node_on_mesh_updated.out";
+   std::ofstream ofsN(ossN.str().c_str());
+   ofsN << nodes.Size() << std::endl;
+   for (int i=0; i<nodes.Size(); i++)
+   {
+      ofsN << i
+           // << " " << nodes[i].vert_refc
+           //  << " " << nodes[i].edge_refc
+           << " " << nodes[i].HasVertex()
+           << " " << nodes[i].HasEdge()
+           << " " << nodes[i].vert_index
+           << " " << nodes[i].edge_index
+           << " " << nodes[i].p1
+           << " " << nodes[i].p2
+           << " " << nodes[i].next << std::endl;
+   }
+   ofsN.close();
+
+   NEdges = mesh->GetNEdges();
+   NFaces = mesh->GetNumFaces();
+
+   std::cout << "Leaving NCMesh::OnMeshUpdated" << std::endl;
 }
 
 
@@ -3306,12 +3353,15 @@ const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
 {
    int mid = nodes.FindId(v0, v1);
-   if (mid >= 0 && nodes[mid].HasVertex())
+   if (mid >= 0)
    {
-      indices.Append(mid);
+      if (nodes[mid].HasVertex())
+      {
+         indices.Append(mid);
 
-      CollectEdgeVertices(v0, mid, indices);
-      CollectEdgeVertices(mid, v1, indices);
+         CollectEdgeVertices(v0, mid, indices);
+         CollectEdgeVertices(mid, v1, indices);
+      }
    }
 }
 
@@ -3373,6 +3423,78 @@ void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
    }
 }
 
+void NCMesh::CollectElementVertices(int elem_id, Array<int> &indices)
+{
+   Element &el = elements[elem_id];
+
+   if (el.ref_type != 0)
+   {
+      // This element has been refined so recurse into its children
+      for (int i = 0; i < 8; i++)
+      {
+         if (el.child[i] >= 0 && el.child[i] < elements.Size())
+         {
+            CollectElementVertices(el.child[i], indices);
+         }
+      }
+   }
+   else
+   {
+      // This element has not been refined so add its vertices
+      for (int i=0; i<8; i++)
+      {
+         if (el.node[i] >= 0 && el.node[i] < nodes.Size())
+         {
+            indices.Append(el.node[i]);
+         }
+      }
+   }
+}
+
+void NCMesh::CollectElementEdges(int elem_id, Array<int> &indices)
+{
+   Element &el = elements[elem_id];
+
+   if (el.ref_type != 0)
+   {
+      // This element has been refined so recurse into its children
+      for (int i = 0; i < 8; i++)
+      {
+         if (el.child[i] >= 0 && el.child[i] < elements.Size())
+         {
+            CollectElementEdges(el.child[i], indices);
+         }
+      }
+   }
+   else
+   {
+      int* node = el.node;
+      GeomInfo& gi = GI[(int) el.geom];
+
+      for (int i = 0; i < gi.nv; i++)
+      {
+         if (nodes[node[i]].HasEdge())
+         {
+            indices.Append(node[i]);
+         }
+      }
+
+      for (int i = 0; i < gi.ne; i++)
+      {
+         const int* ev = gi.edges[i];
+         int index = nodes.FindId(node[ev[0]], node[ev[1]]);
+
+         if (index >= 0)
+         {
+            if (nodes[index].HasEdge())
+            {
+               indices.Append(index);
+            }
+         }
+      }
+   }
+}
+
 void NCMesh::BuildElementToVertexTable()
 {
    int nrows = leaf_elements.Size();
@@ -4874,6 +4996,107 @@ int NCMesh::GetElementDepth(int i) const
    return depth;
 }
 
+void NCMesh::GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edges)
+{
+   std::cout << "entering NCMesh::GetRefinedEdges "
+             <<"searching for edge with vertices: " << vn0 << " and " << vn1
+             << std::endl;
+   int mid = nodes.FindId(vn0, vn1);
+   if (mid < 0) { return; }
+
+   Node &nd = nodes[mid];
+
+   // if ( nd.edge_index < 0 ) { return; }
+
+   // edges.Append(nd.edge_index);
+   if ( nd.HasEdge() )
+   {
+      std::cout << " found node " << mid << std::endl;
+      edges.Append(mid);
+   }
+
+   GetRefinedEdges(vn0, mid, edges);
+   GetRefinedEdges(mid, vn1, edges);
+}
+
+void NCMesh::GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
+                             BlockArray<int> & face_ids)
+{
+   // Face* fa = faces.Find(vn0, vn1, vn2, vn3);
+   int face = faces.FindId(vn0, vn1, vn2, vn3);
+   /*
+   if (fa)
+   {
+      if ( fa->index >= 0 )
+      {
+         face_ids.Append(fa->index);
+      }
+      return;
+   }
+   */
+   if (face>=0)
+   {
+      if ( faces[face].index >= 0 )
+      {
+         face_ids.Append(face);
+      }
+      return;
+   }
+
+   // we need to recurse deeper
+   int mid[4];
+   int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, mid);
+
+   if (split == 1) // "X" split face
+   {
+      GetRefinedFaces(vn0, mid[0], mid[2], vn3, face_ids);
+      GetRefinedFaces(mid[0], vn1, vn2, mid[2], face_ids);
+   }
+   else if (split == 2) // "Y" split face
+   {
+      GetRefinedFaces(vn0, vn1, mid[1], mid[3], face_ids);
+      GetRefinedFaces(mid[3], mid[1], vn2, vn3, face_ids);
+   }
+}
+
+void NCMesh::GetRefinedElements(int elem_id, BlockArray<int> & elem_ids)
+{
+   // std::cout << "entering NCMesh::GetRefinedElements searching for element id: "
+   //             << elem_id << std::endl;
+   Element &el = elements[elem_id];
+   /*
+   if (el.index >= 0 && el.rank >= 0)
+   {
+      elem_ids.Append(el.index);
+      return;
+   }
+
+   for (int i = 0; i < 8; i++)
+   {
+      if (el.child[i] >= 0 && el.child[i] < elements.Size() )
+      {
+         GetRefinedElements(el.child[i], elem_ids);
+      }
+   }
+   */
+   if (el.ref_type != 0)
+   {
+      // This element has been refined so recurse into its children
+      for (int i = 0; i < 8; i++)
+      {
+         if (el.child[i] >= 0 && el.child[i] < elements.Size() )
+         {
+            GetRefinedElements(el.child[i], elem_ids);
+         }
+      }
+   }
+   else
+   {
+      // This element has not been refined so add it
+      elem_ids.Append(elem_id);
+   }
+}
+
 int NCMesh::GetElementSizeReduction(int i) const
 {
    int elem = leaf_elements[i];
@@ -4993,6 +5216,183 @@ void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
    bdr_edges.Unique();
 }
 
+void NCMesh::GetEntitySetClosure(EntitySets::EntityType type,
+                                 int set_index,
+                                 Array<int> &es_vertices,
+                                 Array<int> &es_edges,
+                                 Array<int> &es_faces)
+{
+   es_vertices.SetSize(0);
+   es_edges.SetSize(0);
+   es_faces.SetSize(0);
+
+   MFEM_VERIFY(ncent_sets != NULL, "NCMesh object contains no "
+               "entity set information");
+   if (!ncent_sets->SetExists(type, set_index))
+   {
+      std::ostringstream oss; oss << "Entity set of type \""
+                                  << EntitySets::GetTypeName(type)
+                                  << "\" and index " << set_index
+                                  << " was not found.";
+
+      MFEM_VERIFY(false, oss.str().c_str());
+   }
+
+   int ni = ncent_sets->GetNumEntities(type ,set_index);
+   Array<int> inds;
+   Array<int> coll_inds;
+
+   switch (type)
+   {
+      case EntitySets::VERTEX:
+      {
+         /// Do nothing because vertices cannot hide
+      }
+      break;
+      case EntitySets::EDGE:
+      {
+         for (int i=0; i<ni; i++)
+         {
+            ncent_sets->GetEntityIndex(type, set_index, i, inds);
+
+            // collect vertices
+            inds.Copy(coll_inds);
+            this->CollectEdgeVertices(inds[0], inds[1], coll_inds);
+            for (int j=0; j<coll_inds.Size(); j++)
+            {
+               int index = nodes[coll_inds[j]].vert_index;
+               if (index >= 0)
+               {
+                  es_vertices.Append(index);
+               }
+            }
+         }
+      }
+      break;
+      case EntitySets::FACE:
+      {
+         for (int i=0; i<ni; i++)
+         {
+            ncent_sets->GetEntityIndex(type, set_index, i, inds);
+
+            // collect vertices
+            inds.Copy(coll_inds);
+            if (inds.Size() == 4)
+            {
+               this->CollectQuadFaceVertices(inds[0], inds[1], inds[2], inds[3],
+                                             coll_inds);
+            }
+            else
+            {
+               this->CollectTriFaceVertices(inds[0], inds[1], inds[2],
+                                            coll_inds);
+            }
+            for (int j=0; j<coll_inds.Size(); j++)
+            {
+               int index = nodes[coll_inds[j]].vert_index;
+               if (index >= 0)
+               {
+                  es_vertices.Append(index);
+               }
+            }
+         }
+      }
+      break;
+      case EntitySets::ELEMENT:
+      {
+         for (int i=0; i<ni; i++)
+         {
+            int elem_id = (*ncent_sets)(type, set_index, i);
+            std::cout << "examining element " << elem_id << std::endl;
+
+            // collect vertices
+            coll_inds.SetSize(0);
+            this->CollectElementVertices(elem_id, coll_inds);
+            for (int j=0; j<coll_inds.Size(); j++)
+            {
+               int index = nodes[coll_inds[j]].vert_index;
+               if (index >= 0)
+               {
+                  es_vertices.Append(index);
+               }
+            }
+
+            // collect edges
+            coll_inds.SetSize(0);
+            this->CollectElementEdges(elem_id, coll_inds);
+            for (int j=0; j<coll_inds.Size(); j++)
+            {
+               int index = nodes[coll_inds[j]].edge_index;
+               if (index >= 0)
+               {
+                  es_edges.Append(index);
+               }
+            }
+
+         }
+      }
+      break;
+      default:
+         MFEM_ABORT("GetEnitySetClosure - Unknown entity set type: \""
+                    << EntitySets::GetTypeName(type) << "\"");
+   }
+   /*
+   if (Dim == 3)
+   {
+      GetFaceList(); // make sure 'boundary_faces' is up to date
+
+      for (int i = 0; i < boundary_faces.Size(); i++)
+      {
+         int face = boundary_faces[i];
+         if (bdr_attr_is_ess[faces[face].attribute - 1])
+         {
+            int node[4];
+            FindFaceNodes(face, node);
+
+            for (int j = 0; j < 4; j++)
+            {
+               bdr_vertices.Append(nodes[node[j]].vert_index);
+
+               int enode = nodes.FindId(node[j], node[(j+1) % 4]);
+               MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
+               bdr_edges.Append(nodes[enode].edge_index);
+
+               while ((enode = GetEdgeMaster(enode)) >= 0)
+               {
+                  // append master edges that may not be accessible from any
+                  // boundary element, this happens in 3D in re-entrant corners
+                  bdr_edges.Append(nodes[enode].edge_index);
+               }
+            }
+         }
+      }
+   }
+   else if (Dim == 2)
+   {
+      GetEdgeList(); // make sure 'boundary_faces' is up to date
+
+      for (int i = 0; i < boundary_faces.Size(); i++)
+      {
+         int face = boundary_faces[i];
+         Face &fc = faces[face];
+         if (bdr_attr_is_ess[fc.attribute - 1])
+         {
+            bdr_vertices.Append(nodes[fc.p1].vert_index);
+            bdr_vertices.Append(nodes[fc.p3].vert_index);
+         }
+      }
+   }
+   */
+   es_vertices.Sort();
+   es_vertices.Unique();
+
+   es_edges.Sort();
+   es_edges.Unique();
+
+   es_faces.Sort();
+   es_faces.Unique();
+}
+
 static int max4(int a, int b, int c, int d)
 {
    return std::max(std::max(a, b), std::max(c, d));

mesh/ncmesh.hpp

@@ -19,6 +19,7 @@
 #include "../linalg/densemat.hpp"
 #include "element.hpp"
 #include "vertex.hpp"
+#include "entsets.hpp"
 #include "../fem/geom.hpp"
 
 #include <vector>
@@ -117,6 +118,9 @@ struct MatrixMap; // for internal use
  */
 class NCMesh
 {
+   friend class EntitySets;
+   friend class NCEntitySets;
+
 public:
    //// Initialize with elements from an existing 'mesh'.
    explicit NCMesh(const Mesh *mesh);
@@ -343,6 +347,16 @@ public:
                                    Array<int> &bdr_vertices,
                                    Array<int> &bdr_edges);
 
+   /** Get a list of vertices (2D/3D), edges (2D/3D), and faces (3D) that
+       coincide with members of the specified entity set. In 3D this function
+       also reveals "hidden" edges or faces. In parallel it helps identifying
+       vertices/edges/faces affected by non-local entities. */
+   virtual void GetEntitySetClosure(EntitySets::EntityType t,
+                                    int set_index,
+                                    Array<int> &es_vertices,
+                                    Array<int> &es_edges,
+                                    Array<int> &es_faces);
+
    /// Return element geometry type. @a index is the Mesh element number.
    Geometry::Type GetElementGeometry(int index) const
    { return elements[leaf_elements[index]].Geom(); }
@@ -357,6 +371,19 @@ public:
    /// Return the distance of leaf 'i' from the root.
    int GetElementDepth(int i) const;
 
+   /** Collect edge indices of all refined edges which are children of
+       the coarse edge defined by the given vertices. */
+   void GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edge_ids);
+
+   /** Collect face indices of all refined faces which are children of
+       the coarse face defined by the given vertices. */
+   void GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
+                        BlockArray<int> & face_ids);
+
+   /** Collect element indices of all refined elements which are children of
+       the coarse element defined by the given element index. */
+   void GetRefinedElements(int elem_id, BlockArray<int> & elem_ids);
+
    /** Return the size reduction compared to the root element (ignoring local
        stretching and curvature). */
    int GetElementSizeReduction(int i) const;
@@ -501,6 +528,7 @@ protected: // implementation
    Array<double> coordinates;
 
 
+
    // secondary data
 
    /** Apart from the primary data structure, which is the element/node/face
@@ -530,6 +558,8 @@ protected: // implementation
 
    Table element_vertex; ///< leaf-element to vertex table, see FindSetNeighbors
 
+   // Node/edge/Face/Element sets defined on the coarse mesh
+   NCEntitySets * ncent_sets;
 
    void UpdateLeafElements();
    void UpdateVertices(); ///< update Vertex::index and vertex_nodeId
@@ -711,6 +741,10 @@ protected: // implementation
    void CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices);
    void CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
                                 Array<int> &indices);
+   void CollectElementVertices(int elem_id, Array<int> &indices);
+
+   void CollectElementEdges(int elem_id, Array<int> &indices);
+
    void BuildElementToVertexTable();
 
    void UpdateElementToVertexTable()
@@ -926,6 +960,7 @@ public:
 #endif
 
    friend class ParNCMesh; // for ParNCMesh::ElementSet
+   friend class ParNCEntitySets;
    friend struct MatrixMap;
    friend struct PointMatrixHash;
 };

mesh/pentsets.cpp

@@ -0,0 +1,392 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "pentsets.hpp"
+#include "pmesh.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+ParEntitySets::ParEntitySets(const ParEntitySets & ent_sets)
+   : EntitySets(ent_sets),
+     pmesh_(ent_sets.GetParMesh())
+{
+   MPI_Comm_size(pmesh_->GetComm(), &NRanks_);
+   MPI_Comm_rank(pmesh_->GetComm(), &MyRank_);
+   cout << MyRank_ << ": Entering ParEntitySets copy c'tor" << endl;
+   cout << MyRank_ << ": Leaving ParEntitySets copy c'tor" << endl;
+}
+
+ParEntitySets::ParEntitySets(ParMesh & pmesh, const EntitySets & ent_sets,
+                             int * partitioning,
+                             const Array<int> & vert_global_local)
+   : EntitySets(ent_sets),
+     pmesh_(&pmesh)
+{
+   // The copy constructor for EntitySets will initialize this object's
+   // data with the correct set names, and numbers of sets.  However,
+   // the set entries themselves will need to be recomputed based on
+   // local numberings and the paritioning.
+   //
+   // The EntitySets object will be a copy of the serial object.  This
+   // constructor will have to prune and renumber the data.  Once this
+   // is done the mesh pointer stored in the EntitySets object can be
+   // replaced with the local portion of the parallel mesh.
+
+   MPI_Comm MyComm = pmesh_->GetComm();
+
+   MPI_Comm_size(MyComm, &NRanks_);
+   MPI_Comm_rank(MyComm, &MyRank_);
+   cout << MyRank_ << ": Entering ParEntitySets(ParMesh, EntitySets, ...) c'tor" <<
+        endl;
+
+   int nelem = mesh_->GetNE();
+
+   DSTable v_to_v(vert_global_local.Size());
+   pmesh_->GetVertexToVertexTable(v_to_v);
+
+   STable3D * faces_tbl = NULL;
+
+   const Table * serial_edge_vertex = NULL;
+   const Table * serial_face_vertex = NULL;
+
+   if ( ent_sets.GetNumSets(EDGE) > 0 )
+   {
+      serial_edge_vertex = ent_sets.GetEdgeVertexTable();
+   }
+   if ( ent_sets.GetNumSets(FACE) > 0 )
+   {
+      serial_face_vertex = ent_sets.GetFaceVertexTable();
+      faces_tbl = pmesh_->GetFacesTable();
+   }
+
+   Array<int> elem_global_local(nelem);
+   elem_global_local = -1;
+   int elem_counter = 0;
+   for (int i=0; i<nelem; i++)
+   {
+      if ( partitioning[i] == MyRank_ )
+      {
+         elem_global_local[i] = elem_counter;
+         elem_counter++;
+      }
+   }
+
+   EntityType t;
+   unsigned int ns;
+
+   t = VERTEX;
+   ns = ent_sets.GetNumSets(t);
+   for (unsigned int s=0; s<ns; s++)
+   {
+      set<int>::iterator it;
+      sets_[t][s].clear();
+      for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
+      {
+         int v0 = vert_global_local[*it];
+         if ( v0 >= 0 )
+         {
+            sets_[t][s].insert(v0);
+         }
+      }
+   }
+
+   if ( pmesh_->Dimension() > 1 )
+   {
+      t = EDGE;
+      ns = ent_sets.GetNumSets(t);
+      for (unsigned int s=0; s<ns; s++)
+      {
+         set<int>::iterator it;
+         sets_[t][s].clear();
+         for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
+         {
+            int old_edge = *it;
+            const int *v = serial_edge_vertex->GetRow(old_edge);
+            int v0 = vert_global_local[v[0]];
+            int v1 = vert_global_local[v[1]];
+            if ( v0 >= 0 && v1 >= 0 )
+            {
+               int new_edge = v_to_v(v0,v1);
+               if ( new_edge >= 0 )
+               {
+                  sets_[t][s].insert(new_edge);
+               }
+            }
+         }
+      }
+   }
+
+   if ( pmesh_->Dimension() > 2 )
+   {
+      Array<int> v;
+      t = FACE;
+      ns = ent_sets.GetNumSets(t);
+      for (unsigned int s=0; s<ns; s++)
+      {
+         set<int>::iterator it;
+         sets_[t][s].clear();
+         for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
+         {
+            int old_face = *it;
+            int numv = serial_face_vertex->RowSize(old_face);
+            const int *v = serial_face_vertex->GetRow(old_face);
+            if ( vert_global_local[v[0]] >= 0 &&
+                 vert_global_local[v[1]] >= 0 &&
+                 vert_global_local[v[2]] >= 0 )
+            {
+               int new_face = -1;
+               if ( numv == 3 )
+               {
+                  new_face = (*faces_tbl)(vert_global_local[v[0]],
+                                          vert_global_local[v[1]],
+                                          vert_global_local[v[2]]);
+               }
+               else
+               {
+                  new_face = (*faces_tbl)(vert_global_local[v[0]],
+                                          vert_global_local[v[1]],
+                                          vert_global_local[v[2]],
+                                          vert_global_local[v[3]]);
+               }
+               if ( new_face >= 0 )
+               {
+                  sets_[t][s].insert(new_face);
+               }
+            }
+         }
+      }
+      delete faces_tbl;
+   }
+
+   t = ELEMENT;
+   ns = ent_sets.GetNumSets(t);
+   for (unsigned int s=0; s<ns; s++)
+   {
+      set<int>::iterator it;
+      sets_[t][s].clear();
+      for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
+      {
+         if ( partitioning[*it] == MyRank_ )
+         {
+            sets_[t][s].insert(elem_global_local[*it]);
+         }
+      }
+   }
+
+   this->mesh_ = (Mesh*)this->pmesh_;
+
+   this->CopyMeshTables();
+   cout << MyRank_ << ": Leaving ParEntitySets(ParMesh, EntitySets, ...) c'tor" <<
+        endl;
+}
+
+ParEntitySets::ParEntitySets(ParMesh & pmesh, ParNCMesh &pncmesh)
+   : EntitySets(pmesh),
+     pmesh_(&pmesh)
+{
+   MPI_Comm MyComm = pmesh_->GetComm();
+
+   MPI_Comm_size(MyComm, &NRanks_);
+   MPI_Comm_rank(MyComm, &MyRank_);
+   cout << MyRank_ << ": Entering ParEntitySets(ParMesh, ParNCMesh) c'tor" << endl;
+
+   this->BuildEntitySets(pncmesh, VERTEX);
+   this->BuildEntitySets(pncmesh, EDGE);
+   this->BuildEntitySets(pncmesh, FACE);
+   this->BuildEntitySets(pncmesh, ELEMENT);
+   cout << MyRank_ << ": Leaving ParEntitySets(ParMesh, ParNCMesh) c'tor" << endl;
+}
+
+ParEntitySets::~ParEntitySets()
+{
+   cout << MyRank_ << ": Entering ParEntitySets d'tor" << endl;
+   cout << MyRank_ << ": Leaving ParEntitySets d'tor" << endl;
+}
+
+void
+ParEntitySets::PrintSetInfo(std::ostream & output) const
+{
+   if ( MyRank_ == 0 &&
+        ( GetNumSets(VERTEX) > 0 || GetNumSets(EDGE)    > 0 ||
+          GetNumSets(FACE)   > 0 || GetNumSets(ELEMENT) > 0 ) )
+   {
+      output << "\nMFEM Parallel Entity Sets:\n";
+   }
+   this->PrintEntitySetInfo(output, VERTEX,  "Vertex");
+   this->PrintEntitySetInfo(output, EDGE,    "Edge");
+   this->PrintEntitySetInfo(output, FACE,    "Face");
+   this->PrintEntitySetInfo(output, ELEMENT, "Element");
+}
+
+void
+ParEntitySets::PrintEntitySetInfo(std::ostream & output, EntityType t,
+                                  const string & ent_name) const
+{
+   if ( sets_[t].size() > 0 )
+   {
+      if ( MyRank_ == 0 )
+      {
+         output << "  " << ent_name
+                << " Sets (Index, Set Name, Global Size):\n";
+      }
+      for (unsigned int s=0; s<sets_[t].size(); s++)
+      {
+         int loc_size = sets_[t][s].size();
+         int glb_size = -1;
+         MPI_Reduce(&loc_size, &glb_size, 1, MPI_INT, MPI_SUM, 0,
+                    pmesh_->GetComm());
+         if ( MyRank_ == 0 )
+         {
+            output << '\t' << s
+                   << '\t' << set_names_[t][s]
+                   << '\t' << glb_size
+                   << '\n';
+         }
+      }
+      if ( MyRank_ == 0 )
+      {
+         output << '\n';
+      }
+   }
+}
+
+void
+ParEntitySets::BuildEntitySets(ParNCMesh &pncmesh, EntityType t)
+{
+   cout << MyRank_ << ": BuildEntitySets for type " << GetTypeName(t) << endl;
+   int es = pncmesh.pncent_sets->GetEntitySize(t);
+   unsigned int ns = pncmesh.pncent_sets->GetNumSets(t);
+   cout << MyRank_ << ": num sets " << ns << endl;
+
+   Array<int> inds(es);
+
+   sets_[t].resize(ns);
+   set_names_[t].resize(ns);
+   for (unsigned int s=0; s<ns; s++)
+   {
+      int ni = pncmesh.pncent_sets->GetNumEntities(t, s);
+      set_names_[t][s] = pncmesh.pncent_sets->GetSetName(t, s);
+      set_index_by_name_[t][set_names_[t][s]] = s;
+
+      switch (t)
+      {
+         case VERTEX:
+            for (int i=0; i<ni; i++)
+            {
+               int node = (*pncmesh.pncent_sets)(t, s, i);
+               int index = pncmesh.nodes[node].vert_index;
+               if (!pncmesh.IsGhost(0,index))
+               {
+                  sets_[t][s].insert(index);
+               }
+            }
+            break;
+         case EDGE:
+            for (int i=0; i<ni; i++)
+            {
+               pncmesh.pncent_sets->GetEntityIndex(t, s, i, inds);
+               BlockArray<int> ind_coll;
+               pncmesh.GetRefinedEdges(inds[0], inds[1],
+                                       ind_coll);
+
+               for (int j=0; j<ind_coll.Size(); j++)
+               {
+                  int edge = ind_coll[j];
+                  int index = pncmesh.nodes[edge].edge_index;
+                  if (index >= 0 && !pncmesh.IsGhost(1, index))
+                  {
+                     sets_[t][s].insert(index);
+                  }
+               }
+            }
+            break;
+         case FACE:
+            for (int i=0; i<ni; i++)
+            {
+               pncmesh.pncent_sets->GetEntityIndex(t, s, i, inds);
+               BlockArray<int> ind_coll;
+               pncmesh.GetRefinedFaces(inds[0], inds[1], inds[2], inds[3],
+                                       ind_coll);
+
+               for (int j=0; j<ind_coll.Size(); j++)
+               {
+                  int face = ind_coll[j];
+                  int index = pncmesh.faces[face].index;
+                  if (index >= 0 && !pncmesh.IsGhost(2, index))
+                  {
+                     sets_[t][s].insert(index);
+                  }
+               }
+            }
+            break;
+
+         case ELEMENT:
+            for (int i=0; i<ni; i++)
+            {
+               int elem = (*pncmesh.pncent_sets)(t, s, i);
+               BlockArray<int> ind_coll;
+               pncmesh.GetRefinedElements(elem, ind_coll);
+
+               for (int j=0; j<ind_coll.Size(); j++)
+               {
+                  sets_[t][s].insert(pncmesh.elements[ind_coll[j]].index);
+               }
+            }
+            break;
+         default:
+            MFEM_ABORT("Unknown entity set type: \"" << GetTypeName(t) << "\"");
+      }
+      cout << MyRank_ << ": " << set_names_[t][s] << " " << s << " set size " <<
+           sets_[t][s].size() << "{";
+      for (set<int>::iterator it=sets_[t][s].begin(); it!=sets_[t][s].end(); it++)
+      {
+         cout << " " << *it;
+      }
+      cout << "}" << endl;
+   }
+   map<string,int>::iterator it;
+   cout << MyRank_ << ": set index by name ";
+   for (it=set_index_by_name_[t].begin(); it != set_index_by_name_[t].end(); it++)
+   {
+      cout << " " << it->first << "->" << it->second;
+   }
+   cout << endl;
+   cout << MyRank_ << ": done BuildEntitySets for type " << GetTypeName(t) << endl;
+}
+
+ParNCEntitySets::ParNCEntitySets(MPI_Comm comm, const NCMesh &ncmesh)
+   : NCEntitySets(*ncmesh.ncent_sets)
+{
+   MyComm_ = comm;
+   MPI_Comm_size(MyComm_, &NRanks_);
+   MPI_Comm_rank(MyComm_, &MyRank_);
+
+   if ( MyRank_ == 0 )
+   {
+      cout << "Entering ParNCEntitySets(NCMesh) c'tor" << endl;
+   }
+
+   if ( MyRank_ == 0 )
+   {
+      cout << "Leaving ParNCEntitySets(NCMesh) c'tor" << endl;
+   }
+}
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI

mesh/pentsets.hpp

@@ -0,0 +1,74 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_PAR_ENTITY_SETS
+#define MFEM_PAR_ENTITY_SETS
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "entsets.hpp"
+#include "../general/communication.hpp"
+
+namespace mfem
+{
+
+class ParMesh;
+class ParNCMesh;
+
+class ParEntitySets : public EntitySets
+{
+   friend class ParMesh;
+
+public:
+   ParEntitySets(const ParEntitySets & ent_sets);
+   ParEntitySets(ParMesh & _mesh, const EntitySets & ent_sets, int * part,
+                 const Array<int> & vert_global_local);
+   ParEntitySets(ParMesh & mesh, ParNCMesh &ncmesh);
+
+   virtual ~ParEntitySets();
+
+   virtual void PrintSetInfo(std::ostream &output) const;
+
+   inline ParMesh *GetParMesh() const { return pmesh_; }
+
+private:
+
+   void PrintEntitySetInfo(std::ostream & output, EntityType t,
+                           const std::string & ent_name) const;
+
+   void BuildEntitySets(ParNCMesh &pncmesh, EntityType t);
+
+   ParMesh  * pmesh_;
+   int        NRanks_;
+   int        MyRank_;
+};
+
+class ParNCEntitySets : public NCEntitySets
+{
+public:
+   // ParNCEntitySets(MPI_Comm comm, EntitySets &ent_sets, NCMesh &ncmesh);
+   ParNCEntitySets(MPI_Comm comm, const NCMesh &ncmesh);
+   // ParNCEntitySets(const ParMesh & pmesh, const ParNCMesh &pncmesh);
+   // ParNCEntitySets(const ParNCEntitySets & pncent_sets);
+
+private:
+   MPI_Comm MyComm_;
+   int      NRanks_;
+   int      MyRank_;
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_PAR_ENTITY_SETS

mesh/pmesh.cpp

@@ -91,6 +91,10 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
       *Nodes = *pmesh.Nodes;
       own_nodes = 1;
    }
+
+   // Copy entity sets if present in the input mesh
+   ent_sets = pent_sets =
+                 (pmesh.pent_sets) ? new ParEntitySets(*pmesh.pent_sets) : NULL;
 }
 
 ParMesh::ParMesh(ParMesh &&mesh) : ParMesh()
@@ -110,6 +114,7 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
    , glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
+   , pent_sets(NULL)
 {
    int *partitioning = NULL;
    Array<bool> activeBdrElem;
@@ -118,6 +123,8 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
    MPI_Comm_size(MyComm, &NRanks);
    MPI_Comm_rank(MyComm, &MyRank);
 
+   Array<int> vert_global_local;
+
    if (mesh.Nonconforming())
    {
       if (partitioning_)
@@ -148,6 +155,10 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
       mesh.bdr_attributes.Copy(bdr_attributes);
 
       GenerateNCFaceInfo();
+
+      //      if (mesh.ent_sets)
+      // NumOfVertices = BuildLocalVertices(mesh, partitioning,
+      //                                    vert_global_local);
    }
    else // mesh.Conforming()
    {
@@ -168,7 +179,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
       // re-enumerate the partitions to better map to actual processor
       // interconnect topology !?
 
-      Array<int> vert_global_local;
       NumOfVertices = BuildLocalVertices(mesh, partitioning, vert_global_local);
       NumOfElements = BuildLocalElements(mesh, partitioning, vert_global_local);
 
@@ -240,6 +250,12 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
 
       SetMeshGen();
       meshgen = mesh.meshgen; // copy the global 'meshgen'
+
+      ent_sets = pent_sets =
+                    (mesh.ent_sets) ? new ParEntitySets(*this, *mesh.ent_sets,
+                                                        partitioning,
+                                                        vert_global_local)
+                    : NULL;
    }
 
    if (mesh.NURBSext)
@@ -289,7 +305,12 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
       // for compatibility (e.g., Mesh::GetVertex())
       SetVerticesFromNodes(Nodes);
    }
-
+   /*
+   ent_sets = pent_sets =
+                 (mesh.ent_sets) ? new ParEntitySets(*this, *mesh.ent_sets,
+                                                     partitioning,
+                                                     vert_global_local) : NULL;
+   */
    if (partitioning != partitioning_)
    {
       delete [] partitioning;
@@ -859,6 +880,7 @@ ParMesh::ParMesh(const ParNCMesh &pncmesh)
    , glob_offset_sequence(-1)
    , gtopo(MyComm)
    , pncmesh(NULL)
+   , pent_sets(NULL)
 {
    Mesh::InitFromNCMesh(pncmesh);
    ReduceMeshGen();
@@ -925,6 +947,7 @@ ParMesh::ParMesh(MPI_Comm comm, istream &input, bool refine)
    , glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
+   , pent_sets(NULL)
 {
    MyComm = comm;
    MPI_Comm_size(MyComm, &NRanks);
@@ -1139,7 +1162,8 @@ void ParMesh::MakeRefined_(ParMesh &orig_mesh, int ref_factor, int ref_type)
    gtopo = orig_mesh.gtopo;
    have_face_nbr_data = false;
    pncmesh = NULL;
-
+   pent_sets = NULL;
+   
    Array<int> ref_factors(orig_mesh.GetNE());
    ref_factors = ref_factor;
    Mesh::MakeRefined_(orig_mesh, ref_factors, ref_type);
@@ -3768,6 +3792,13 @@ void ParMesh::NonconformingRefinement(const Array<Refinement> &refinements,
    // and this mesh will be the new fine mesh
    Mesh::Swap(*pmesh2, false);
 
+   // swap entity set information if present
+   mfem::Swap(pmesh2->pent_sets, this->pent_sets);
+   if (this->pent_sets)
+   {
+      this->pent_sets->pmesh_ = this;
+   }
+
    delete pmesh2; // NOTE: old face neighbors destroyed here
 
    pncmesh->GetConformingSharedStructures(*this);
@@ -6171,6 +6202,9 @@ void ParMesh::Destroy()
    delete pncmesh;
    ncmesh = pncmesh = NULL;
 
+   delete pent_sets;
+   ent_sets = pent_sets = NULL;
+
    DeleteFaceNbrData();
 
    for (int i = 0; i < shared_edges.Size(); i++)

mesh/pmesh.hpp

@@ -20,6 +20,7 @@
 #include "../general/globals.hpp"
 #include "mesh.hpp"
 #include "pncmesh.hpp"
+#include "pentsets.hpp"
 #include <iostream>
 
 namespace mfem
@@ -320,6 +321,7 @@ public:
    Table            send_face_nbr_vertices;
 
    ParNCMesh* pncmesh;
+   ParEntitySets* pent_sets;
 
    int GetNGroups() const { return gtopo.NGroups(); }
 

mesh/pncmesh.cpp

@@ -20,6 +20,8 @@
 #include <map>
 #include <climits> // INT_MIN, INT_MAX
 
+#include <fstream> // MLS Debugging
+
 namespace mfem
 {
 
@@ -27,6 +29,7 @@ using namespace bin_io;
 
 ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh, int *part)
    : NCMesh(ncmesh)
+   , pncent_sets(NULL)
 {
    MyComm = comm;
    MPI_Comm_size(MyComm, &NRanks);
@@ -41,6 +44,40 @@ ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh, int *part)
 
    Update();
 
+   std::ostringstream oss; oss << "elements_" << MyRank << ".out";
+   std::ofstream ofs(oss.str().c_str());
+
+   for (int i=0; i<elements.Size(); i++)
+   {
+      ofs << i
+          << '\t' << elements[i].index
+          << '\t' << elements[i].rank
+          << '\t' << elements[i].attribute
+          << '\t' << elements[i].parent;
+      if ( elements[i].ref_type == 0 )
+      {
+         ofs << " nodes {";
+         for (int j=0; j<8; j++)
+         {
+            ofs << " " << elements[i].node[j];
+         }
+         ofs << "}";
+      }
+      else
+      {
+         ofs << " children {";
+         for (int j=0; j<8; j++)
+         {
+            ofs << " " << elements[i].child[j];
+         }
+         ofs << "}";
+      }
+      ofs << std::endl;
+   }
+
+   ncent_sets = pncent_sets =
+                   (ncmesh.ncent_sets) ? new ParNCEntitySets(comm, ncmesh) : NULL;
+
    // note that at this point all processors still have all the leaf elements;
    // we however may now start pruning the refinement tree to get rid of
    // branches that only contain someone else's leaves (see Prune())
@@ -85,6 +122,9 @@ ParNCMesh::ParNCMesh(const ParNCMesh &other)
 ParNCMesh::~ParNCMesh()
 {
    ClearAuxPM();
+
+   delete pncent_sets;
+   ncent_sets = pncent_sets = NULL;
 }
 
 void ParNCMesh::Update()
@@ -115,6 +155,386 @@ void ParNCMesh::Update()
    boundary_layer.SetSize(0);
 }
 
+  /*
+void ParNCMesh::AssignLeafIndices()
+{
+   // This is an override of NCMesh::AssignLeafIndices(). The difference is
+   // that we shift all elements we own to the beginning of the array
+   // 'leaf_elements' and assign all ghost elements indices >= NElements.
+
+   // Also note that the ordering of ghosts and non-ghosts is preserved here,
+   // which is important for ParNCMesh::GetFaceNeighbors.
+
+   // We store the original leaf ordering in 'leaf_glob_order'. This is later
+   // used (and deleted) in GetConformingSharedStructures
+
+   NCMesh::AssignLeafIndices(); // original numbering, for 'leaf_glob_order'
+
+   int nleafs = leaf_elements.Size();
+
+   Array<int> ghosts;
+   ghosts.Reserve(nleafs);
+
+   NElements = 0;
+   for (int i = 0; i < nleafs; i++)
+   {
+      int elem = leaf_elements[i];
+      if (elements[elem].rank == MyRank)
+      {
+         leaf_elements[NElements++] = elem;
+      }
+      else
+      {
+         ghosts.Append(elem);
+      }
+   }
+   NGhostElements = ghosts.Size();
+
+   leaf_elements.SetSize(NElements);
+   leaf_elements.Append(ghosts);
+
+   // store original (globally consistent) numbering in 'leaf_glob_order'
+   leaf_glob_order.SetSize(nleafs);
+   for (int i = 0; i < nleafs; i++)
+   {
+      leaf_glob_order[i] = elements[leaf_elements[i]].index;
+   }
+
+   // new numbering with ghost shifted to the back
+   NCMesh::AssignLeafIndices();
+}
+
+void ParNCMesh::UpdateVertices()
+{
+   // This is an override of NCMesh::UpdateVertices. This version first
+   // assigns vert_index to vertices of elements of our rank. Only these
+   // vertices then make it to the Mesh in NCMesh::GetMeshComponents.
+   // The remaining (ghost) vertices are assigned indices greater or equal to
+   // Mesh::GetNV().
+
+   for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
+   {
+      if (node->HasVertex()) { node->vert_index = -1; }
+   }
+
+   NVertices = 0;
+   for (int i = 0; i < leaf_elements.Size(); i++)
+   {
+      Element &el = elements[leaf_elements[i]];
+      if (el.rank == MyRank)
+      {
+         for (int j = 0; j < GI[el.Geom()].nv; j++)
+         {
+            int &vindex = nodes[el.node[j]].vert_index;
+            if (vindex < 0) { vindex = NVertices++; }
+         }
+      }
+   }
+
+   vertex_nodeId.SetSize(NVertices);
+   for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
+   {
+      if (node->HasVertex() && node->vert_index >= 0)
+      {
+         vertex_nodeId[node->vert_index] = node.index();
+      }
+   }
+
+   NGhostVertices = 0;
+   for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
+   {
+      if (node->HasVertex() && node->vert_index < 0)
+      {
+         node->vert_index = NVertices + (NGhostVertices++);
+      }
+   }
+}
+
+void ParNCMesh::OnMeshUpdated(Mesh *mesh)
+{
+   std::cout << MyRank << ": Entering ParNCMesh::OnMeshUpdated" << std::endl;
+   // This is an override (or extension of) NCMesh::OnMeshUpdated().
+   // In addition to getting edge/face indices from 'mesh', we also
+   // assign indices to ghost edges/faces that don't exist in the 'mesh'.
+
+   // clear edge_index and Face::index
+   for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
+   {
+      if (node->HasEdge()) { node->edge_index = -1; }
+   }
+   for (face_iterator face = faces.begin(); face != faces.end(); ++face)
+   {
+      face->index = -1;
+   }
+
+   // go assign existing edge/face indices
+   NCMesh::OnMeshUpdated(mesh);
+
+   std::cout << MyRank << ": NVertices = " << NVertices << std::endl;
+
+   std::ostringstream ossN;
+   ossN << "node_on_mesh_updated_" << MyRank << ".out";
+   std::ofstream ofsN(ossN.str().c_str());
+   ofsN << nodes.Size() << std::endl;
+   for (int i=0; i<nodes.Size(); i++)
+   {
+      ofsN << i
+           // << " " << nodes[i].vert_refc
+           //  << " " << nodes[i].edge_refc
+           << " " << nodes[i].HasVertex()
+           << " " << nodes[i].HasEdge()
+           << " " << nodes[i].vert_index
+           << " " << nodes[i].edge_index
+           << " " << nodes[i].p1
+           << " " << nodes[i].p2
+           << " " << nodes[i].next << std::endl;
+   }
+   ofsN.close();
+
+   // count ghost edges and assign their indices
+   NEdges = mesh->GetNEdges();
+   NGhostEdges = 0;
+   for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
+   {
+      if (node->HasEdge() && node->edge_index < 0)
+      {
+         node->edge_index = NEdges + (NGhostEdges++);
+      }
+   }
+
+   // count ghost faces
+   NFaces = mesh->GetNumFaces();
+   NGhostFaces = 0;
+   for (face_iterator face = faces.begin(); face != faces.end(); ++face)
+   {
+      if (face->index < 0) { NGhostFaces++; }
+   }
+
+   if (Dim == 2)
+   {
+      // in 2D we have fake faces because of DG
+      MFEM_ASSERT(NFaces == NEdges, "");
+      MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
+   }
+
+   // resize face_geom (default_geom is for slave faces beyond the ghost layer)
+   Geometry::Type default_geom = Geometry::SQUARE;
+   face_geom.SetSize(NFaces + NGhostFaces, default_geom);
+
+   // update 'face_geom' for ghost faces, assign ghost face indices
+   int nghosts = 0;
+   for (int i = 0; i < NGhostElements; i++)
+   {
+      Element &el = elements[leaf_elements[NElements + i]]; // ghost element
+      GeomInfo &gi = GI[el.Geom()];
+
+      for (int j = 0; j < gi.nf; j++)
+      {
+         const int *fv = gi.faces[j];
+         Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
+                                 el.node[fv[2]], el.node[fv[3]]);
+         MFEM_ASSERT(face, "face not found!");
+
+         if (face->index < 0)
+         {
+            face->index = NFaces + (nghosts++);
+
+            // store the face geometry
+            static const Geometry::Type types[5] =
+            {
+               Geometry::INVALID, Geometry::INVALID,
+               Geometry::SEGMENT, Geometry::TRIANGLE, Geometry::SQUARE
+            };
+            face_geom[face->index] = types[gi.nfv[j]];
+         }
+      }
+   }
+
+   // assign valid indices also to faces beyond the ghost layer
+   for (face_iterator face = faces.begin(); face != faces.end(); ++face)
+   {
+      if (face->index < 0) { face->index = NFaces + (nghosts++); }
+   }
+   MFEM_ASSERT(nghosts == NGhostFaces, "");
+
+   {
+      /// Debugging output
+      std::ostringstream oss; oss << "elements_on_mesh_updated_"
+                                  << MyRank << ".out";
+      std::ofstream ofs(oss.str().c_str());
+
+      for (int i=0; i<elements.Size(); i++)
+      {
+         ofs << i
+             << '\t' << elements[i].index
+             << '\t' << elements[i].rank
+             << '\t' << elements[i].attribute
+             << '\t' << elements[i].parent;
+         if ( elements[i].ref_type == 0 )
+         {
+            ofs << " nodes {";
+            for (int j=0; j<8; j++)
+            {
+               ofs << " " << elements[i].node[j];
+            }
+            ofs << "}";
+         }
+         else
+         {
+            ofs << " children {";
+            for (int j=0; j<8; j++)
+            {
+               ofs << " " << elements[i].child[j];
+            }
+            ofs << "}";
+         }
+         ofs << std::endl;
+      }
+
+      if (pncent_sets)
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated pncent_sets is non NULL" << std::endl;
+      }
+      else
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated pncent_sets is NULL" << std::endl;
+      }
+      if (ncent_sets)
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated ncent_sets is non NULL" << std::endl;
+      }
+      else
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated ncent_sets is NULL" << std::endl;
+      }
+      if (mesh->ent_sets)
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated mesh->ent_sets is non NULL" << std::endl;
+      }
+      else
+      {
+         std::cout << "ParNCMesh::OnMeshUpdated mesh->ent_sets is NULL" << std::endl;
+      }
+      ParMesh * pmesh = dynamic_cast<ParMesh*>(mesh);
+      if (pmesh)
+      {
+         std::cout << "dynamic cast succeeded: mesh is a ParMesh" << std::endl;
+
+         if (pmesh->pent_sets != NULL)
+         {
+            std::cout << "ParNCMesh::OnMeshUpdated deleting ParEntitySets object in ParMesh"
+                      << std::endl;
+            delete pmesh->pent_sets;
+         }
+         else if (pmesh->ent_sets != NULL)
+         {
+            std::cout << "ParNCMesh::OnMeshUpdated deleting EntitySets object in ParMesh" <<
+                      std::endl;
+            delete pmesh->ent_sets;
+         }
+         std::cout << "ParNCMesh::OnMeshUpdated creating ParEntitySets object in ParMesh"
+                   << std::endl;
+         pmesh->ent_sets = pmesh->pent_sets =
+                              (pncent_sets) ? new ParEntitySets(*pmesh, *this): NULL;
+	 */
+         /*
+         if (pmesh->ent_sets)
+         {
+         std::cout << MyRank << ": ParNCMesh::OnMeshUpdated pmesh->ent_sets is non NULL" << std::endl;
+         pmesh->ent_sets->PrintSetInfo(std::cout);
+
+         std::ostringstream oss; oss << "ent_sets_" << MyRank << ".out";
+         std::ofstream ofs(oss.str().c_str());
+         pmesh->ent_sets->Print(ofs);
+         MPI_Barrier(MyComm);
+
+         std::cout << MyRank << ": testing " << NElements << std::endl;
+         //pmesh->ent_sets->Prune(NElements);
+         }
+         else
+         {
+         std::cout << "ParNCMesh::OnMeshUpdated pmesh->ent_sets is NULL" << std::endl;
+         }
+         */
+	 /*
+         if (pmesh->pent_sets)
+         {
+            std::cout << "ParNCMesh::OnMeshUpdated pmesh->pent_sets is non NULL" <<
+                      std::endl;
+         }
+         else
+         {
+            std::cout << "ParNCMesh::OnMeshUpdated pmesh->pent_sets is NULL" << std::endl;
+         }
+      }
+      else
+      {
+         std::cout << "dynamic cast failed: mesh is not a ParMesh" << std::endl;
+      }
+      */
+      /*
+      if (pncent_sets)
+      {
+        if (!pmesh->pent_sets)
+        {
+          pmesh->pent_sets = new ParEntitySets(*pmesh, *this);
+        }
+      }
+      */
+      /*
+      // Prune the Entity Sets
+      if ( entity_sets )
+      {
+         EntitySets::EntityType t = EntitySets::INVALID;
+         unsigned int ns = -1;
+
+         std::cout << "Processing node sets" << std::endl;
+
+         t = EntitySets::VERTEX;
+         ns = entity_sets->GetNumSets(t);
+         for (unsigned int s=0; s<ns; s++)
+         {
+            unsigned int ni = entity_sets->GetNumEntities(t, s);
+       int e = 0;
+       for (unsigned int i=0; i<ni; i++)
+         {
+           if ( (*mesh->ent_sets)(t, s, i) < NVertices )
+             {
+          (*mesh->ent_sets)(t, s, e) = (*mesh->ent_sets)(t, s, i);
+          e++;
+             }
+         }
+       (*mesh->ent_sets)(t, s).resize(e);
+         }
+
+         t = EntitySets::EDGE;
+         ns = entity_sets->GetNumSets(t);
+         for (unsigned int s=0; s<ns; s++)
+         {
+            unsigned int ni = entity_sets->GetNumEntities(t, s);
+       BlockArray<int> ids;
+
+            for (unsigned int i=0; i<ni; i++)
+            {
+         if ( (*mesh->ent_sets)(t, s, i) < NEdges )
+         {
+           ids.Append((*mesh->ent_sets)(t, s, i));
+         }
+            }
+       (*mesh->ent_sets)(t, s).resize(ids.Size());
+       for (int i=0; i<ids.Size(); i++)
+         {
+           (*mesh->ent_sets)(t, s, i) = ids[i];
+         }
+         }
+      }
+      */
+	/*
+      std::cout << MyRank << ": Leaving ParNCMesh::OnMeshUpdated" << std::endl;
+   }
+}
+*/
+
 void ParNCMesh::ElementSharesFace(int elem, int local, int face)
 {
    // Analogous to ElementSharesEdge.
@@ -2732,6 +3152,94 @@ void ParNCMesh::GetDebugMesh(Mesh &debug_mesh) const
    debug_mesh.ncmesh = copy;
 }
 
+void ParNCMesh::GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edges)
+{
+   std::cout << MyRank
+             << ": entering ParNCMesh::GetRefinedEdges "
+             <<"searching for edge with vertices: " << vn0 << " and " << vn1
+             << std::endl;
+   return this->NCMesh::GetRefinedEdges(vn0, vn1, edges);
+
+   int mid = nodes.FindId(vn0, vn1);
+   if (mid < 0) { return; }
+
+   /*
+   Node &nd = nodes[mid];
+
+   if ( nd.edge_index < 0 ) { return; }
+
+   edges.Append(nd.edge_index);
+
+   GetRefinedEdges(vn0, mid, edges);
+   GetRefinedEdges(mid, vn1, edges);
+   */
+   edges.Append(mid);
+
+   GetRefinedEdges(vn0, mid, edges);
+   GetRefinedEdges(mid, vn1, edges);
+}
+
+void ParNCMesh::GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
+                                BlockArray<int> & face_ids)
+{
+   return this->NCMesh::GetRefinedFaces(vn0, vn1, vn2, vn3, face_ids);
+   /*
+   Face* fa = faces.Find(vn0, vn1, vn2, vn3);
+
+   if (fa)
+   {
+      if ( fa->index >= 0 )
+      {
+         face_ids.Append(fa->index);
+      }
+      return;
+   }
+
+   // we need to recurse deeper
+   int mid[4];
+   int split = FaceSplitType(vn0, vn1, vn2, vn3, mid);
+
+   if (split == 1) // "X" split face
+   {
+      GetRefinedFaces(vn0, mid[0], mid[2], vn3, face_ids);
+      GetRefinedFaces(mid[0], vn1, vn2, mid[2], face_ids);
+   }
+   else if (split == 2) // "Y" split face
+   {
+      GetRefinedFaces(vn0, vn1, mid[1], mid[3], face_ids);
+      GetRefinedFaces(mid[3], mid[1], vn2, vn3, face_ids);
+   }
+   */
+}
+
+void ParNCMesh::GetRefinedElements(int elem_id, BlockArray<int> & elem_ids)
+{
+   // std::cout << MyRank
+   //           << ": entering ParNCMesh::GetRefinedElements "
+   //           <<"searching for element id: " << elem_id << std::endl;
+   Element &el = elements[elem_id];
+
+   if (el.ref_type != 0)
+   {
+      // This element has been refined so recurse into its children
+      for (int i = 0; i < 8; i++)
+      {
+         if (el.child[i] >= 0 && el.child[i] < elements.Size() )
+         {
+            GetRefinedElements(el.child[i], elem_ids);
+         }
+      }
+   }
+   else
+   {
+      // This element has not been refined so add it if it's a local element
+      if (el.rank == MyRank)
+      {
+         elem_ids.Append(elem_id);
+      }
+   }
+}
+
 void ParNCMesh::Trim()
 {
    NCMesh::Trim();

mesh/pncmesh.hpp

@@ -20,6 +20,7 @@
 #include <set>
 
 #include "ncmesh.hpp"
+#include "pentsets.hpp"
 #include "../general/communication.hpp"
 #include "../general/sort_pairs.hpp"
 
@@ -248,9 +249,29 @@ public:
        The debug mesh will have element attributes set to element rank + 1. */
    void GetDebugMesh(Mesh &debug_mesh) const;
 
+   /** Collect edge indices of all refined edges which are children of
+       the coarse edge defined by the given vertices.  This method
+       overrides a method in NCMesh and only returns locally owned
+       edges. */
+   void GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edge_ids);
+
+   /** Collect face indices of all refined faces which are children of
+       the coarse face defined by the given vertices.  This method
+       overrides a method in NCMesh and only returns locally owned
+       faces. */
+   void GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
+                        BlockArray<int> & face_ids);
+
+   /** Collect element indices of all refined elements which are
+       children of the coarse element defined by the given element
+       index. This method overrides a method in NCMesh and only
+       returns locally owned elements. */
+   void GetRefinedElements(int elem_id, BlockArray<int> & elem_ids);
+
 protected: // interface for ParMesh
 
    friend class ParMesh;
+   friend class ParEntitySets;
 
    /** For compatibility with conforming code in ParMesh and ParFESpace.
        Initializes shared structures in ParMesh: gtopo, shared_*, group_s*, s*_l*.
@@ -540,6 +561,8 @@ protected: // implementation
    Array<DenseMatrix*> aux_pm_store;
    void ClearAuxPM();
 
+   ParNCEntitySets * pncent_sets;
+
    long GroupsMemoryUsage() const;
 
    friend class NeighborRowMessage;