Some additions and tweaks to get ex14p to work with empty partitions

fem/fespace.cpp

@@ -3201,6 +3201,16 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
    return FE;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalFE() const
+{
+   if (mesh->GetNE() > 0) { return GetFE(0); }
+
+   Geometry::Type geom = mesh->GetTypicalElementGeometry();
+   const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
+   MFEM_VERIFY(fe != nullptr, "Could not determine a typical FE!");
+   return fe;
+}
+
 const FiniteElement *FiniteElementSpace::GetBE(int i) const
 {
    int order = fec->GetOrder();

fem/fespace.hpp

@@ -1168,6 +1168,13 @@ public:
         an empty partition. */
    virtual const FiniteElement *GetFE(int i) const;
 
+   /** @brief Return GetFE(0) if the local mesh is not empty; otherwise return a
+       typical FE based on the Geometry types in the global mesh.
+
+       This method can be used as a replacement for GetFE(0) that will be valid
+       even if the local mesh is empty. */
+   const FiniteElement *GetTypicalFE() const;
+
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th boundary face in the mesh object. */
    const FiniteElement *GetBE(int i) const;
@@ -1341,18 +1348,19 @@ public:
    virtual ~FiniteElementSpace();
 };
 
-/// @brief Return true if the mesh contains only one topology and the elements are tensor elements.
+/// @brief Return true if the mesh contains only one topology and the elements
+/// are tensor elements.
 inline bool UsesTensorBasis(const FiniteElementSpace& fes)
 {
    Mesh & mesh = *fes.GetMesh();
    const bool mixed = mesh.GetNumGeometries(mesh.Dimension()) > 1;
-   // Potential issue: empty local mesh --> no element 0.
    return !mixed &&
-          dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
+          dynamic_cast<const mfem::TensorBasisElement *>(
+             fes.GetTypicalFE()) != nullptr;
 }
 
-/// @brief Return LEXICOGRAPHIC if mesh contains only one topology and the elements are tensor
-/// elements, otherwise, return NATIVE.
+/// @brief Return LEXICOGRAPHIC if mesh contains only one topology and the
+/// elements are tensor elements, otherwise, return NATIVE.
 ElementDofOrdering GetEVectorOrdering(const FiniteElementSpace& fes);
 
 }

fem/integ/bilininteg_dgdiffusion_pa.cpp

@@ -453,6 +453,7 @@ void DGDiffusionIntegrator::SetupPA(const FiniteElementSpace &fes,
 
    const int ne = fes.GetNE();
    nf = fes.GetNFbyType(type);
+   if (nf == 0) { return; }
 
    // Assumes tensor-product elements
    Mesh &mesh = *fes.GetMesh();

fem/integ/bilininteg_diffusion_pa.cpp

@@ -26,8 +26,7 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
-   if (mesh->GetNE() == 0) { return; }
-   const FiniteElement &el = *fes.GetFE(0);
+   const FiniteElement &el = *fes.GetTypicalFE();
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
    if (DeviceCanUseCeed())
    {

fem/qspace.cpp

@@ -139,9 +139,9 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
 }
 
 QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
-   : QuadratureSpaceBase(mesh_, mesh_.GetElementGeometry(0), ir)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
 {
-   MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) == 1,
+   MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
    ConstructOffsets();
 }

fem/restriction.cpp

@@ -939,7 +939,7 @@ L2FaceRestriction::L2FaceRestriction(const FiniteElementSpace &fes,
      face_dofs(nf > 0 ?
                fes.GetTraceElement(0, fes.GetMesh()->GetFaceGeometry(0))->GetDof()
                : 0),
-     elem_dofs(fes.GetFE(0)->GetDof()),
+     elem_dofs(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
      nfdofs(nf*face_dofs),
      ndofs(fes.GetNDofs()),
      type(type),
@@ -1235,7 +1235,7 @@ void L2FaceRestriction::CheckFESpace()
 
 #ifdef MFEM_DEBUG
    // If fespace == L2
-   const FiniteElement *fe0 = fes.GetFE(0);
+   const FiniteElement *fe0 = fes.GetTypicalFE();
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe0);
    MFEM_VERIFY(tfe != NULL &&
                (tfe->GetBasisType()==BasisType::GaussLobatto ||

mesh/mesh.cpp

@@ -1484,6 +1484,34 @@ Array<int> Mesh::GetFaceToBdrElMap() const
    return face_to_be;
 }
 
+Geometry::Type Mesh::GetTypicalElementGeometry() const
+{
+   if (GetNE() > 0) { return GetElementGeometry(0); }
+
+   const int dim = Dimension();
+   if (dim == 1)
+   {
+      return Geometry::SEGMENT;
+   }
+   const int mesh_gen = MeshGenerator();
+   Geometry::Type geom = Geometry::INVALID;
+   if (dim == 2)
+   {
+      geom = ((mesh_gen & 1) ? Geometry::TRIANGLE :
+              ((mesh_gen & 2) ? Geometry::SQUARE : Geometry::INVALID));
+   }
+   else if (dim == 3)
+   {
+      geom = ((mesh_gen & 1) ? Geometry::TETRAHEDRON :
+              ((mesh_gen & 2) ? Geometry::CUBE :
+               ((mesh_gen & 4) ? Geometry::PRISM :
+                ((mesh_gen & 8) ? Geometry::PYRAMID : Geometry::INVALID))));
+   }
+   MFEM_VERIFY(geom != Geometry::INVALID,
+               "Could not determine a typical element Geometry!");
+   return geom;
+}
+
 void Mesh::Init()
 {
    // in order of declaration:
@@ -13329,7 +13357,7 @@ void GeometricFactors::Compute(const GridFunction &nodes,
 {
 
    const FiniteElementSpace *fespace = nodes.FESpace();
-   const FiniteElement *fe = fespace->GetFE(0);
+   const FiniteElement *fe = fespace->GetTypicalFE();
    const int dim  = fe->GetDim();
    const int vdim = fespace->GetVDim();
    const int NE   = fespace->GetNE();

mesh/mesh.hpp

@@ -1370,6 +1370,13 @@ public:
       return elements[i]->GetGeometryType();
    }
 
+   /** @brief If the local mesh is not empty, return GetElementGeometry(0);
+       otherwise, return a typical Geometry present in the global mesh.
+
+       This method can be used to replace calls like GetElementGeometry(0) in
+       order to handle empty local meshes better. */
+   Geometry::Type GetTypicalElementGeometry() const;
+
    Geometry::Type GetBdrElementGeometry(int i) const
    {
       return boundary[i]->GetGeometryType();