Applying potential fix for boundary element index being used on face index call.

fem/bilinearform.cpp

@@ -278,6 +278,13 @@ void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
    boundary_face_integs_marker.Append(&bdr_marker);
 }
 
+void BilinearForm::AddInteriorBoundaryFaceIntegrator(BilinearFormIntegrator* bfi,
+                                        Array<int>& int_bdr_marker)
+{
+    interior_boundary_face_integs.Append(bfi);
+    interior_boundary_face_integs_marker.Append(&int_bdr_marker);
+}
+
 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
 {
    if (element_matrices)
@@ -628,6 +635,62 @@ void BilinearForm::Assemble(int skip_zeros)
             }
          }
       }
+
+      if (interior_boundary_face_integs.Size())
+      {
+          FaceElementTransformations* tr;
+          Array<int> vdofs2;
+          const FiniteElement* fe1, * fe2;
+
+          // Which interior boundary attributes need to be processed?
+          Array<int> int_bdr_attr_marker(mesh->bdr_attributes.Size() ?
+              mesh->bdr_attributes.Max() : 0);
+          int_bdr_attr_marker = 0;
+          for (int k = 0; k < interior_boundary_face_integs.Size(); k++)
+          {
+              if (interior_boundary_face_integs_marker[k] == NULL)
+              {
+                  int_bdr_attr_marker = 1;
+                  break;
+              }
+              Array<int>& int_bdr_marker = *interior_boundary_face_integs_marker[k];
+              MFEM_ASSERT(int_bdr_marker.Size() == int_bdr_attr_marker.Size(),
+                  "invalid boundary marker for boundary face integrator #"
+                  << k << ", counting from zero");
+              for (int i = 0; i < int_bdr_attr_marker.Size(); i++)
+              {
+                  int_bdr_attr_marker[i] |= int_bdr_marker[i];
+              }
+          }
+
+          for (int i = 0; i < fes->GetNBE(); i++)
+          {
+              const int mesh_bdr_attr = mesh->GetBdrAttribute(i);
+              if (int_bdr_attr_marker[mesh_bdr_attr - 1] == 0) { continue; }
+
+              tr = mesh->GetInteriorFaceTransformations(mesh->GetBdrFace(i));
+              if (tr != NULL)
+              {
+                  fes->GetElementVDofs(tr->Elem1No, vdofs);
+                  fes->GetElementVDofs(tr->Elem2No, vdofs2);
+                  vdofs.Append(vdofs2);
+                  fe1 = fes->GetFE(tr->Elem1No);
+                  fe2 = fes->GetFE(tr->Elem2No);
+                  for (int k = 0; k < interior_boundary_face_integs.Size(); k++)
+                  {
+                      if (interior_boundary_face_integs_marker[k] &&
+                          (*interior_boundary_face_integs_marker[k])[mesh_bdr_attr - 1] == 0)
+                      {
+                          continue;
+                      }
+
+                      interior_boundary_face_integs[k]->AssembleFaceMatrix(*fe1, *fe2, *tr,
+                          elemmat);
+                      mat->AddSubMatrix(vdofs, vdofs, elemmat, skip_zeros);
+                  }
+              }
+          }
+   }
    }
 
 #ifdef MFEM_USE_LEGACY_OPENMP

fem/bilinearform.hpp

@@ -113,6 +113,10 @@ protected:
    Array<BilinearFormIntegrator*> boundary_face_integs;
    Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
 
+    /// Set of interior boundary face Integrators to be applied.
+   Array<BilinearFormIntegrator*> interior_boundary_face_integs;
+   Array<Array<int>*> interior_boundary_face_integs_marker; ///< Entries are not owned.
+
    DenseMatrix elemmat;
    Array<int>  vdofs;
 
@@ -272,6 +276,11 @@ public:
    Array<Array<int>*> *GetBFBFI_Marker()
    { return &boundary_face_integs_marker; }
 
+   /// Access all the integrators added with AddInteriorBoundaryFaceIntegrator().
+   Array<BilinearFormIntegrator*>* GetIBFI() { return &interior_boundary_face_integs; }
+   /** @brief Access all boundary markers added with AddInteriorBoundaryFaceIntegrator().*/
+   Array<Array<int>*>* GetIBFI_Marker() { return &interior_boundary_face_integs_marker; }
+
    /// Returns a reference to: \f$ M_{ij} \f$
    const double &operator()(int i, int j) { return (*mat)(i,j); }
 
@@ -416,6 +425,11 @@ public:
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
                              Array<int> &bdr_marker);
 
+   /// Adds new Interior Boundary Integrator restricted to certain faces specified by
+   /// the @a int_bdr_marker.
+   void AddInteriorBoundaryFaceIntegrator(BilinearFormIntegrator* bfi,
+       Array<int>& int_bdr_marker);
+
    /// Sets all sparse values of \f$ M \f$ and \f$ M_e \f$ to 'a'.
    void operator=(const double a)
    {