Merge remote-tracking branch 'origin/hypre-runtime-compute-policy' into apps-dev

fem/coefficient.cpp

@@ -1592,6 +1592,10 @@ void VectorQuadratureFunctionCoefficient::Eval(Vector &V,
    QuadF.HostRead();
 
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { V = 0.0; return; }
+
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
 
    if (index == 0 && vdim == QuadF.GetVDim())
@@ -1629,6 +1633,9 @@ double QuadratureFunctionCoefficient::Eval(ElementTransformation &T,
    QuadF.HostRead();
    Vector temp(1);
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { return 0.0; }
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
    QuadF.GetValues(el_idx, ip_idx, temp);
    return temp[0];

fem/integ/bilininteg_mass_pa.cpp

@@ -128,7 +128,7 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
 
    int map_type = el.GetMapType();
    dim = el.GetDim(); // Dimension of the boundary element, *not* the mesh
-   ne = fes.GetMesh()->GetNBE();
+   ne = fes.GetMesh()->GetNFbyType(FaceType::Boundary);
    nq = ir->GetNPoints();
    face_geom = mesh->GetFaceGeometricFactors(*ir, GeometricFactors::DETERMINANTS,
                                              FaceType::Boundary, mt);

fem/qspace.cpp

@@ -164,13 +164,20 @@ int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const
 
 int FaceQuadratureSpace::GetEntityIndex(const ElementTransformation &T) const
 {
+   auto get_face_index = [this](const int idx)
+   {
+      const auto it = face_indices_inv.find(idx);
+      if (it == face_indices_inv.end()) { return -1; }
+      else { return it->second; }
+   };
+
    switch (T.ElementType)
    {
       case ElementTransformation::FACE:
-         return face_indices_inv.at(T.ElementNo);
+         return get_face_index(T.ElementNo);
       case ElementTransformation::BDR_ELEMENT:
       case ElementTransformation::BDR_FACE:
-         return face_indices_inv.at(mesh.GetBdrElementEdgeIndex(T.ElementNo));
+         return get_face_index(mesh.GetBdrElementEdgeIndex(T.ElementNo));
       default:
          MFEM_ABORT("Invalid element type.");
          return -1;

fem/qspace.hpp

@@ -87,7 +87,8 @@ public:
    ///
    /// For a QuadratureSpace defined on elements, this just returns the element
    /// index. For FaceQuadratureSpace, the returned index depends on the chosen
-   /// FaceType.
+   /// FaceType. If the entity is not found (for example, if @a T represents an
+   /// interior face, and the space has FaceType::Boundary) then -1 is returned.
    virtual int GetEntityIndex(const ElementTransformation &T) const = 0;
 
    /// Write the QuadratureSpace to the stream @a out.