merge with master and resolve conflicts

fem/bilinearform.cpp

@@ -99,6 +99,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    boundary_integs_marker = bf->boundary_integs_marker;
 
    interior_face_integs = bf->interior_face_integs;
+   interior_face_integs_marker = bf->interior_face_integs_marker;
 
    boundary_face_integs = bf->boundary_face_integs;
    boundary_face_integs_marker = bf->boundary_face_integs_marker;
@@ -254,9 +255,17 @@ void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
    boundary_integs_marker.Append(&bdr_marker);
 }
 
+void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi,
+                                             Array<int> *face_marker)
+{
+   interior_face_integs.Append (bfi);
+   interior_face_integs_marker.Append(face_marker);
+}
+
 void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi)
 {
    interior_face_integs.Append (bfi);
+   interior_face_integs_marker.Append(nullptr);
 }
 
 void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
@@ -685,6 +694,11 @@ void BilinearForm::Assemble(int skip_zeros)
             vdofs.Append (vdofs2);
             for (int k = 0; k < interior_face_integs.Size(); k++)
             {
+               // Skip the face if it's not in the integrator's attribute list.
+               if (interior_face_integs_marker[k] &&
+                   interior_face_integs_marker[k]->Find(tr->Attribute) == -1)
+               { continue; }
+
                interior_face_integs[k]->
                AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
                                   *fes->GetFE(tr->Elem2No),
@@ -1502,6 +1516,12 @@ void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
    trace_face_integs.Append (bfi);
 }
 
+
+void MixedBilinearForm::AddFaceIntegrator (BilinearFormIntegrator *bfi)
+{
+   face_integs.Append(bfi);
+}
+
 void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
 {
    boundary_trace_face_integs.Append(bfi);
@@ -1670,7 +1690,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
    if (boundary_face_integs.Size())
    {
       FaceElementTransformations *ftr;
-      Array<int> tr_vdofs2, te_vdofs2;
+      Array<int> trial_vdofs2, test_vdofs2;
       const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
 
       // Which boundary attributes need to be processed?
@@ -1761,10 +1781,39 @@ void MixedBilinearForm::Assemble(int skip_zeros)
       }
    }
 
+   if (face_integs.Size())
+   {
+      FaceElementTransformations *ftr;
+      Array<int> trial_vdofs2, test_vdofs2;
+
+      int nfaces = mesh->GetNumFaces();
+      for (int i = 0; i < nfaces; i++)
+      {
+         ftr = mesh->GetFaceElementTransformations(i);
+         trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
+         test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
+         if (ftr->Elem2No >= 0)
+         {
+            trial_fes->GetElementVDofs(ftr->Elem2No, trial_vdofs2);
+            test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
+            trial_vdofs.Append(trial_vdofs2);
+            test_vdofs.Append(test_vdofs2);
+         }
+
+         for (int k = 0; k < face_integs.Size(); k++)
+         {
+            face_integs[k]->AssembleFaceMatrix(*trial_fes->GetFE(ftr->Elem1No),
+                                               *test_fes->GetFE(ftr->Elem1No),
+                                               *ftr, elemmat);
+            mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
+         }
+      }
+   }
+
    if (boundary_trace_face_integs.Size())
    {
       FaceElementTransformations *ftr;
-      Array<int> te_vdofs2;
+      Array<int> test_vdofs2;
       const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
 
       // Which boundary attributes need to be processed?
@@ -2355,6 +2404,8 @@ MixedBilinearForm::~MixedBilinearForm()
       { delete boundary_face_integs[i]; }
       for (i = 0; i < trace_face_integs.Size(); i++)
       { delete trace_face_integs[i]; }
+      for (i = 0; i < face_integs.Size(); i++)
+      { delete face_integs[i]; }
       for (i = 0; i < boundary_trace_face_integs.Size(); i++)
       { delete boundary_trace_face_integs[i]; }
    }

fem/bilinearform.hpp

@@ -114,6 +114,10 @@ protected:
 
    /// Set of interior face Integrators to be applied.
    Array<BilinearFormIntegrator*> interior_face_integs;
+   /// List of attributes for each integrator. The integrator is applie only on
+   /// faces that have such attributes. Corresponds to Mesh::GetFaceAttribute().
+   /// Note: it is a list; it's not a marker over all existing face attributes.
+   Array<Array<int>*> interior_face_integs_marker;  ///< Entries are not owned.
 
    /// Set of boundary face Integrators to be applied.
    Array<BilinearFormIntegrator*> boundary_face_integs;
@@ -278,6 +282,7 @@ public:
 
    /// Access all integrators added with AddInteriorFaceIntegrator().
    Array<BilinearFormIntegrator*> *GetFBFI() { return &interior_face_integs; }
+   Array<Array<int>*> *GetFBFI_Marker() { return &interior_face_integs_marker; }
 
    /// Access all integrators added with AddBdrFaceIntegrator().
    Array<BilinearFormIntegrator*> *GetBFBFI() { return &boundary_face_integs; }
@@ -424,6 +429,9 @@ public:
                               Array<int> &bdr_marker);
 
    /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi,
+                                  Array<int> *face_marker);
+
    void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
@@ -796,6 +804,9 @@ protected:
    /// Trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> trace_face_integs;
 
+   /// face integrators.
+   Array<BilinearFormIntegrator*> face_integs;
+
    /// Boundary trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> boundary_trace_face_integs;
    /// Entries are not owned.
@@ -910,6 +921,9 @@ public:
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
    void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);
 
+   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddFaceIntegrator(BilinearFormIntegrator *bfi);
+
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
    void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
                               Array<int> &bdr_marker);

fem/fespace.hpp

@@ -676,6 +676,20 @@ public:
    /// Copy assignment not supported
    FiniteElementSpace& operator=(const FiniteElementSpace&) = delete;
 
+   void ReplaceElemDofTable(const Table &new_elem_dof, int ndofs_new)
+   {
+      *elem_dof = new_elem_dof;
+      ndofs = ndofs_new;
+   }
+   void ReplaceBdrElemDofTable(const Table &new_bdrElem_dof)
+   {
+      *bdr_elem_dof = new_bdrElem_dof;
+   }
+   void ReplaceFaceDofTable(const Table &new_face_dof)
+   {
+      *face_dof = new_face_dof;
+   }
+
    /// Returns the mesh
    inline Mesh *GetMesh() const { return mesh; }
 

fem/linearform.cpp

@@ -105,6 +105,13 @@ void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
    interior_face_integs.Append(lfi);
 }
 
+void LinearForm::AddTraceFaceIntegrator(LinearFormIntegrator *tfi,
+                                        Array<int> &attr_list)
+{
+   trace_face_integs.Append(tfi);
+   trace_face_integs_attributes.Append(&attr_list);
+}
+
 bool LinearForm::SupportsDevice() const
 {
    // return false for NURBS meshes, so we don’t convert it to non-NURBS
@@ -313,6 +320,39 @@ void LinearForm::Assemble()
       }
    }
 
+   if (trace_face_integs.Size())
+   {
+      Mesh *mesh = fes->GetMesh();
+      FaceElementTransformations *tr;
+      const FiniteElement *fe_1, *fe_2;
+      Array<int> vdofs2;
+
+      const int nfaces = mesh->GetNumFaces();
+      for (int f = 0; f < nfaces; f++)
+      {
+         const int attr = mesh->GetFace(f)->GetAttribute();
+         tr   = mesh->GetFaceElementTransformations(f);
+         fe_1 = fes->GetFE(tr->Elem1No);
+         fes->GetElementVDofs(tr->Elem1No, vdofs);
+         if (tr->Elem2No >= 0)
+         {
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+            fe_2 = fes->GetFE(tr->Elem2No);
+         }
+         else
+         {
+            fe_2 = fe_1;
+         }
+
+         for (int k = 0; k < trace_face_integs.Size(); k++)
+         {
+            trace_face_integs[k]->AssembleRHSElementVect(*fe_1, *fe_2, *tr, elemvect);
+            AddElementVector(vdofs, elemvect);
+         }
+      }
+   }
+
    if (interior_face_integs.Size())
    {
       Mesh *mesh = fes->GetMesh();

fem/linearform.hpp

@@ -65,6 +65,11 @@ protected:
    /// Set of Internal Face Integrators to be applied.
    Array<LinearFormIntegrator*> interior_face_integs;
 
+   /// Set of trace (all faces - both interior and boundary) integrators.
+   Array<LinearFormIntegrator *> trace_face_integs;
+   Array<Array<int> *>
+   trace_face_integs_attributes; ///< Entries are not owned.
+
    /// The element ids where the centers of the delta functions lie
    Array<int> domain_delta_integs_elem_id;
 
@@ -151,6 +156,13 @@ public:
    /// Adds new Boundary Face Integrator. Assumes ownership of @a lfi.
    void AddBdrFaceIntegrator(LinearFormIntegrator *lfi);
 
+   /** @brief Add new Trace Face Integrator, restricted to the given face
+       attributes.
+       Assumes ownership of @a tfi. The array @a attr_list is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddTraceFaceIntegrator(LinearFormIntegrator *tfi,
+                               Array<int> &attr_list);
+
    /** @brief Add new Boundary Face Integrator, restricted to the given boundary
        attributes.
 
@@ -182,6 +194,9 @@ public:
    /// Access all integrators added with AddBdrFaceIntegrator().
    Array<LinearFormIntegrator*> *GetFLFI() { return &boundary_face_integs; }
 
+   /// Access all integrators added with AddTraceFaceIntegrator().
+   Array<LinearFormIntegrator*> *GetTLFI() { return &trace_face_integs; }
+
    /// Access all integrators added with AddInteriorFaceIntegrator().
    Array<LinearFormIntegrator*> *GetIFLFI() { return &interior_face_integs; }
 

fem/lor/lor.cpp

@@ -61,6 +61,51 @@ void LORBase::AddIntegratorsAndMarkers(BilinearForm &a_from,
    }
 }
 
+void LORBase::AddIntegratorsAndMarkers2(BilinearForm &a_from,
+                                       BilinearForm &a_to,
+                                       GetIntegratorsFn get_integrators,
+                                       GetMarkersFn get_markers,
+                                       AddIntegratorMarkersFn add_integrator_marker,
+                                       AddIntegratorFn add_integrator,
+                                       const IntegrationRule *ir)
+{
+   Array<BilinearFormIntegrator*> *integrators = (a_from.*get_integrators)();
+   Array<Array<int>*> &markers = *(a_from.*get_markers)();
+
+   for (int i=0; i<integrators->Size(); ++i)
+   {
+      BilinearFormIntegrator *integrator = (*integrators)[i];
+      if (markers[i] != nullptr)
+      {
+         (a_to.*add_integrator_marker)(integrator, *markers[i]);
+      }
+      else
+      {
+         (a_to.*add_integrator)(integrator);
+      }
+      ir_map[integrator] = integrator->GetIntRule();
+      if (ir) { integrator->SetIntegrationRule(*ir); }
+   }
+}
+
+void LORBase::AddIntegratorsAndMarkers(BilinearForm &a_from,
+                                       BilinearForm &a_to,
+                                       GetIntegratorsFn get_integrators,
+                                       GetMarkersFn get_markers,
+                                       AddIntegratorMarkersPtrFn add_integrator,
+                                       const IntegrationRule *ir)
+{
+   Array<BilinearFormIntegrator*> *integrators = (a_from.*get_integrators)();
+   Array<Array<int>*> *markers = (a_from.*get_markers)();
+
+   for (int i=0; i<integrators->Size(); ++i)
+   {
+      (a_to.*add_integrator)((*integrators)[i], *markers[i]);
+      ir_map[(*integrators)[i]] = ((*integrators)[i])->GetIntegrationRule();
+      if (ir) { ((*integrators)[i])->SetIntegrationRule(*ir); }
+   }
+}
+
 void LORBase::ResetIntegrationRules(GetIntegratorsFn get_integrators)
 {
    Array<BilinearFormIntegrator*> *integrators = (a->*get_integrators)();

fem/lor/lor.hpp

@@ -27,6 +27,8 @@ private:
    using AddIntegratorFn = void (BilinearForm::*)(BilinearFormIntegrator*);
    using AddIntegratorMarkersFn =
       void (BilinearForm::*)(BilinearFormIntegrator*, Array<int>&);
+   using AddIntegratorMarkersPtrFn =
+      void (BilinearForm::*)(BilinearFormIntegrator*, Array<int>*);
 
    IntegrationRules irs;
    const IntegrationRule *ir_el, *ir_face;
@@ -52,6 +54,19 @@ private:
                                  AddIntegratorMarkersFn add_integrator_marker,
                                  AddIntegratorFn add_integrator,
                                  const IntegrationRule *ir);
+   void AddIntegratorsAndMarkers2(BilinearForm &a_from,
+                                 BilinearForm &a_to,
+                                 GetIntegratorsFn get_integrators,
+                                 GetMarkersFn get_markers,
+                                 AddIntegratorMarkersFn add_integrator_marker,
+                                 AddIntegratorFn add_integrator,
+                                 const IntegrationRule *ir);
+   void AddIntegratorsAndMarkers(BilinearForm &a_from,
+                                 BilinearForm &a_to,
+                                 GetIntegratorsFn get_integrators,
+                                 GetMarkersFn get_markers,
+                                 AddIntegratorMarkersPtrFn add_integrator,
+                                 const IntegrationRule *ir);
 
    /// Resets the integration rules of the integrators of @a a to their original
    /// values (after temporarily changing them for LOR assembly).

general/table.cpp

@@ -193,6 +193,22 @@ void Table::ShiftUpI()
    I[0] = 0;
 }
 
+void Table::ReplaceConnection(int r, int c_old, int c_new)
+{
+   for (int j=I[r]; j<I[r+1]; j++)
+   {
+      if ( J[j] == c_old) { J[j] = c_new; }
+   }
+}
+
+void Table::RemoveRow(int r)
+{
+   for (int j=I[r]; j<I[r+1]; j++)
+   {
+      J[j] = -1;
+   }
+}
+
 void Table::SetSize(int dim, int connections_per_row)
 {
    SetDims (dim, dim * connections_per_row);

general/table.hpp

@@ -89,6 +89,9 @@ public:
    void AddConnections (int r, const int *c, int nc);
    void ShiftUpI();
 
+   void ReplaceConnection(int r, int c_old, int c_new);
+   void RemoveRow(int r);
+
    /// Set the size and the number of connections for the table.
    void SetSize(int dim, int connections_per_row);
 

mesh/mesh.cpp

@@ -540,6 +540,7 @@ void Mesh::GetFaceTransformation(int FaceNo,
                                  IsoparametricTransformation *FTr) const
 {
    FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
+   FTr->Attribute = faces[FaceNo]->GetAttribute();
    FTr->ElementNo = FaceNo;
    FTr->ElementType = ElementTransformation::FACE;
    FTr->mesh = this;

mesh/mesh.hpp

@@ -2244,6 +2244,12 @@ public:
 
    /// @}
 
+   /// Return the attribute of element i.
+   int GetFaceAttribute(int i) const { return faces[i]->GetAttribute(); }
+
+   /// Set the attribute of face element i.
+   void SetFaceAttribute(int i, int attr) { faces[i]->SetAttribute(attr); }
+
    /// @name Methods related to mesh partitioning
    /// @{
 

miniapps/shifted/cutH1.cpp

@@ -0,0 +1,376 @@
+//            MFEM Example for cutting an H1 space along select faces.
+//
+// Compile with: make cutH1
+//
+// Sample runs:
+//   ./cutH1 -m ../../data/star.mesh -rs 1
+//
+
+#include "mfem.hpp"
+#include "../common/mfem-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace common;
+
+// Used for debugging the elem-to-dof tables when the elements' attributes
+// are associated with materials. Options for lvl:
+// 0 - only duplicated materials per DOF.
+// 1 - all materials per DOF .
+// 2 - full element/material output per DOF.
+void PrintDofElemTable(const Table &elem_dof, const ParMesh &pmesh,
+                       int lvl, bool boundary)
+{
+   Table dof_elem;
+   Transpose(elem_dof, dof_elem);
+
+   const int nrows = dof_elem.Size();
+   if (boundary == false)
+   {
+      std::cout << "--- Dof-to-Elem. Total elem DOFs: " << nrows << std::endl;
+   }
+   else
+   {
+      std::cout << "--- Dof-to-Bdr. Total bndry DOFs: " << nrows << std::endl;
+   }
+
+   Array<int> dof_elements;
+   for (int dof = 0; dof < nrows; dof++)
+   {
+      // Find the materials that share the current dof.
+      std::set<int> dof_materials;
+      dof_elem.GetRow(dof, dof_elements);
+      if (lvl == 2) { std::cout << "Elements for DOF " << dof << ": \n"; }
+      for (int e = 0; e < dof_elements.Size(); e++)
+      {
+         int mat_id;
+         if (boundary == false)
+         {
+            mat_id = pmesh.GetAttribute(dof_elements[e]);
+         }
+         else
+         {
+            int face_id = pmesh.GetBdrFace(dof_elements[e]);
+            int elem_id, tmp;
+            pmesh.GetFaceElements(face_id, &elem_id, &tmp);
+            mat_id = pmesh.GetAttribute(elem_id);
+         }
+
+         if (lvl == 2) { cout << dof_elements[e] << "(" << mat_id << ") "; }
+
+         dof_materials.insert(mat_id);
+      }
+      if (lvl == 2) { std::cout << std::endl; }
+
+      if (lvl == 2) { continue; }
+      if (lvl == 0 && dof_materials.size() < 2) { continue; }
+
+      std::cout << "Materials for DOF " << dof << ": " << std::endl;
+      for (auto it = dof_materials.cbegin(); it != dof_materials.cend(); it++)
+      { std::cout << *it << ' '; }
+      std::cout << std::endl;
+   }
+   std::cout << "--- End of Table" << std::endl;
+}
+
+void VisualizeL2(ParGridFunction &gf, int size, int x, int y)
+{
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+
+   ParMesh *pmesh = gf.ParFESpace()->GetParMesh();
+   const int order = gf.ParFESpace()->GetOrder(0);
+   L2_FECollection fec(order, pmesh->Dimension());
+   ParFiniteElementSpace pfes(pmesh, &fec);
+   ParGridFunction gf_l2(&pfes);
+   gf_l2.ProjectGridFunction(gf);
+
+   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 << gf_l2;
+   sol_sock << "window_geometry " << x << " " << y << " "
+            << size << " " << size << "\n"
+            << "window_title '" << "Y" << "'\n"
+            << "keys mRjlc\n" << flush;
+}
+
+void cutH1Space(ParFiniteElementSpace &pfes, bool vis, bool print)
+{
+   ParMesh &pmesh = *pfes.GetParMesh();
+   ParGridFunction x_vis(&pfes);
+
+   // Duplicate DOFs on the material interface.
+   // That is, the DOF touches different element attributes.
+   const Table &elem_dof = pfes.GetElementToDofTable(),
+                &bdre_dof = pfes.GetBdrElementToDofTable();
+   Table dof_elem, dof_bdre;
+   Table new_elem_dof(elem_dof), new_bdre_dof(bdre_dof);
+   Transpose(elem_dof, dof_elem);
+   Transpose(bdre_dof, dof_bdre);
+   const int nrows = dof_elem.Size(), n_bdr_dofs = dof_bdre.Size();
+   int ndofs = nrows;
+   Array<int> dof_elements, dof_boundaries;
+   if (print)
+   {
+      PrintDofElemTable(elem_dof, pmesh, 0, false);
+      PrintDofElemTable(bdre_dof, pmesh, 2, true);
+   }
+   for (int dof = 0; dof < nrows; dof++)
+   {
+      // Check which materials share the current dof.
+      std::set<int> dof_materials;
+      dof_elem.GetRow(dof, dof_elements);
+      for (int e = 0; e < dof_elements.Size(); e++)
+      {
+         const int mat_id = pmesh.GetAttribute(dof_elements[e]);
+         dof_materials.insert(mat_id);
+      }
+      // Count the materials for the current DOF.
+      const int dof_mat_cnt = dof_materials.size();
+
+      // Duplicate the dof if it is shared between materials.
+      if (dof_mat_cnt > 1)
+      {
+         // The material with the lowest index keeps the old DOF id.
+         // All other materials duplicate the dof.
+         auto mat = dof_materials.cbegin();
+         mat++;
+         while (mat != dof_materials.cend())
+         {
+            // Replace in all elements with material mat.
+            const int new_dof_id = ndofs;
+            for (int e = 0; e < dof_elements.Size(); e++)
+            {
+               if (pmesh.GetAttribute(dof_elements[e]) == *mat)
+               {
+                  if (print)
+                  {
+                     std::cout << "Replacing DOF (for element) : "
+                               << dof << " -> " << new_dof_id
+                               << " in EL " << dof_elements[e] << std::endl;
+                  }
+                  new_elem_dof.ReplaceConnection(dof_elements[e],
+                                                 dof, new_dof_id);
+               }
+            }
+
+            // Replace in all boundary elements with material mat.
+            int dof_bdr_cnt = 0;
+            if (dof < n_bdr_dofs)
+            {
+               dof_bdre.GetRow(dof, dof_boundaries);
+               dof_bdr_cnt = dof_boundaries.Size();
+            }
+            for (int b = 0; b < dof_bdr_cnt; b++)
+            {
+               int face_id = pmesh.GetBdrFace(dof_boundaries[b]);
+               int elem_id, tmp;
+               pmesh.GetFaceElements(face_id, &elem_id, &tmp);
+               if (pmesh.GetAttribute(elem_id) == *mat)
+               {
+                  std::cout << "Replacing DOF (for boundary): "
+                            << dof << " -> " << new_dof_id
+                            << " in BE " << dof_boundaries[b] << std::endl;
+                  new_bdre_dof.ReplaceConnection(dof_boundaries[b],
+                                                 dof, new_dof_id);
+               }
+            }
+
+            // TODO go over faces (in face_dof) that have the replaced dof (the
+            // old id), and check if they have the higher el-attributes on
+            // noth sides. For such faces, the face_dof table should be updated
+            // with the  new_dof_id.
+            // These are faces that touch the interface at a point or an edge.
+
+            ndofs++;
+            mat++;
+         }
+      }
+
+      // Used only for visualization.
+      // Must be visualized before the space update.
+      x_vis(dof) = dof_mat_cnt;
+   }
+
+   // Send the solution by socket to a GLVis server.
+   if (vis)
+   {
+      int size = 500;
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      const int myid = pfes.GetMyRank(), num_procs = pfes.GetNRanks();
+
+      socketstream sol_sock_x(vishost, visport);
+      sol_sock_x << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_x.precision(8);
+      sol_sock_x << "solution\n" << pmesh << x_vis;
+      sol_sock_x << "window_geometry " << 0 << " " << 0 << " "
+                 << size << " " << size << "\n"
+                 << "window_title '" << "X" << "'\n"
+                 << "keys mRjlc\n" << flush;
+   }
+
+   if (print)
+   {
+      PrintDofElemTable(elem_dof, pmesh, 0, false);
+      PrintDofElemTable(new_elem_dof, pmesh, 0, false);
+   }
+
+   // Remove face dofs for cut faces.
+   const Table &face_dof = pfes.GetFaceToDofTable();
+   Table new_face_dof(face_dof);
+   for (int f = 0; f < pmesh.GetNumFaces(); f++)
+   {
+      auto *ftr = pmesh.GetFaceElementTransformations(f, 3);
+      if (ftr->Elem2No > 0 &&
+          pmesh.GetAttribute(ftr->Elem1No) != pmesh.GetAttribute(ftr->Elem2No))
+      {
+         if (print)
+         {
+            std::cout << ftr->Elem1No << " " << ftr->Elem2No << std::endl;
+            std::cout << pmesh.GetAttribute(ftr->Elem1No) << " "
+                      << pmesh.GetAttribute(ftr->Elem2No) << std::endl;
+            std::cout << "Removing face dofs for face " << f << std::endl;
+         }
+         new_face_dof.RemoveRow(f);
+      }
+   }
+   new_face_dof.Finalize();
+
+   // Cut the space.
+   pfes.ReplaceElemDofTable(new_elem_dof, ndofs);
+   pfes.ReplaceBdrElemDofTable(new_bdre_dof);
+   pfes.ReplaceFaceDofTable(new_face_dof);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   Mpi::Init(argc, argv);
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int rs_levels = 0;
+   int order = 2;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   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())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0) { args.PrintOptions(cout); }
+
+   // 3. 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);
+   if (myid == 0) { device.Print(); }
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace pfes(&pmesh, &fec);
+
+   // Assign material indices to the element attributes.
+   const int NE = pmesh.GetNE();
+   for (int i = 0; i < NE; i++)
+   {
+      Vector center;
+      pmesh.GetElementCenter(i, center);
+      Element *el = pmesh.GetElement(i);
+      if (center(0) <= 0.5 || center(1) >= 0.5)
+      {
+         el->SetAttribute(0);
+      }
+      else { el->SetAttribute(1); }
+   }
+
+   cutH1Space(pfes, true, true);
+
+   // Set face_attribute = 77 to faces that are on the material interface.
+   // Remove face dofs for cut faces.
+   for (int f = 0; f < pmesh.GetNumFaces(); f++)
+   {
+      auto *ftr = pmesh.GetFaceElementTransformations(f, 3);
+      if (ftr->Elem2No > 0 &&
+          pmesh.GetAttribute(ftr->Elem1No) != pmesh.GetAttribute(ftr->Elem2No))
+      {
+         pmesh.SetFaceAttribute(f, 77);
+      }
+   }
+
+   // Simple Dirichlet BC.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      pfes.FiniteElementSpace::GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // RHS.
+   ParLinearForm b(&pfes);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // LHS.
+   BilinearForm a(&pfes);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   Array<int> cut_face_attributes(1);
+   cut_face_attributes[0] = 77;
+   const double sigma = -1.0, kappa = -1.0;
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                               &cut_face_attributes);
+   a.Assemble();
+
+   // Form the system.
+   ParGridFunction u(&pfes);
+   u = 0.0;
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
+
+   // Solve.
+   GSSmoother M((SparseMatrix&)(*A));
+   //PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   CG(*A, B, X, 1, 1000, 1e-12, 0.0);
+   a.RecoverFEMSolution(X, b, u);
+
+   VisualizeL2(u, 500, 500, 0);
+
+   const double norm = u.Norml2();
+   std::cout << "Norm: " << norm << std::endl;
+   return 0;
+}

miniapps/shifted/makefile

@@ -32,9 +32,11 @@ EXTRAPOLATE_SRC = extrapolate.cpp extrapolator.cpp marking.cpp
 EXTRAPOLATE_OBJ = $(EXTRAPOLATE_SRC:.cpp=.o)
 ALGOIM_SRC =  lsf_integral.cpp
 ALGOIM_OBJ = $(ALGOIM_SRC:.cpp=.o)
+CUTH1_SRC  = cutH1.cpp
+CUTH1_OBJ = $(CUTH1_SRC:.cpp=.o)
 
 SEQ_MINIAPPS = lsf_integral
-PAR_MINIAPPS = distance diffusion extrapolate
+PAR_MINIAPPS = distance diffusion extrapolate cutH1
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -75,6 +77,9 @@ diffusion: $(DIFFUSION_OBJ)
 extrapolate: $(EXTRAPOLATE_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(EXTRAPOLATE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
 
+cutH1: $(CUTH1_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CUTH1_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
 # Rule for building lib-common
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
@@ -99,7 +104,7 @@ $(MFEM_LIB_FILE):
 clean: clean-build clean-exec
 
 clean-build:
-	rm -f *.o *~ distance diffusion extrapolate lsf_integral
+	rm -f *.o *~ distance diffusion extrapolate lsf_integral cutH1
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec: