Parallel P construction now works with one process.

examples/CMakeLists.txt

@@ -11,6 +11,7 @@
 
 list(APPEND ALL_EXE_SRCS
   ex1.cpp
+  ex1-subdomain.cpp
   ex2.cpp
   ex3.cpp
   ex4.cpp
@@ -39,6 +40,7 @@ list(APPEND ALL_EXE_SRCS
 if (MFEM_USE_MPI)
   list(APPEND ALL_EXE_SRCS
     ex1p.cpp
+    ex1p-subdomain.cpp
     ex2p.cpp
     ex3p.cpp
     ex4p.cpp

examples/ex1-subdomain.cpp

@@ -0,0 +1,283 @@
+//                                MFEM Example 1
+//
+// Compile with: make ex1
+//
+// Sample runs:  ex1 -m ../data/square-disc.mesh
+//               ex1 -m ../data/star.mesh
+//               ex1 -m ../data/star-mixed.mesh
+//               ex1 -m ../data/escher.mesh
+//               ex1 -m ../data/fichera.mesh
+//               ex1 -m ../data/fichera-mixed.mesh
+//               ex1 -m ../data/toroid-wedge.mesh
+//               ex1 -m ../data/periodic-annulus-sector.msh
+//               ex1 -m ../data/periodic-torus-sector.msh
+//               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/star-mixed-p2.mesh -o 2
+//               ex1 -m ../data/disc-nurbs.mesh -o -1
+//               ex1 -m ../data/pipe-nurbs.mesh -o -1
+//               ex1 -m ../data/fichera-mixed-p2.mesh -o 2
+//               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
+//
+// Device sample runs:
+//               ex1 -pa -d cuda
+//               ex1 -pa -d raja-cuda
+//             * ex1 -pa -d raja-hip
+//               ex1 -pa -d occa-cuda
+//               ex1 -pa -d raja-omp
+//               ex1 -pa -d occa-omp
+//               ex1 -pa -d ceed-cpu
+//             * ex1 -pa -d ceed-cuda
+//             * ex1 -pa -d ceed-hip
+//               ex1 -pa -d ceed-cuda:/gpu/cuda/shared
+//               ex1 -m ../data/beam-hex.mesh -pa -d cuda
+//               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cuda:/gpu/cuda/ref
+//
+// 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 = "../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   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) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   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. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. 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(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. 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 =
+         (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 5. 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;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (mesh.GetNodes())
+   {
+      fec = mesh.GetNodes()->OwnFEC();
+      delete_fec = false;
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+
+   Array<int> subdomain_attr(1);
+   subdomain_attr = 1;
+   for(int e=0;e<mesh.GetNE();e++)
+   {
+      Vector center(mesh.SpaceDimension());
+      mesh.GetElementCenter(e, center);
+
+      if(center(0) > 0.0)
+      {
+         mesh.SetAttribute(e, 2);
+      }
+   }
+   auto subdomain = SubdomainFromAttributes(&mesh, subdomain_attr);
+
+   FiniteElementSpace fespace(&mesh, fec, 1, Ordering::byNODES, subdomain);
+   std::cout << "Subdomain with " << fespace.GetNE() << "/" << mesh.GetNE() << std::endl;
+   cout << "Number of finite element unknowns: "
+        << fespace.GetTrueVSize() << endl;
+
+   // 6. 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 (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 7. 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(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 8. 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;
+
+   // 9. 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(&fespace);
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 10. 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();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   cout << "Size of linear system: " << A->Height() << endl;
+
+   // 11. Solve the linear system A X = B.
+   if (!pa)
+   {
+#ifndef MFEM_USE_SUITESPARSE
+      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+#else
+      // 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
+   }
+   else // Jacobi preconditioning in partial assembly mode
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
+      }
+      else
+      {
+         CG(*A, B, X, 1, 400, 1e-12, 0.0);
+      }
+   }
+
+   // 12. Recover the solution as a finite element grid function.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 13. 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);
+
+   // 14. 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);
+
+      // Transfer from subspace to full space for visualization.
+      FiniteElementSpace fespace_full(&mesh, fec, 1, Ordering::byNODES);
+      GridFunction x_full(&fespace_full);
+      x_full = 0.0;
+      Array<int> dofs_full, dofs;
+      Vector elemvect;
+      for(int ei = 0; ei < fespace.GetNE(); ei++)
+      {
+         fespace.GetElementDofs(ei, dofs);
+         x.GetSubVector(dofs, elemvect);
+         fespace_full.GetElementDofs(fespace.MapElement(ei), dofs_full);
+         x_full.SetSubVector(dofs_full, elemvect);
+      }
+
+      sol_sock << "solution\n" << mesh << x_full << flush;
+   }
+
+   // 15. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+
+   return 0;
+}

fem/bilinearform.cpp

@@ -460,7 +460,7 @@ void BilinearForm::Assemble(int skip_zeros)
 
       for (int i = 0; i < fes -> GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
          const FiniteElement &be = *fes->GetBE(i);
@@ -504,10 +504,10 @@ void BilinearForm::Assemble(int skip_zeros)
       FaceElementTransformations *tr;
       Array<int> vdofs2;
 
-      int nfaces = mesh->GetNumFaces();
+      int nfaces = fes->GetNF();
       for (int i = 0; i < nfaces; i++)
       {
-         tr = mesh -> GetInteriorFaceTransformations (i);
+         tr = fes -> GetInteriorFaceTransformations (i);
          if (tr != NULL)
          {
             fes -> GetElementVDofs (tr -> Elem1No, vdofs);
@@ -552,10 +552,10 @@ void BilinearForm::Assemble(int skip_zeros)
 
       for (int i = 0; i < fes -> GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
-         tr = mesh -> GetBdrFaceTransformations (i);
+         tr = fes -> GetBdrFaceTransformations (i);
          if (tr != NULL)
          {
             fes -> GetElementVDofs (tr -> Elem1No, vdofs);
@@ -1111,6 +1111,9 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
    extern_bfs = 0;
    assembly = AssemblyLevel::LEGACYFULL;
    ext = NULL;
+
+   MFEM_ASSERT(tr_fes->Subdomain() == te_fes->Subdomain(),
+   "Test and trial space must live on the same domain.");
 }
 
 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
@@ -1347,7 +1350,7 @@ void MixedBilinearForm::Assemble (int skip_zeros)
 
       for (int i = 0; i < test_fes -> GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = trial_fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
          trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
@@ -1372,10 +1375,10 @@ void MixedBilinearForm::Assemble (int skip_zeros)
       Array<int> te_vdofs2;
       const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
 
-      int nfaces = mesh->GetNumFaces();
+      int nfaces = trial_fes->GetNF();
       for (int i = 0; i < nfaces; i++)
       {
-         ftr = mesh->GetFaceElementTransformations(i);
+         ftr = trial_fes->GetFaceElementTransformations(i);
          trial_fes->GetFaceVDofs(i, tr_vdofs);
          test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
          trial_face_fe = trial_fes->GetFaceElement(i);
@@ -1431,10 +1434,10 @@ void MixedBilinearForm::Assemble (int skip_zeros)
 
       for (int i = 0; i < trial_fes -> GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = trial_fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
-         ftr = mesh->GetBdrFaceTransformations(i);
+         ftr = trial_fes->GetBdrFaceTransformations(i);
          if (ftr)
          {
             trial_fes->GetFaceVDofs(ftr->ElementNo, tr_vdofs);
@@ -1837,12 +1840,12 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
 
    if (tfbfi.Size())
    {
-      const int nfaces = test_fes->GetMesh()->GetNumFaces();
+      const int nfaces = test_fes->GetNF();
       for (int i = 0; i < nfaces; i++)
       {
          trial_fes->GetFaceVDofs(i, dom_vdofs);
          test_fes->GetFaceVDofs(i, ran_vdofs);
-         T = test_fes->GetMesh()->GetFaceTransformation(i);
+         T = test_fes->GetFaceTransformation(i);
          dom_fe = trial_fes->GetFaceElement(i);
          ran_fe = test_fes->GetFaceElement(i);
 

fem/bilinearform_ext.cpp

@@ -505,7 +505,7 @@ void EABilinearFormExtension::Assemble()
 {
    SetupRestrictionOperators(L2FaceValues::SingleValued);
 
-   ne = trialFes->GetMesh()->GetNE();
+   ne = trialFes->GetNE();
    elemDofs = trialFes->GetFE(0)->GetDof();
 
    ea_data.SetSize(ne*elemDofs*elemDofs, Device::GetMemoryType());

fem/fespace.cpp

@@ -74,6 +74,7 @@ FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
 {
    mesh = mesh ? mesh : orig.mesh;
    fec = fec ? fec : orig.fec;
+   subdomain = orig.subdomain;
    NURBSExtension *NURBSext = NULL;
    if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
    {
@@ -95,13 +96,13 @@ FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
 
 int FiniteElementSpace::GetOrder(int i) const
 {
-   Geometry::Type GeomType = mesh->GetElementBaseGeometry(i);
+   Geometry::Type GeomType = mesh->GetElementBaseGeometry(MapElement(i));
    return fec->FiniteElementForGeometry(GeomType)->GetOrder();
 }
 
 int FiniteElementSpace::GetFaceOrder(int i) const
 {
-   Geometry::Type GeomType = mesh->GetFaceBaseGeometry(i);
+   Geometry::Type GeomType = mesh->GetFaceBaseGeometry(MapElement(i));
    return fec->FiniteElementForGeometry(GeomType)->GetOrder();
 }
 
@@ -218,14 +219,14 @@ void FiniteElementSpace::BuildElementToDofTable() const
 
    Table *el_dof = new Table;
    Array<int> dofs;
-   el_dof -> MakeI (mesh -> GetNE());
-   for (int i = 0; i < mesh -> GetNE(); i++)
+   el_dof -> MakeI (this -> GetNE());
+   for (int i = 0; i < this -> GetNE(); i++)
    {
       GetElementDofs (i, dofs);
       el_dof -> AddColumnsInRow (i, dofs.Size());
    }
    el_dof -> MakeJ();
-   for (int i = 0; i < mesh -> GetNE(); i++)
+   for (int i = 0; i < this -> GetNE(); i++)
    {
       GetElementDofs (i, dofs);
       el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
@@ -240,14 +241,14 @@ void FiniteElementSpace::BuildBdrElementToDofTable() const
 
    Table *bel_dof = new Table;
    Array<int> dofs;
-   bel_dof->MakeI(mesh->GetNBE());
-   for (int i = 0; i < mesh->GetNBE(); i++)
+   bel_dof->MakeI(this->GetNBE());
+   for (int i = 0; i < this->GetNBE(); i++)
    {
       GetBdrElementDofs(i, dofs);
       bel_dof->AddColumnsInRow(i, dofs.Size());
    }
    bel_dof->MakeJ();
-   for (int i = 0; i < mesh->GetNBE(); i++)
+   for (int i = 0; i < this->GetNBE(); i++)
    {
       GetBdrElementDofs(i, dofs);
       bel_dof->AddConnections(i, (int *)dofs, dofs.Size());
@@ -266,7 +267,7 @@ void FiniteElementSpace::BuildFaceToDofTable() const
 
    Table *fc_dof = new Table;
    Array<int> dofs;
-   fc_dof->MakeI(mesh->GetNumFaces());
+   fc_dof->MakeI(this->GetNF());
    for (int i = 0; i < fc_dof->Size(); i++)
    {
       GetFaceDofs(i, dofs);
@@ -318,7 +319,7 @@ void FiniteElementSpace::BuildDofToArrays()
    dof_elem_array.SetSize (ndofs);
    dof_ldof_array.SetSize (ndofs);
    dof_elem_array = -1;
-   for (int i = 0; i < mesh -> GetNE(); i++)
+   for (int i = 0; i < this -> GetNE(); i++)
    {
       const int *dofs = elem_dof -> GetRow(i);
       const int n = elem_dof -> RowSize(i);
@@ -352,7 +353,7 @@ void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
    ess_vdofs.SetSize(GetVSize());
    ess_vdofs = 0;
 
-   for (int i = 0; i < GetNBE(); i++)
+   for (int i = 0; i < this->GetNBE(); i++)
    {
       if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
       {
@@ -376,39 +377,137 @@ void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
    // local DOFs affected by boundary elements on other processors
    if (Nonconforming())
    {
-      Array<int> bdr_verts, bdr_edges;
-      mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
+      Array<int> bdr_vertices, bdr_edges;
+      mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_vertices, bdr_edges);
 
-      for (int i = 0; i < bdr_verts.Size(); i++)
+      int ignored_vertices = 0;
+      int ignored_edges = 0;
+
+      for (const int vertex_in_mesh : bdr_vertices)
       {
+         const int vertex = MapVertexBack(vertex_in_mesh);
+         if(vertex<0)
+         {
+            ignored_vertices++;
+            continue;
+         }
+
          if (component < 0)
          {
-            GetVertexVDofs(bdr_verts[i], vdofs);
+            GetVertexVDofs(vertex, vdofs);
             mark_dofs(vdofs, ess_vdofs);
          }
          else
          {
-            GetVertexDofs(bdr_verts[i], dofs);
+            GetVertexDofs(vertex, 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 < bdr_edges.Size(); i++)
+
+      // Some debug info.
+      if(ignored_vertices>0)
       {
+         MFEM_WARNING(ignored_vertices << "/" << bdr_vertices.Size() << " closure vertices ignored. Are your boundaries okay?" << std::endl);
+#ifdef MFEM_DEBUG
+         auto v2e = mesh->GetVertexToElementTable();
+         Array<int> elements;
+         Array<int> ignored_vertices;
+         for (const int vertex_in_mesh : bdr_vertices)
+         {
+            // These are fine.
+            if(MapVertexBack(vertex_in_mesh)>= 0)
+            {
+               continue;
+            }
+
+            v2e->GetRow(vertex_in_mesh, elements);
+
+            bool has_element_in_subspace = false;
+            for(const int e : elements)
+            {
+               // Dirty hack to check if an element is part of our subspace
+               if(MapElementBack(e) >= 0)
+               {
+                  has_element_in_subspace = true;
+               }
+            }
+
+            if(has_element_in_subspace)
+            {
+               ignored_vertices.Append(vertex_in_mesh);
+            }
+         }
+         delete v2e;
+
+         MFEM_ASSERT(ignored_vertices.Size() == 0, "Something went wrong when applying DBCs. " << ignored_vertices.Size() << " vertices were ignored." << ignored_vertices);
+#endif
+      }
+
+      for (const int edge_in_mesh : bdr_edges)
+      {
+         const int edge = MapEdgeBack(edge_in_mesh);
+         if(edge<0)
+         {
+            ignored_edges++;
+            continue;
+         }
+
          if (component < 0)
          {
-            GetEdgeVDofs(bdr_edges[i], vdofs);
+            GetEdgeVDofs(edge, vdofs);
             mark_dofs(vdofs, ess_vdofs);
          }
          else
          {
-            GetEdgeDofs(bdr_edges[i], dofs);
-            for (int d = 0; d < dofs.Size(); d++)
-            { dofs[d] = DofToVDof(dofs[d], component); }
+            GetEdgeDofs(edge, dofs);
+            for (int& dof : dofs)
+            {
+               dof = DofToVDof(dof, component);
+            }
             mark_dofs(dofs, ess_vdofs);
          }
       }
+      // Some debug info.
+      if(ignored_edges>0)
+      {
+         MFEM_WARNING(ignored_edges << "/" << bdr_edges.Size() << " closure edges ignored. Are your boundaries okay?" << std::endl);
+#ifdef MFEM_DEBUG
+         auto edge2element = mesh->GetEdgeToElementTable();
+         Array<int> elements;
+         Array<int> ignored_edges;
+         for (const int edge_in_mesh : bdr_edges)
+         {
+            // These are fine.
+            if(MapEdgeBack(edge_in_mesh)>= 0)
+            {
+               continue;
+            }
+
+            bool has_element_in_subspace = false;
+
+            edge2element->GetRow(edge_in_mesh, elements);
+
+            for(const int e : elements)
+            {
+               // Dirty hack to check if an element is part of our subspace
+               if(MapElementBack(e) >= 0)
+               {
+                  has_element_in_subspace = true;
+               }
+            }
+
+            if(has_element_in_subspace)
+            {
+               ignored_edges.Append(edge_in_mesh);
+            }
+         }
+         delete edge2element;
+
+         MFEM_ASSERT(ignored_edges.Size() == 0, "Something went wrong when applying DBCs. " << ignored_edges.Size() << " edges were ignored." << ignored_edges);
+#endif
+      }
    }
 }
 
@@ -488,7 +587,7 @@ FiniteElementSpace::D2C_GlobalRestrictionMatrix (FiniteElementSpace *cfes)
 
    R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
 
-   for (i = 0; i < mesh -> GetNE(); i++)
+   for (i = 0; i < this -> GetNE(); i++)
    {
       this -> GetElementVDofs (i, d_vdofs);
       cfes -> GetElementVDofs (i, c_vdofs);
@@ -520,7 +619,7 @@ FiniteElementSpace::D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
 
    R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
 
-   for (i = 0; i < mesh -> GetNE(); i++)
+   for (i = 0; i < this -> GetNE(); i++)
    {
       this -> GetElementDofs (i, d_dofs);
       cfes -> GetElementDofs (i, c_dofs);
@@ -557,13 +656,13 @@ FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
    const FiniteElement *l_fe = NULL;
    IsoparametricTransformation T;
 
-   for (int i = 0; i < mesh -> GetNE(); i++)
+   for (int i = 0; i < this -> GetNE(); i++)
    {
       this -> GetElementDofs (i, h_dofs);
       lfes -> GetElementDofs (i, l_dofs);
 
       // Assuming 'loc_restr' depends only on the Geometry::Type.
-      const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
+      const Geometry::Type geom = mesh->GetElementBaseGeometry(MapElement(i));
       if (geom != cached_geom)
       {
          h_fe = this -> GetFE (i);
@@ -711,7 +810,17 @@ void FiniteElementSpace::BuildConformingInterpolation() const
       {
          const NCMesh::Master &master = list.masters[mi];
 
-         GetEntityDofs(entity, master.index, master_dofs);
+         // Skip DoF outside of subdomain
+         if(MapElementBack(mesh->ncmesh->GetLeafElement(master.element)) == -1)
+         {
+            continue;
+         }
+
+         // FIXME: This does not work for NC tets.
+         const int master_index = MapEntityBack(entity, master.index);
+         if(master_index < 0) continue;
+
+         GetEntityDofs(entity, master_index, master_dofs);
          if (!master_dofs.Size()) { continue; }
 
          const FiniteElement* fe = fec->FiniteElementForGeometry(master.Geom());
@@ -728,12 +837,22 @@ void FiniteElementSpace::BuildConformingInterpolation() const
          for (int si = master.slaves_begin; si < master.slaves_end; si++)
          {
             const NCMesh::Slave &slave = list.slaves[si];
-            GetEntityDofs(entity, slave.index, slave_dofs, master.Geom());
+
+            // Skip DoF outside of subdomain
+            if(MapElementBack(mesh->ncmesh->GetLeafElement(slave.element)) == -1)
+            {
+               continue;
+            }
+
+            // FIXME: This does not work for NC tets.
+            const int slave_index = MapEntityBack(entity, slave.index);
+            if(slave_index < 0) continue;
+
+            GetEntityDofs(entity, slave_index, slave_dofs, master.Geom());
             if (!slave_dofs.Size()) { continue; }
 
             list.OrientedPointMatrix(slave, T.GetPointMat());
             fe->GetLocalInterpolation(T, I);
-
             // make each slave DOF dependent on all master DOFs
             AddDependencies(deps, master_dofs, slave_dofs, I);
          }
@@ -1034,10 +1153,10 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
 
    const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
 
-   for (int k = 0; k < mesh->GetNE(); k++)
+   for (int k = 0; k < this->GetNE(); k++)
    {
-      const Embedding &emb = rtrans.embeddings[k];
-      const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
+      const Embedding &emb = rtrans.embeddings[MapElement(k)];
+      const Geometry::Type geom = mesh->GetElementBaseGeometry(MapElement(k));
       const DenseMatrix &lP = localP[geom](emb.matrix);
       const int fine_ldof = localP[geom].SizeI();
 
@@ -1163,10 +1282,11 @@ void FiniteElementSpace::RefinementOperator
 
    Vector subY, subX;
 
-   for (int k = 0; k < mesh->GetNE(); k++)
+   for (int k = 0; k < fespace->GetNE(); k++)
    {
-      const Embedding &emb = rtrans.embeddings[k];
-      const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
+      const auto k_mesh = fespace->MapElement(k);
+      const Embedding &emb = rtrans.embeddings[k_mesh];
+      const Geometry::Type geom = mesh->GetElementBaseGeometry(k_mesh);
       const DenseMatrix &lP = localP[geom](emb.matrix);
 
       subY.SetSize(lP.Height());
@@ -1205,10 +1325,11 @@ void FiniteElementSpace::RefinementOperator
 
    Vector subY, subX;
 
-   for (int k = 0; k < mesh->GetNE(); k++)
+   for (int k = 0; k < fespace->GetNE(); k++)
    {
-      const Embedding &emb = rtrans.embeddings[k];
-      const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
+      const auto k_mesh = fespace->MapElement(k);
+      const Embedding &emb = rtrans.embeddings[k_mesh];
+      const Geometry::Type geom = mesh->GetElementBaseGeometry(k_mesh);
       const DenseMatrix &lP = localP[geom](emb.matrix);
 
       fespace->GetElementDofs(k, f_dofs);
@@ -1432,6 +1553,7 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
    MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
 
    int num_marked = 0;
+   // TODO (fe.subdomain): this needs some investigation.
    for (int k = 0; k < dtrans.embeddings.Size(); k++)
    {
       const Embedding &emb = dtrans.embeddings[k];
@@ -1589,7 +1711,7 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
    face_to_be = -1;
    for (int b = 0; b < GetNBE(); b++)
    {
-      int f = mesh->GetBdrElementEdgeIndex(b);
+      int f = mesh->GetBdrElementEdgeIndex(MapBoundary(b));
       face_to_be[f] = b;
    }
 
@@ -1605,8 +1727,8 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
       //        same orientation.
       if (dim > 1)
       {
-         const Element *fe = mesh->GetFace(f);
-         const Element *be = mesh->GetBdrElement(b);
+         const Element *fe = mesh->GetFace(MapFace(f));
+         const Element *be = mesh->GetBdrElement(MapBoundary(b));
          const int nv = be->GetNVertices();
          const int *fv = fe->GetVertices();
          const int *bv = be->GetVertices();
@@ -1645,14 +1767,14 @@ void FiniteElementSpace::Construct()
    cP_is_set = false;
    // 'Th' is initialized/destroyed before this method is called.
 
-   nvdofs = mesh->GetNV() * fec->DofForGeometry(Geometry::POINT);
+   nvdofs = this->GetNV() * fec->DofForGeometry(Geometry::POINT);
 
    if (mesh->Dimension() > 1)
    {
-      nedofs = mesh->GetNEdges() * fec->DofForGeometry(Geometry::SEGMENT);
+      nedofs = this->GetNEdges() * fec->DofForGeometry(Geometry::SEGMENT);
    }
 
-   if (mesh->GetNFaces() > 0)
+   if (this->GetNFaces() > 0)
    {
       bool have_face_dofs = false;
       for (int g = Geometry::DimStart[2]; g < Geometry::DimStart[3]; g++)
@@ -1666,11 +1788,11 @@ void FiniteElementSpace::Construct()
       }
       if (have_face_dofs)
       {
-         fdofs = new int[mesh->GetNFaces()+1];
+         fdofs = new int[this->GetNFaces()+1];
          fdofs[0] = 0;
-         for (int i = 0; i < mesh->GetNFaces(); i++)
+         for (int i = 0; i < this->GetNFaces(); i++)
          {
-            nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
+            nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(MapFace(i)));
             fdofs[i+1] = nfdofs;
          }
       }
@@ -1678,11 +1800,11 @@ void FiniteElementSpace::Construct()
 
    if (mesh->Dimension() > 0)
    {
-      bdofs = new int[mesh->GetNE()+1];
+      bdofs = new int[this->GetNE()+1];
       bdofs[0] = 0;
-      for (int i = 0; i < mesh->GetNE(); i++)
+      for (int i = 0; i < this->GetNE(); i++)
       {
-         nbdofs += fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
+         nbdofs += fec->DofForGeometry(mesh->GetElementBaseGeometry(MapElement(i)));
          bdofs[i+1] = nbdofs;
       }
    }
@@ -1708,25 +1830,48 @@ void FiniteElementSpace::GetElementDofs(int i, Array<int> &dofs) const
       dim = mesh->Dimension();
       nv = fec->DofForGeometry(Geometry::POINT);
       ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
-      nb = (dim > 0) ? fec->DofForGeometry(mesh->GetElementBaseGeometry(i)) : 0;
+      nb = (dim > 0) ? fec->DofForGeometry(mesh->GetElementBaseGeometry(MapElement(i))) : 0;
       if (nv > 0)
       {
-         mesh->GetElementVertices(i, V);
+         mesh->GetElementVertices(MapElement(i), V);
+         if(subdomain) // We have to map back to subdomain
+         {
+            for(int vi = 0; vi < V.Size(); vi++)
+            {
+               V[vi] = MapVertexBack(V[vi]);
+            }
+         }
       }
       if (ne > 0)
       {
-         mesh->GetElementEdges(i, E, Eo);
+         mesh->GetElementEdges(MapElement(i), E, Eo);
+         if(subdomain) // We have to map back to subdomain
+         {
+            for(int ei = 0; ei < E.Size(); ei++)
+            {
+               E[ei] = MapEdgeBack(E[ei]);
+            }
+         }
       }
       nfd = 0;
+
       if (dim == 3)
       {
-         if (fec->HasFaceDofs(mesh->GetElementBaseGeometry(i)))
+         if (fec->HasFaceDofs(mesh->GetElementBaseGeometry(MapElement(i))))
          {
-            mesh->GetElementFaces(i, F, Fo);
+            mesh->GetElementFaces(MapElement(i), F, Fo);
             for (k = 0; k < F.Size(); k++)
             {
                nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
             }
+
+            if(subdomain) // We have to map back to subdomain
+            {
+               for(int fi = 0; fi < F.Size(); fi++)
+               {
+                  F[fi] = MapFaceBack(F[fi]);
+               }
+            }
          }
       }
       nd = V.Size() * nv + E.Size() * ne + nfd + nb;
@@ -1767,9 +1912,9 @@ void FiniteElementSpace::GetElementDofs(int i, Array<int> &dofs) const
       {
          for (k = 0; k < F.Size(); k++)
          {
-            ind = fec->DofOrderForOrientation(mesh->GetFaceBaseGeometry(F[k]),
+            ind = fec->DofOrderForOrientation(mesh->GetFaceBaseGeometry(MapFace(F[k])),
                                               Fo[k]);
-            nf = fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
+            nf = fec->DofForGeometry(mesh->GetFaceBaseGeometry(MapFace(F[k])));
             for (j = 0; j < nf; j++)
             {
                if (ind[j] < 0)
@@ -1802,7 +1947,7 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
                "Invalid element id " << i << ", maximum allowed " << mesh->GetNE()-1);
 
    const FiniteElement *FE =
-      fec->FiniteElementForGeometry(mesh->GetElementBaseGeometry(i));
+      fec->FiniteElementForGeometry(mesh->GetElementBaseGeometry(MapElement(i)));
 
    if (NURBSext)
    {
@@ -1828,20 +1973,35 @@ void FiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
       nv = fec->DofForGeometry(Geometry::POINT);
       if (nv > 0)
       {
-         mesh->GetBdrElementVertices(i, V);
+         mesh->GetBdrElementVertices(MapBoundary(i), V);
+         if(subdomain) // We have to map back to subdomain
+         {
+            for(int vi = 0; vi < V.Size(); vi++)
+            {
+               V[vi] = MapVertexBack(V[vi]);
+            }
+         }
       }
       ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
       if (ne > 0)
       {
-         mesh->GetBdrElementEdges(i, E, Eo);
+         mesh->GetBdrElementEdges(MapBoundary(i), E, Eo);
+         if(subdomain) // We have to map back to subdomain
+         {
+            for(int ei = 0; ei < E.Size(); ei++)
+            {
+               E[ei] = MapEdgeBack(E[ei]);
+            }
+         }
       }
       nd = V.Size() * nv + E.Size() * ne;
       nf = (dim == 3) ? (fec->DofForGeometry(
-                            mesh->GetBdrElementBaseGeometry(i))) : (0);
+                            mesh->GetBdrElementBaseGeometry(MapBoundary(i)))) : (0);
       if (nf > 0)
       {
          nd += nf;
-         mesh->GetBdrElementFace(i, &iF, &oF);
+         mesh->GetBdrElementFace(MapBoundary(i), &iF, &oF);
+         iF = MapFaceBack(iF);
       }
       dofs.SetSize(nd);
       if (nv > 0)
@@ -1879,7 +2039,7 @@ void FiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
       {
          ne = nv + ne * E.Size();
          ind = fec->DofOrderForOrientation(
-                  mesh->GetBdrElementBaseGeometry(i), oF);
+                  mesh->GetBdrElementBaseGeometry(MapBoundary(i)), oF);
          for (j = 0; j < nf; j++)
          {
             if (ind[j] < 0)
@@ -1914,11 +2074,19 @@ void FiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs) const
       ne = (dim > 1) ? fec->DofForGeometry(Geometry::SEGMENT) : 0;
       if (nv > 0)
       {
-         mesh->GetFaceVertices(i, V);
+         mesh->GetFaceVertices(MapFace(i), V);
+         for(int vi = 0; vi < V.Size(); vi++)
+         {
+            V[vi] = MapVertexBack(V[vi]);
+         }
       }
       if (ne > 0)
       {
-         mesh->GetFaceEdges(i, E, Eo);
+         mesh->GetFaceEdges(MapFace(i), E, Eo);
+         for(int ei = 0; ei < E.Size(); ei++)
+         {
+            E[ei] = MapEdgeBack(E[ei]);
+         }
       }
       nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
       nd = V.Size() * nv + E.Size() * ne + nf;
@@ -1971,7 +2139,11 @@ void FiniteElementSpace::GetEdgeDofs(int i, Array<int> &dofs) const
    nv = fec->DofForGeometry(Geometry::POINT);
    if (nv > 0)
    {
-      mesh->GetEdgeVertices(i, V);
+      mesh->GetEdgeVertices(MapEdge(i), V);
+      for(int vi = 0; vi < V.Size(); vi++)
+      {
+         V[vi] = MapVertexBack(V[vi]);
+      }
    }
    ne = fec->DofForGeometry(Geometry::SEGMENT);
    dofs.SetSize(2*nv+ne);
@@ -1994,11 +2166,9 @@ void FiniteElementSpace::GetEdgeDofs(int i, Array<int> &dofs) const
 
 void FiniteElementSpace::GetVertexDofs(int i, Array<int> &dofs) const
 {
-   int j, nv;
-
-   nv = fec->DofForGeometry(Geometry::POINT);
+   const int nv = fec->DofForGeometry(Geometry::POINT);
    dofs.SetSize(nv);
-   for (j = 0; j < nv; j++)
+   for (int j = 0; j < nv; j++)
    {
       dofs[j] = i*nv+j;
    }
@@ -2006,12 +2176,11 @@ void FiniteElementSpace::GetVertexDofs(int i, Array<int> &dofs) const
 
 void FiniteElementSpace::GetElementInteriorDofs (int i, Array<int> &dofs) const
 {
-   int j, k, nb;
    if (mesh->Dimension() == 0) { dofs.SetSize(0); return; }
-   nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
+   const int nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (MapElement(i)));
    dofs.SetSize (nb);
-   k = nvdofs + nedofs + nfdofs + bdofs[i];
-   for (j = 0; j < nb; j++)
+   const int k = nvdofs + nedofs + nfdofs + bdofs[i];
+   for (int j = 0; j < nb; j++)
    {
       dofs[j] = k + j;
    }
@@ -2019,11 +2188,9 @@ void FiniteElementSpace::GetElementInteriorDofs (int i, Array<int> &dofs) const
 
 void FiniteElementSpace::GetEdgeInteriorDofs (int i, Array<int> &dofs) const
 {
-   int j, k, ne;
-
-   ne = fec -> DofForGeometry (Geometry::SEGMENT);
+   const int ne = fec -> DofForGeometry (Geometry::SEGMENT);
    dofs.SetSize (ne);
-   for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
+   for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
    {
       dofs[j] = k;
    }
@@ -2031,13 +2198,11 @@ void FiniteElementSpace::GetEdgeInteriorDofs (int i, Array<int> &dofs) const
 
 void FiniteElementSpace::GetFaceInteriorDofs (int i, Array<int> &dofs) const
 {
-   int j, k, nf;
-
-   nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
+   const int nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
    dofs.SetSize (nf);
    if (nf > 0)
    {
-      for (j = 0, k = nvdofs+nedofs+fdofs[i]; j < nf; j++, k++)
+      for (int j = 0, k = nvdofs+nedofs+fdofs[i]; j < nf; j++, k++)
       {
          dofs[j] = k;
       }
@@ -2059,7 +2224,7 @@ const FiniteElement *FiniteElementSpace::GetBE (int i) const
       case 3:
       default:
          BE = fec->FiniteElementForGeometry(
-                 mesh->GetBdrElementBaseGeometry(i));
+                 mesh->GetBdrElementBaseGeometry(MapElement(i)));
    }
 
    if (NURBSext)
@@ -2084,7 +2249,7 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
          break;
       case 3:
       default:
-         fe = fec->FiniteElementForGeometry(mesh->GetFaceBaseGeometry(i));
+         fe = fec->FiniteElementForGeometry(mesh->GetFaceBaseGeometry(MapElement(i)));
    }
 
    if (NURBSext)
@@ -2162,6 +2327,7 @@ void FiniteElementSpace::Destroy()
       delete [] fdofs;
    }
    RemoveCeedBasisAndRestriction(this);
+   if(own_subdomain) delete subdomain;
 }
 
 void FiniteElementSpace::GetTransferOperator(
@@ -2374,6 +2540,8 @@ void FiniteElementSpace::Save(std::ostream &out) const
       }
       out << "End: MFEM FiniteElementSpace v1.0\n";
    }
+
+   /// TODO (fe.subspace) store subspace information
 }
 
 FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
@@ -2472,6 +2640,8 @@ FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
       }
    }
 
+   /// TODO (fe.subspace) load subspace information
+
    Constructor(m, NURBSext, r_fec, vdim, ord);
 
    return r_fec;
@@ -2482,7 +2652,7 @@ void QuadratureSpace::Construct()
 {
    // protected method
    int offset = 0;
-   const int num_elem = mesh->GetNE();
+   const int num_elem = subdomain ? subdomain->element_map.Size() : mesh->GetNE();
    element_offsets = new int[num_elem + 1];
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
@@ -2490,8 +2660,9 @@ void QuadratureSpace::Construct()
    }
    for (int i = 0; i < num_elem; i++)
    {
+      const auto i_mesh = subdomain ? subdomain->element_map[i] : i;
       element_offsets[i] = offset;
-      int geom = mesh->GetElementBaseGeometry(i);
+      int geom = mesh->GetElementBaseGeometry(i_mesh);
       if (int_rule[geom] == NULL)
       {
          int_rule[geom] = &IntRules.Get(geom, order);
@@ -2529,6 +2700,8 @@ void QuadratureSpace::Save(std::ostream &out) const
    out << "QuadratureSpace\n"
        << "Type: default_quadrature\n"
        << "Order: " << order << '\n';
+
+   // TODO (fe.subspace) store subspace information
 }
 
 
@@ -2753,8 +2926,8 @@ L2ProjectionGridTransfer::L2Projection::L2Projection(
    MFEM_VERIFY(mesh_ho->GetNumGeometries(mesh_ho->Dimension()) <= 1,
                "mixed meshes are not supported");
 
-   // If the local mesh is empty, skip all computations
-   if (mesh_ho->GetNE() == 0) { return; }
+   // If the local space is empty, skip all computations
+   if (fes_ho_.GetNE() == 0) { return; }
 
    const FiniteElement *fe_lor = fes_lor.GetFE(0);
    const FiniteElement *fe_ho = fes_ho.GetFE(0);

fem/fespace.hpp

@@ -82,6 +82,101 @@ class QuadratureInterpolator;
 class FaceQuadratureInterpolator;
 
 
+// Holds all necessary information to map stuff to a subdomain of the mesh.
+struct SubdomainExtension
+{
+    // The mappings from the 1..N elements on subdomain back to the mesh
+    Array<int> vertex_map, edge_map, face_map, element_map;
+
+    Array<int> boundary_map;
+};
+
+// This construction conserves the original ordering of entities.
+inline SubdomainExtension* SubdomainFromAttributes(Mesh* mesh, Array<int> attributes)
+{
+   attributes.Sort();
+
+   auto subdomain = new SubdomainExtension();
+
+   // Add elements
+   for(int e = 0; e < mesh->GetNE(); e++)
+   {
+      if(attributes.FindSorted(mesh->GetAttribute(e)) != -1)
+      {
+         subdomain->element_map.Append(e);
+      }
+   }
+
+   // Add boundary
+   for(int be = 0; be < mesh->GetNBE(); be++)
+   {
+      int e, info;
+      mesh->GetBdrElementAdjacentElement(be, e, info);
+      if(attributes.FindSorted(mesh->GetAttribute(e)) != -1)
+      {
+         subdomain->boundary_map.Append(be);
+      }
+   }
+
+   Array<int> elements;
+   {
+      // Add faces
+      auto face2el = mesh->GetFaceToElementTable();
+      for(int f = 0; f < mesh->GetNFaces(); f++)
+      {
+         face2el->GetRow(f, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(mesh->GetAttribute(e)) != -1)
+            {
+               subdomain->face_map.Append(f);
+               break;
+            }
+         }
+      }
+      delete face2el;
+   }
+
+   {
+      auto edge2el = mesh->GetEdgeToElementTable();
+
+      // Add edges
+      for(int edge = 0; edge < mesh->GetNEdges(); edge++)
+      {
+         edge2el->GetRow(edge, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(mesh->GetAttribute(e)) != -1)
+            {
+               subdomain->edge_map.Append(edge);
+               break;
+            }
+         }
+      }
+      delete edge2el;
+   }
+
+   {
+      auto vertex2el = mesh->GetVertexToElementTable();
+      // Add vertices
+      for(int v = 0; v < mesh->GetNV(); v++)
+      {
+         vertex2el->GetRow(v, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(mesh->GetAttribute(e)) != -1)
+            {
+               subdomain->vertex_map.Append(v);
+               break;
+            }
+         }
+      }
+
+      delete vertex2el;
+   }
+   return subdomain;
+}
+
 /** @brief Class FiniteElementSpace - responsible for providing FEM view of the
     mesh, mainly managing the set of degrees of freedom. */
 class FiniteElementSpace
@@ -120,6 +215,10 @@ protected:
    NURBSExtension *NURBSext;
    int own_ext;
 
+   // In case that the FiniteElementSpace is just definied on a part of the mesh
+   SubdomainExtension *subdomain = nullptr;
+   bool own_subdomain = false;
+
    /** Matrix representing the prolongation from the global conforming dofs to
        a set of intermediate partially conforming dofs, e.g. the dofs associated
        with a "cut" space on a non-conforming mesh. */
@@ -298,7 +397,9 @@ public:
 
    FiniteElementSpace(Mesh *mesh,
                       const FiniteElementCollection *fec,
-                      int vdim = 1, int ordering = Ordering::byNODES)
+                      int vdim = 1, int ordering = Ordering::byNODES,
+                      SubdomainExtension* subdomain_ = NULL)
+   : subdomain(subdomain_)
    { Constructor(mesh, NULL, fec, vdim, ordering); }
 
    /// Construct a NURBS FE space based on the given NURBSExtension, @a ext.
@@ -313,6 +414,142 @@ public:
    /// Returns the mesh
    inline Mesh *GetMesh() const { return mesh; }
 
+   /// Some helpers to get the correct indices on the full mesh.
+   inline int MapElement(int i) const
+   { return subdomain ? subdomain->element_map[i] : i; }
+   inline int MapElementBack(int i) const
+   {
+      if(!subdomain)
+      {
+         return i;
+      }
+      for(int bi = 0; bi < subdomain->element_map.Size(); bi++)
+      {
+         if(subdomain->element_map[bi] == i)
+         {
+            return bi;
+         }
+      }
+      return -1;
+   }
+   inline int MapFace(int i) const
+   { return subdomain ? subdomain->face_map[i] : i; }
+   inline int MapFaceBack(int i) const
+   {
+      if(!subdomain)
+      {
+         return i;
+      }
+      for(int bi = 0; bi < subdomain->face_map.Size(); bi++)
+      {
+         if(subdomain->face_map[bi] == i)
+         {
+            return bi;
+         }
+      }
+      return -1;
+   }
+   inline int MapEdge(int i) const
+   { return subdomain ? subdomain->edge_map[i] : i; }
+   inline int MapEdgeBack(int i) const
+   {
+      if(!subdomain)
+      {
+         return i;
+      }
+      for(int bi = 0; bi < subdomain->edge_map.Size(); bi++)
+      {
+         if(subdomain->edge_map[bi] == i)
+         {
+            return bi;
+         }
+      }
+      return -1;
+   }
+   inline int MapVertex(int i) const
+   { return subdomain ? subdomain->vertex_map[i] : i; }
+   inline int MapVertexBack(int i) const
+   {
+      if(!subdomain)
+      {
+         return i;
+      }
+      for(int vi = 0; vi < subdomain->vertex_map.Size(); vi++)
+      {
+         if(subdomain->vertex_map[vi] == i)
+         {
+            return vi;
+         }
+      }
+      return -1;
+   }
+   inline int MapBoundary(int i) const
+   { return subdomain ? subdomain->boundary_map[i] : i; }
+   inline int MapBoundaryBack(int i) const
+   {
+      if(!subdomain)
+      {
+         return i;
+      }
+      for(int bi = 0; bi < subdomain->boundary_map.Size(); bi++)
+      {
+         if(subdomain->boundary_map[bi] == i)
+         {
+            return bi;
+         }
+      }
+      return -1;
+   }
+   inline int MapEntityBack(int entity, int i) const
+   {
+      switch(entity)
+      {
+         case 0:
+         return MapVertexBack(i);
+
+         case 1:
+         return MapEdgeBack(i);
+
+         case 2:
+         return MapFaceBack(i);
+
+         case 3:
+         return MapElementBack(i);
+
+         default:
+         MFEM_ABORT("Unknown entity " << i);
+      }
+      return -1;
+   }
+
+   inline int MapEntity(int entity, int i) const
+   {
+      switch(entity)
+      {
+         case 0:
+         return MapVertex(i);
+
+         case 1:
+         return MapEdge(i);
+
+         case 2:
+         return MapFace(i);
+
+         case 3:
+         return MapElement(i);
+
+         default:
+         MFEM_ABORT("Unknown entity " << i);
+      }
+      return -1;
+   }
+
+   inline bool UsesSubdomain() const
+   { return subdomain != nullptr; }
+
+   inline SubdomainExtension* Subdomain() const
+   { return subdomain; }
+
    const NURBSExtension *GetNURBSext() const { return NURBSext; }
    NURBSExtension *GetNURBSext() { return NURBSext; }
    NURBSExtension *StealNURBSext();
@@ -432,20 +669,32 @@ public:
    /// Number of all scalar face-interior dofs
    int GetNFDofs() const { return nfdofs; }
 
-   /// Returns number of vertices in the mesh.
-   inline int GetNV() const { return mesh->GetNV(); }
+   /// Returns number of vertices in the space.
+   inline int GetNV() const
+   { return subdomain ? subdomain->vertex_map.Size() :  mesh->GetNV(); }
 
-   /// Returns number of elements in the mesh.
-   inline int GetNE() const { return mesh->GetNE(); }
+   /// Returns number of elements in the space.
+   inline int GetNE() const
+   { return subdomain ? subdomain->element_map.Size() : mesh->GetNE(); }
+
+   /// Returns number of edges in the space.
+   inline int GetNEdges() const
+   { return subdomain ? subdomain->edge_map.Size() : mesh->GetNEdges(); }
+
+   /// Returns number of faces in the space.
+   inline int GetNFaces() const
+   { return subdomain ? subdomain->face_map.Size() : mesh->GetNFaces(); }
 
    /// Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
    /** The co-dimension 1 entities are those that have dimension 1 less than the
        mesh dimension, e.g. for a 2D mesh, the faces are the 1D entities, i.e.
        the edges. */
-   inline int GetNF() const { return mesh->GetNumFaces(); }
+   inline int GetNF() const
+   { return subdomain ? subdomain->face_map.Size() : mesh->GetNumFaces(); }
 
    /// Returns number of boundary elements in the mesh.
-   inline int GetNBE() const { return mesh->GetNBE(); }
+   inline int GetNBE() const
+   { return subdomain ? subdomain->boundary_map.Size() : mesh->GetNBE(); }
 
    /// Returns the number of faces according to the requested type.
    /** If type==Boundary returns only the "true" number of boundary faces
@@ -458,32 +707,44 @@ public:
 
    /// Returns the type of element i.
    inline int GetElementType(int i) const
-   { return mesh->GetElementType(i); }
+   { return mesh->GetElementType(MapElement(i)); }
 
    /// Returns the vertices of element i.
    inline void GetElementVertices(int i, Array<int> &vertices) const
-   { mesh->GetElementVertices(i, vertices); }
+   { mesh->GetElementVertices(MapElement(i), vertices); }
 
    /// Returns the type of boundary element i.
    inline int GetBdrElementType(int i) const
-   { return mesh->GetBdrElementType(i); }
+   { return mesh->GetBdrElementType(MapBoundary(i)); }
 
    /// Returns ElementTransformation for the @a i-th element.
    ElementTransformation *GetElementTransformation(int i) const
-   { return mesh->GetElementTransformation(i); }
+   { return mesh->GetElementTransformation(MapElement(i)); }
 
    /** @brief Returns the transformation defining the @a i-th element in the
        user-defined variable @a ElTr. */
    void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
-   { mesh->GetElementTransformation(i, ElTr); }
+   { mesh->GetElementTransformation(MapElement(i), ElTr); }
+
+   ElementTransformation *GetFaceTransformation(int i) const
+   { return mesh->GetFaceTransformation(i); }
 
    /// Returns ElementTransformation for the @a i-th boundary element.
    ElementTransformation *GetBdrElementTransformation(int i) const
-   { return mesh->GetBdrElementTransformation(i); }
+   { return mesh->GetBdrElementTransformation(MapBoundary(i)); }
 
-   int GetAttribute(int i) const { return mesh->GetAttribute(i); }
+   FaceElementTransformations *GetBdrFaceTransformations(int i) const
+   { return mesh->GetBdrFaceTransformations(MapBoundary(i)); }
 
-   int GetBdrAttribute(int i) const { return mesh->GetBdrAttribute(i); }
+   FaceElementTransformations *GetInteriorFaceTransformations(int i) const
+   { return mesh->GetInteriorFaceTransformations(MapFace(i)); }
+
+   FaceElementTransformations *GetFaceElementTransformations(int i) const
+   { return mesh->GetFaceElementTransformations(MapFace(i)); }
+
+   int GetAttribute(int i) const { return mesh->GetAttribute(MapElement(i)); }
+
+   int GetBdrAttribute(int i) const { return mesh->GetBdrAttribute(MapBoundary(i)); }
 
    /// Returns indexes of degrees of freedom in array dofs for i'th element.
    virtual void GetElementDofs(int i, Array<int> &dofs) const;
@@ -506,7 +767,7 @@ public:
    void GetFaceInteriorDofs(int i, Array<int> &dofs) const;
 
    int GetNumElementInteriorDofs(int i) const
-   { return fec->DofForGeometry(mesh->GetElementBaseGeometry(i)); }
+   { return fec->DofForGeometry(mesh->GetElementBaseGeometry(MapElement(i))); }
 
    void GetEdgeInteriorDofs(int i, Array<int> &dofs) const;
 
@@ -755,6 +1016,10 @@ protected:
    // element_offsets, and size.
    void Construct();
 
+   //
+   SubdomainExtension *subdomain = nullptr;
+   bool own_subdomain = false;
+
 public:
    /// Create a QuadratureSpace based on the global rules from #IntRules.
    QuadratureSpace(Mesh *mesh_, int order_)
@@ -763,7 +1028,8 @@ public:
    /// Read a QuadratureSpace from the stream @a in.
    QuadratureSpace(Mesh *mesh_, std::istream &in);
 
-   virtual ~QuadratureSpace() { delete [] element_offsets; }
+   virtual ~QuadratureSpace()
+   { delete [] element_offsets; if(own_subdomain) delete subdomain; }
 
    /// Return the total number of quadrature points.
    int GetSize() const { return size; }
@@ -775,11 +1041,15 @@ public:
    inline Mesh *GetMesh() const { return mesh; }
 
    /// Returns number of elements in the mesh.
-   inline int GetNE() const { return mesh->GetNE(); }
+   inline int GetNE() const
+   { return subdomain ? subdomain->element_map.Size() : mesh->GetNE(); }
 
    /// Get the IntegrationRule associated with mesh element @a idx.
    const IntegrationRule &GetElementIntRule(int idx) const
-   { return *int_rule[mesh->GetElementBaseGeometry(idx)]; }
+   {
+      auto real_idx = subdomain ? subdomain->element_map[idx] : idx;
+      return *int_rule[mesh->GetElementBaseGeometry(real_idx)];
+   }
 
    /// Write the QuadratureSpace to the stream @a out.
    void Save(std::ostream &out) const;

fem/linearform.cpp

@@ -132,7 +132,7 @@ void LinearForm::Assemble()
 
       for (i = 0; i < fes -> GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
          fes -> GetBdrElementVDofs (i, vdofs);
          eltrans = fes -> GetBdrElementTransformation (i);
@@ -173,12 +173,12 @@ void LinearForm::Assemble()
          }
       }
 
-      for (i = 0; i < mesh->GetNBE(); i++)
+      for (i = 0; i < fes->GetNBE(); i++)
       {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
+         const int bdr_attr = fes->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
-         tr = mesh->GetBdrFaceTransformations(i);
+         tr = fes->GetBdrFaceTransformations(i);
          if (tr != NULL)
          {
             fes -> GetElementVDofs (tr -> Elem1No, vdofs);

fem/pfespace.cpp

@@ -60,8 +60,10 @@ ParFiniteElementSpace::ParFiniteElementSpace(
 }
 
 ParFiniteElementSpace::ParFiniteElementSpace(
-   ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering)
-   : FiniteElementSpace(pm, f, dim, ordering)
+   ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering
+   , ParSubdomainExtension *psubdomain_)
+   : FiniteElementSpace(pm, f, dim, ordering, psubdomain_)
+   , psubdomain(psubdomain_)
 {
    ParInit(pm);
 }
@@ -154,19 +156,50 @@ void ParFiniteElementSpace::Construct()
       ngedofs = ngfdofs = 0;
 
       // calculate number of ghost DOFs
-      ngvdofs = pncmesh->GetNGhostVertices()
-                * fec->DofForGeometry(Geometry::POINT);
-
-      if (pmesh->Dimension() > 1)
+      if(!UsesSubdomain())
       {
-         ngedofs = pncmesh->GetNGhostEdges()
-                   * fec->DofForGeometry(Geometry::SEGMENT);
+         ngvdofs = pncmesh->GetNGhostVertices()
+                  * fec->DofForGeometry(Geometry::POINT);
+
+         if (pmesh->Dimension() > 1)
+         {
+            ngedofs = pncmesh->GetNGhostEdges()
+                     * fec->DofForGeometry(Geometry::SEGMENT);
+         }
+
+         if (pmesh->Dimension() > 2)
+         {
+            int stride = fec->DofForGeometry(Geometry::SQUARE);
+            ngfdofs = pncmesh->GetNGhostFaces() * stride;
+         }
       }
-
-      if (pmesh->Dimension() > 2)
+      else // FIXME this does not work.
       {
-         int stride = fec->DofForGeometry(Geometry::SQUARE);
-         ngfdofs = pncmesh->GetNGhostFaces() * stride;
+         ngvdofs = std::count_if(psubdomain->vertex_map.begin(),
+                                 psubdomain->vertex_map.end(),
+                                 [this](int vertex_in_mesh){
+                                    return vertex_in_mesh >= pmesh->GetNV();
+                                 })*fec->DofForGeometry(Geometry::POINT);
+
+         if (pmesh->Dimension() > 1)
+         {
+            ngedofs = std::count_if(psubdomain->edge_map.begin(),
+                                 psubdomain->edge_map.end(),
+                                 [this](int edge_in_mesh){
+                                    return edge_in_mesh >= pmesh->GetNEdges();
+                                 })*fec->DofForGeometry(Geometry::SEGMENT);
+         }
+
+         if (pmesh->Dimension() > 2)
+         {
+            ngfdofs = std::count_if(psubdomain->face_map.begin(),
+                                 psubdomain->face_map.end(),
+                                 [this](int face_in_mesh){
+                                    return face_in_mesh >= pmesh->GetNFaces();
+                                 })*fec->DofForGeometry(Geometry::SQUARE);
+         }
+
+         std::cout << "Ghosts: " << ngvdofs << " " << ngedofs << " " << ngfdofs << std::endl;
       }
 
       // total number of ghost DOFs. Ghost DOFs start at index 'ndofs', i.e.,
@@ -257,12 +290,58 @@ void ParFiniteElementSpace::GetGroupComm(
 
    // count the number of ldofs in all groups (excluding the local group 0)
    group_ldof_counter = 0;
-   for (gr = 1; gr < ng; gr++)
+   // Technically we do not need this differentiation, but the code is slightly
+   // faster this way.
+   if(!UsesSubdomain())
    {
-      group_ldof_counter += nvd * pmesh->GroupNVertices(gr);
-      group_ldof_counter += ned * pmesh->GroupNEdges(gr);
-      group_ldof_counter += ntd * pmesh->GroupNTriangles(gr);
-      group_ldof_counter += nqd * pmesh->GroupNQuadrilaterals(gr);
+      for (gr = 1; gr < ng; gr++)
+      {
+         group_ldof_counter += nvd * pmesh->GroupNVertices(gr);
+         group_ldof_counter += ned * pmesh->GroupNEdges(gr);
+         group_ldof_counter += ntd * pmesh->GroupNTriangles(gr);
+         group_ldof_counter += nqd * pmesh->GroupNQuadrilaterals(gr);
+      }
+   }
+   else
+   // In the subdomain case we just have to account for entities in the subdomain.
+   {
+      for (gr = 1; gr < ng; gr++)
+      {
+         for(int svi = 0; svi < pmesh->GroupNVertices(gr); svi++)
+         {
+            if(MapVertexBack(pmesh->GroupVertex(gr, svi)) >= 0)
+            {
+               group_ldof_counter += nvd;
+            }
+         }
+         for(int sei = 0; sei < pmesh->GroupNEdges(gr); sei++)
+         {
+            int ei, eo;
+            pmesh->GroupEdge(gr, sei, ei, eo);
+            if(MapEdgeBack(ei) >= 0)
+            {
+               group_ldof_counter += ned;
+            }
+         }
+         for(int stfi = 0; stfi < pmesh->GroupNTriangles(gr); stfi++)
+         {
+            int fi, fo;
+            pmesh->GroupTriangle(gr, stfi, fi, fo);
+            if(MapFaceBack(fi) >= 0)
+            {
+               group_ldof_counter += ntd;
+            }
+         }
+         for(int sqfi = 0; sqfi < pmesh->GroupNQuadrilaterals(gr); sqfi++)
+         {
+            int fi, fo;
+            pmesh->GroupQuadrilateral(gr, sqfi, fi, fo);
+            if(MapFaceBack(fi) >= 0)
+            {
+               group_ldof_counter += nqd;
+            }
+         }
+      }
    }
    if (ldof_type)
    {
@@ -289,7 +368,8 @@ void ParFiniteElementSpace::GetGroupComm(
       {
          for (j = 0; j < nv; j++)
          {
-            k = pmesh->GroupVertex(gr, j);
+            k = MapVertexBack(pmesh->GroupVertex(gr, j));
+            if(k < 0) continue;
 
             dofs.SetSize(nvd);
             m = nvd * k;
@@ -316,6 +396,8 @@ void ParFiniteElementSpace::GetGroupComm(
          for (j = 0; j < ne; j++)
          {
             pmesh->GroupEdge(gr, j, k, o);
+            k = MapEdgeBack(k);
+            if(k < 0) continue;
 
             dofs.SetSize(ned);
             m = nvdofs+k*ned;
@@ -354,6 +436,8 @@ void ParFiniteElementSpace::GetGroupComm(
          for (j = 0; j < nt; j++)
          {
             pmesh->GroupTriangle(gr, j, k, o);
+            k = MapFaceBack(k);
+            if(k < 0) continue;
 
             dofs.SetSize(ntd);
             m = nvdofs+nedofs+fdofs[k];
@@ -392,6 +476,8 @@ void ParFiniteElementSpace::GetGroupComm(
          for (j = 0; j < nq; j++)
          {
             pmesh->GroupQuadrilateral(gr, j, k, o);
+            k = MapFaceBack(k);
+            if(k < 0) continue;
 
             dofs.SetSize(nqd);
             m = nvdofs+nedofs+fdofs[k];
@@ -503,7 +589,7 @@ void ParFiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs) const
 
 const FiniteElement *ParFiniteElementSpace::GetFE(int i) const
 {
-   int ne = mesh->GetNE();
+   int ne = this->GetNE();
    if (i >= ne) { return GetFaceNbrFE(i - ne); }
    else { return FiniteElementSpace::GetFE(i); }
 }
@@ -1220,7 +1306,7 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    // Works for NC mesh where 'i' is an index returned by
    // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
    // the index of a ghost face.
-   MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
+   MFEM_ASSERT(Nonconforming() && i >= this->GetNF(), "");
    int el1, el2, inf1, inf2;
    pmesh->GetFaceElements(i, &el1, &el2);
    el2 = -1 - el2;
@@ -1247,7 +1333,7 @@ const FiniteElement *ParFiniteElementSpace::GetFaceNbrFE(int i) const
 {
    const FiniteElement *FE =
       fec->FiniteElementForGeometry(
-         pmesh->face_nbr_elements[i]->GetGeometryType());
+         pmesh->face_nbr_elements[MapFaceBack(i)]->GetGeometryType());
 
    if (NURBSext)
    {
@@ -1265,7 +1351,7 @@ const FiniteElement *ParFiniteElementSpace::GetFaceNbrFaceFE(int i) const
    // Works in tandem with GetFaceNbrFaceVDofs() defined above.
 
    MFEM_ASSERT(Nonconforming() && !NURBSext, "");
-   Geometry::Type face_geom = pmesh->GetFaceGeometryType(i);
+   Geometry::Type face_geom = pmesh->GetFaceGeometryType(MapFaceBack(i));
    return fec->FiniteElementForGeometry(face_geom);
 }
 
@@ -1281,7 +1367,7 @@ void ParFiniteElementSpace::Lose_Dof_TrueDof_Matrix()
 
 void ParFiniteElementSpace::ConstructTrueDofs()
 {
-   int i, gr, n = GetVSize();
+   const int n = GetVSize();
    GroupTopology &gt = pmesh->gtopo;
    gcomm = new GroupCommunicator(gt);
    Table &group_ldof = gcomm->GroupLDofTable();
@@ -1296,18 +1382,18 @@ void ParFiniteElementSpace::ConstructTrueDofs()
    ldof_group = 0;
    ldof_ltdof = -1;
 
-   for (gr = 1; gr < group_ldof.Size(); gr++)
+   for (int gr = 1; gr < group_ldof.Size(); gr++)
    {
       const int *ldofs = group_ldof.GetRow(gr);
       const int nldofs = group_ldof.RowSize(gr);
-      for (i = 0; i < nldofs; i++)
+      for (int i = 0; i < nldofs; i++)
       {
          ldof_group[ldofs[i]] = gr;
       }
 
       if (!gt.IAmMaster(gr)) // we are not the master
       {
-         for (i = 0; i < nldofs; i++)
+         for (int i = 0; i < nldofs; i++)
          {
             ldof_ltdof[ldofs[i]] = -2;
          }
@@ -1316,7 +1402,7 @@ void ParFiniteElementSpace::ConstructTrueDofs()
 
    // count ltdof_size
    ltdof_size = 0;
-   for (i = 0; i < n; i++)
+   for (int i = 0; i < n; i++)
    {
       if (ldof_ltdof[i] == -1)
       {
@@ -2028,6 +2114,8 @@ int ParFiniteElementSpace
 {
    bool dg = (nvdofs == 0 && nedofs == 0 && nfdofs == 0);
 
+   std::cout << "BuildParallelConformingInterpolation" << std::endl;
+
 #ifdef MFEM_PMATRIX_STATS
    n_msgs_sent = n_msgs_recv = 0;
    n_rows_sent = n_rows_recv = n_rows_fwd = 0;
@@ -2043,6 +2131,7 @@ int ParFiniteElementSpace
       Array<int> master_dofs, slave_dofs;
 
       // loop through *all* master edges/faces, constrain their slaves
+      // TODO: Why do we start here at 0 and in serial at 1?
       for (int entity = 0; entity <= 2; entity++)
       {
          const NCMesh::NCList &list = pncmesh->GetNCList(entity);
@@ -2054,19 +2143,34 @@ int ParFiniteElementSpace
          // process masters that we own or that affect our edges/faces
          for (int mi = 0; mi < list.masters.Size(); mi++)
          {
-            const NCMesh::Master &mf = list.masters[mi];
+            const NCMesh::Master &master = list.masters[mi];
+
+            // Skip DoF outside of subdomain
+            if(UsesSubdomain())
+            {
+               // FIXME Why is the element sometimes -1?
+               if(master.element<0) std::cout << "Step 1 master: " << entity  << " " << master.element << std::endl;
+               if(MapElementBack(pncmesh->GetLeafElement(master.element)) == -1)
+               {
+                  continue;
+               }
+            }
+
+            // FIXME: This does not work for NC tets.
+            auto master_index = MapEntityBack(entity, master.index);
+            if(master_index < 0) continue;
 
             // get master DOFs
-            pncmesh->IsGhost(entity, mf.index)
-            ? GetGhostDofs(entity, mf, master_dofs)
-            : GetEntityDofs(entity, mf.index, master_dofs);
+            pncmesh->IsGhost(entity, master.index)
+            ? GetGhostDofs(entity, master, master_dofs)
+            : GetEntityDofs(entity, master_index, master_dofs);
 
             if (!master_dofs.Size()) { continue; }
 
-            const FiniteElement* fe = fec->FiniteElementForGeometry(mf.Geom());
+            const FiniteElement* fe = fec->FiniteElementForGeometry(master.Geom());
             if (!fe) { continue; }
 
-            switch (mf.Geom())
+            switch (master.Geom())
             {
                case Geometry::SQUARE:   T.SetFE(&QuadrilateralFE); break;
                case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
@@ -2075,15 +2179,31 @@ int ParFiniteElementSpace
             }
 
             // constrain slaves that exist in our mesh
-            for (int si = mf.slaves_begin; si < mf.slaves_end; si++)
+            for (int si = master.slaves_begin; si < master.slaves_end; si++)
             {
-               const NCMesh::Slave &sf = list.slaves[si];
-               if (pncmesh->IsGhost(entity, sf.index)) { continue; }
+               const NCMesh::Slave &slave = list.slaves[si];
 
-               GetEntityDofs(entity, sf.index, slave_dofs, mf.Geom());
+               // Skip DoF outside of subdomain
+               if(UsesSubdomain())
+               {
+                  // FIXME Why is the element sometimes -1?
+                  if(slave.element<0) std::cout << "Step 1 slave: " << entity  << " " << slave.element << std::endl;
+                  if(MapElementBack(pncmesh->GetLeafElement(slave.element)) == -1)
+                  {
+                     continue;
+                  }
+               }
+
+               // FIXME: This does not work for NC tets.
+               auto slave_index = MapEntityBack(entity, slave.index);
+               if(slave_index < 0) continue;
+
+               if (pncmesh->IsGhost(entity, slave.index)) { continue; }
+
+               GetEntityDofs(entity, slave_index, slave_dofs, master.Geom());
                if (!slave_dofs.Size()) { continue; }
 
-               list.OrientedPointMatrix(sf, T.GetPointMat());
+               list.OrientedPointMatrix(slave, T.GetPointMat());
                fe->GetLocalInterpolation(T, I);
 
                // make each slave DOF dependent on all master DOFs
@@ -2094,7 +2214,7 @@ int ParFiniteElementSpace
 
       deps.Finalize();
    }
-
+std::cout << "1 done" << std::endl;
    // *** STEP 2: initialize group and owner ID for each DOF ***
 
    Array<GroupId> dof_group(total_dofs);
@@ -2125,10 +2245,24 @@ int ParFiniteElementSpace
 
                if (id.index < 0) { continue; }
 
+               // Skip DoF outside of subdomain
+               if(UsesSubdomain())
+               {
+                  // FIXME Why is the element sometimes -1?
+                  if(id.element<0) std::cout << "Step 2: " << entity << " " << id.element << std::endl;
+                  if(MapElementBack(pncmesh->GetLeafElement(id.element)) == -1)
+                  {
+                     continue;
+                  }
+               }
+
+               auto index = MapEntityBack(entity, id.index);
+               if(index < 0) continue;
+
                GroupId owner = pncmesh->GetEntityOwnerId(entity, id.index);
                GroupId group = pncmesh->GetEntityGroupId(entity, id.index);
 
-               GetBareDofs(entity, id.index, dofs);
+               GetBareDofs(entity, index, dofs);
 
                for (int j = 0; j < dofs.Size(); j++)
                {
@@ -2140,7 +2274,7 @@ int ParFiniteElementSpace
          }
       }
    }
-
+std::cout << "2 done" << std::endl;
    // *** STEP 3: count true DOFs and calculate P row/column partitions ***
 
    Array<bool> finalized(total_dofs);
@@ -2175,7 +2309,7 @@ int ParFiniteElementSpace
       dof_tdof->SetSize(ndofs*vdim);
       *dof_tdof = -1;
    }
-
+std::cout << "fine?" << std::endl;
    std::vector<PMatrixRow> pmatrix(total_dofs);
 
    bool bynodes = (ordering == Ordering::byNODES);
@@ -2212,7 +2346,7 @@ int ParFiniteElementSpace
          tdof++;
       }
    }
-
+std::cout << "3 done" << std::endl;
    // send identity rows
    NeighborRowMessage::IsendAll(send_msg.back(), MyComm);
 #ifdef MFEM_PMATRIX_STATS
@@ -2234,6 +2368,7 @@ int ParFiniteElementSpace
 
    while (num_finalized < ndofs)
    {
+      //std::cout << "Finalized "  << num_finalized << "/"  << ndofs << std::endl;
       // prepare a new round of send buffers
       if (send_msg.back().size())
       {
@@ -2254,7 +2389,10 @@ int ParFiniteElementSpace
          for (unsigned i = 0; i < rows.size(); i++)
          {
             const NeighborRowMessage::RowInfo &ri = rows[i];
-            int dof = PackDof(ri.entity, ri.index, ri.edof);
+            const auto index = MapEntityBack(ri.entity, ri.index);
+            if(index == -1) continue;
+
+            int dof = PackDof(ri.entity, index, ri.edof);
             pmatrix[dof] = ri.row;
 
             if (dof < ndofs && !finalized[dof]) { num_finalized++; }

fem/pfespace.hpp

@@ -24,6 +24,129 @@
 namespace mfem
 {
 
+// Holds all necessary information to map stuff to a subdomain of the mesh.
+struct ParSubdomainExtension : public SubdomainExtension
+{
+};
+
+// This construction conserves the original ordering of entities.
+inline ParSubdomainExtension* SubdomainFromAttributes(ParMesh* pmesh, Array<int> attributes)
+{
+   attributes.Sort();
+
+   auto subdomain = new ParSubdomainExtension();
+
+   // Add elements
+   for(int e = 0; e < pmesh->GetNE(); e++)
+   {
+      if(attributes.FindSorted(pmesh->GetAttribute(e)) != -1)
+      {
+         subdomain->element_map.Append(e);
+      }
+   }
+
+   // Add boundary
+   for(int be = 0; be < pmesh->GetNBE(); be++)
+   {
+      int e, info;
+      pmesh->GetBdrElementAdjacentElement(be, e, info);
+      if(attributes.FindSorted(pmesh->GetAttribute(e)) != -1)
+      {
+         subdomain->boundary_map.Append(be);
+      }
+   }
+
+   Array<int> elements;
+   {
+      // Add faces
+      auto face2el = pmesh->GetFaceToElementTable();
+      for(int f = 0; f < pmesh->GetNFaces(); f++)
+      {
+         face2el->GetRow(f, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(pmesh->GetAttribute(e)) != -1)
+            {
+               subdomain->face_map.Append(f);
+               break;
+            }
+         }
+      }
+      delete face2el;
+   }
+
+   {
+      auto edge2el = pmesh->GetEdgeToElementTable();
+
+      // Add edges
+      for(int edge = 0; edge < pmesh->GetNEdges(); edge++)
+      {
+         edge2el->GetRow(edge, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(pmesh->GetAttribute(e)) != -1)
+            {
+               subdomain->edge_map.Append(edge);
+               break;
+            }
+         }
+      }
+      delete edge2el;
+   }
+
+   {
+      auto vertex2el = pmesh->GetVertexToElementTable();
+      // Add vertices
+      for(int v = 0; v < pmesh->GetNV(); v++)
+      {
+         vertex2el->GetRow(v, elements);
+         for(int e : elements)
+         {
+            if(attributes.FindSorted(pmesh->GetAttribute(e)) != -1)
+            {
+               subdomain->vertex_map.Append(v);
+               break;
+            }
+         }
+      }
+
+      delete vertex2el;
+   }
+
+   // Some shared entities may be hidden, because their element is on a different process, so we have search for them.
+   // pmesh->ExchangeFaceNbrData();
+   // Array<int> entity_buf;
+   // Array<int> o_buf;
+   // for (int sf = 0; sf < pmesh->GetNSharedFaces(); sf++)
+   // {
+   //    const auto FT = pmesh->GetSharedFaceTransformations(sf, true);
+   //    if(attributes.FindSorted(FT->Elem2->Attribute) != -1)
+   //    {
+   //       //const auto f = FT->Face->ElementNo;
+   //       const auto f = pmesh->GetSharedFace(sf);
+   //       subdomain->face_map.Append(f);
+   //       if(f >= pmesh->GetNumFaces()) continue;
+   //       // TODO how to get these if the face is a ghost?
+   //       pmesh->GetFaceEdges(f, entity_buf, o_buf);
+   //       subdomain->edge_map.Append(entity_buf);
+   //       pmesh->GetFaceVertices(f, entity_buf);
+   //       subdomain->vertex_map.Append(entity_buf);
+   //    }
+   // }
+
+   // The last step messed up everything, so we have to do some cleaunup to guarantee uniqueness of the entities and that they are properly sorted
+   subdomain->vertex_map.Sort();
+   subdomain->vertex_map.Unique();
+
+   subdomain->edge_map.Sort();
+   subdomain->edge_map.Unique();
+
+   subdomain->face_map.Sort();
+   subdomain->face_map.Unique();
+
+   return subdomain;
+}
+
 /// Abstract parallel finite element space.
 class ParFiniteElementSpace : public FiniteElementSpace
 {
@@ -38,6 +161,9 @@ private:
        Not owned. */
    ParNCMesh *pncmesh;
 
+   // In case that the FiniteElementSpace is just definied on a part of the mesh
+   ParSubdomainExtension *psubdomain = nullptr;
+
    /// GroupCommunicator on the local VDofs. Owned.
    GroupCommunicator *gcomm;
 
@@ -230,7 +356,8 @@ public:
                          const FiniteElementCollection *f = NULL);
 
    ParFiniteElementSpace(ParMesh *pm, const FiniteElementCollection *f,
-                         int dim = 1, int ordering = Ordering::byNODES);
+                         int dim = 1, int ordering = Ordering::byNODES,
+                         ParSubdomainExtension *psubdomain = nullptr);
 
    /// Construct a NURBS FE space based on the given NURBSExtension, @a ext.
    /** The parameter @a ext will be deleted by this constructor, replaced by a

mesh/mesh.cpp

@@ -4923,6 +4923,48 @@ Table *Mesh::GetVertexToElementTable()
    return vert_elem;
 }
 
+
+Table *Mesh::GetEdgeToElementTable()
+{
+   MFEM_ASSERT(SpaceDimension() >= 2, "Edges just exist in dimension 2 and higher.");
+
+   Table *edge_elem = new Table;
+
+   if (el_to_edge == NULL)
+   {
+      el_to_edge = new Table;
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+   }
+
+   edge_elem->MakeI(NumOfEdges);
+
+   Array<int> edges;
+   for (int i = 0; i < NumOfElements; i++)
+   {
+      el_to_edge->GetRow(i, edges);
+      for (auto edge : edges)
+      {
+         edge_elem->AddAColumnInRow(edge);
+      }
+   }
+
+   edge_elem->MakeJ();
+
+   for (int i = 0; i < NumOfElements; i++)
+   {
+      el_to_edge->GetRow(i, edges);
+      for (auto edge : edges)
+      {
+         edge_elem->AddConnection(edge, i);
+      }
+   }
+
+   edge_elem->ShiftUpI();
+
+   return edge_elem;
+}
+
+
 Table *Mesh::GetFaceToElementTable() const
 {
    Table *face_elem = new Table;

mesh/mesh.hpp

@@ -1066,6 +1066,9 @@ public:
    ///  The returned Table must be destroyed by the caller
    Table *GetVertexToElementTable();
 
+   ///  The returned Table must be destroyed by the caller
+   Table *GetEdgeToElementTable();
+
    /** Return the "face"-element Table. Here "face" refers to face (3D),
        edge (2D), or vertex (1D).
        The returned Table must be destroyed by the caller. */

mesh/ncmesh.hpp

@@ -341,6 +341,10 @@ public:
    void GetElementFacesAttributes(int i, Array<int> &faces,
                                   Array<int> &fattr) const;
 
+   ///
+   int GetLeafElement(int i) const
+   { return elements[i].index; }
+
 
    /// I/O: Print the "vertex_parents" section of the mesh file (ver. >= 1.1).
    void PrintVertexParents(std::ostream &out) const;