Better support for PA with empty mesh partitions

examples/ex9.cpp

@@ -77,7 +77,7 @@ class DG_Solver : public Solver
    DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
       : M(M_),
         K(K_),
-        prec(fes.GetFE(0)->GetDof(),
+        prec(fes.GetTypicalFE()->GetDof(),
              BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {

examples/ex9p.cpp

@@ -145,7 +145,7 @@ class DG_Solver : public Solver
         linear_solver(M.GetComm()),
         dt(-1.0)
    {
-      int block_size = fes.GetFE(0)->GetDof();
+      int block_size = fes.GetTypicalFE()->GetDof();
       if (prec_type == PrecType::ILU)
       {
          prec = new BlockILU(block_size,

examples/sundials/ex9.cpp

@@ -74,7 +74,7 @@ class DG_Solver : public Solver
    DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
       : M(M_),
         K(K_),
-        prec(fes.GetFE(0)->GetDof(),
+        prec(fes.GetTypicalFE()->GetDof(),
              BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {

examples/sundials/ex9p.cpp

@@ -141,7 +141,7 @@ class DG_Solver : public Solver
         linear_solver(M.GetComm()),
         dt(-1.0)
    {
-      int block_size = fes.GetFE(0)->GetDof();
+      int block_size = fes.GetTypicalFE()->GetDof();
       if (prec_type == PrecType::ILU)
       {
          prec = new BlockILU(block_size,

fem/bilinearform.cpp

@@ -931,7 +931,7 @@ void BilinearForm::ComputeElementMatrices()
    }
 
    int num_elements = fes->GetNE();
-   int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
+   int num_dofs_per_el = fes->GetTypicalFE()->GetDof() * fes->GetVDim();
 
    element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
                                       num_elements);

fem/bilinearform_ext.cpp

@@ -862,7 +862,7 @@ void EABilinearFormExtension::Assemble()
    SetupRestrictionOperators(L2FaceValues::SingleValued);
 
    ne = trial_fes->GetMesh()->GetNE();
-   elemDofs = trial_fes->GetFE(0)->GetDof();
+   elemDofs = trial_fes->GetTypicalFE()->GetDof();
 
    ea_data.SetSize(ne*elemDofs*elemDofs, Device::GetMemoryType());
    ea_data.UseDevice(true);
@@ -878,9 +878,7 @@ void EABilinearFormExtension::Assemble()
       integrators[i]->AssembleEA(*a->FESpace(), ea_data, i);
    }
 
-   faceDofs = trial_fes ->
-              GetTraceElement(0, trial_fes->GetMesh()->GetFaceGeometry(0)) ->
-              GetDof();
+   faceDofs = trial_fes->GetTypicalTraceElement()->GetDof();
 
    MFEM_VERIFY(a->GetBBFI()->Size() == 0,
                "Element assembly does not support AddBoundaryIntegrator yet.");

fem/bilininteg.cpp

@@ -3474,7 +3474,7 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
    {
       const int order = (ndof2) ? max(el1.GetOrder(),
                                       el2.GetOrder()) : el1.GetOrder();
-      ir = &GetRule(order, Trans);
+      ir = &GetRule(order, Trans.GetGeometryType());
    }
 
    // assemble: < {(Q \nabla u).n},[v] >      --> elmat
@@ -3653,11 +3653,11 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
 }
 
 const IntegrationRule &DGDiffusionIntegrator::GetRule(
-   int order, FaceElementTransformations &T)
+   int order, Geometry::Type geom)
 {
    // order is typically the maximum of the order of the left and right elements
    // neighboring the given face.
-   return IntRules.Get(T.GetGeometryType(), 2*order);
+   return IntRules.Get(geom, 2*order);
 }
 
 // static method

fem/bilininteg.hpp

@@ -3293,7 +3293,7 @@ class DGDiffusionIntegrator : public BilinearFormIntegrator
    void AddMultPAFaceNormalDerivatives(const Vector &x, const Vector &dxdn,
                                        Vector &y, Vector &dydn) const override;
 
-   const IntegrationRule &GetRule(int order, FaceElementTransformations &T);
+   const IntegrationRule &GetRule(int order, Geometry::Type geom);
 
 private:
    void SetupPA(const FiniteElementSpace &fes, FaceType type);

fem/ceed/interface/basis.cpp

@@ -135,7 +135,7 @@ void InitBasis(const FiniteElementSpace &fes,
                const IntegrationRule &ir,
                Ceed ceed, CeedBasis *basis)
 {
-   const mfem::FiniteElement &fe = *fes.GetFE(0);
+   const mfem::FiniteElement &fe = *fes.GetTypicalFE();
    InitBasisImpl(fes, fe, ir, ceed, basis);
 }
 

fem/ceed/interface/restriction.cpp

@@ -23,7 +23,7 @@ namespace ceed
 static void InitNativeRestr(const mfem::FiniteElementSpace &fes,
                             Ceed ceed, CeedElemRestriction *restr)
 {
-   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const mfem::FiniteElement *fe = fes.GetTypicalFE();
    const int P = fe->GetDof();
    CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
    const mfem::Table &el_dof = fes.GetElementToDofTable();
@@ -51,7 +51,7 @@ static void InitNativeRestr(const mfem::FiniteElementSpace &fes,
 static void InitLexicoRestr(const mfem::FiniteElementSpace &fes,
                             Ceed ceed, CeedElemRestriction *restr)
 {
-   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const mfem::FiniteElement *fe = fes.GetTypicalFE();
    const int P = fe->GetDof();
    CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
    const mfem::Table &el_dof = fes.GetElementToDofTable();
@@ -75,7 +75,7 @@ static void InitLexicoRestr(const mfem::FiniteElementSpace &fes,
 static void InitRestrictionImpl(const mfem::FiniteElementSpace &fes,
                                 Ceed ceed, CeedElemRestriction *restr)
 {
-   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const mfem::FiniteElement *fe = fes.GetTypicalFE();
    const mfem::TensorBasisElement * tfe =
       dynamic_cast<const mfem::TensorBasisElement *>(fe);
    if ( tfe && tfe->GetDofMap().Size()>0 ) // Native ordering using dof_map
@@ -225,7 +225,7 @@ void InitRestriction(const FiniteElementSpace &fes,
                      CeedElemRestriction *restr)
 {
    // Check for FES -> basis, restriction in hash tables
-   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const mfem::FiniteElement *fe = fes.GetTypicalFE();
    const int P = fe->GetDof();
    const int nelem = fes.GetNE();
    const int ncomp = fes.GetVDim();

fem/conduitdatacollection.cpp

@@ -719,7 +719,7 @@ ConduitDataCollection::MeshToBlueprintMesh(Mesh *mesh,
    // copy out. Some other cases (sidre) may actually have contig
    // allocation but I am  not sure how to detect this case from mfem
    int num_ele = mesh->GetNE();
-   int geom = mesh->GetElementBaseGeometry(0);
+   int geom = mesh->GetTypicalElementGeometry();
    int idxs_per_ele = Geometry::NumVerts[geom];
    int num_conn_idxs =  num_ele * idxs_per_ele;
 

fem/dgmassinv.cpp

@@ -24,7 +24,7 @@ DGMassInverse::DGMassInverse(FiniteElementSpace &fes_orig, Coefficient *coeff,
      fec(fes_orig.GetMaxElementOrder(),
          fes_orig.GetMesh()->Dimension(),
          btype,
-         fes_orig.GetFE(0)->GetMapType()),
+         fes_orig.GetTypicalFE()->GetMapType()),
      fes(fes_orig.GetMesh(), &fec)
 {
    MFEM_VERIFY(fes.IsDGSpace(), "Space must be DG.");
@@ -42,7 +42,9 @@ DGMassInverse::DGMassInverse(FiniteElementSpace &fes_orig, Coefficient *coeff,
    {
       // original basis to solver basis
       const auto mode = DofToQuad::TENSOR;
-      d2q = &fes_orig.GetFE(0)->GetDofToQuad(fes.GetFE(0)->GetNodes(), mode);
+      const FiniteElement &fe_orig = *fes_orig.GetTypicalFE();
+      const FiniteElement &fe = *fes.GetTypicalFE();
+      d2q = &fe_orig.GetDofToQuad(fe.GetNodes(), mode);
 
       int n = d2q->ndof;
       Array<real_t> B_inv = d2q->B; // deep copy

fem/fe_coll.hpp

@@ -195,6 +195,17 @@ class FiniteElementCollection
       return var_orders[p]->FiniteElementForGeometry(geom);
    }
 
+   /// Variable order version of TraceFiniteElementForGeometry().
+   /** The order parameter @a p represents the order of the highest-dimensional
+       FiniteElement%s the fixed-order collection we want to query. In general,
+       this order is different from the order of the returned FiniteElement. */
+   const FiniteElement *GetTraceFE(Geometry::Type geom, int p) const
+   {
+      if (p == base_p) { return TraceFiniteElementForGeometry(geom); }
+      if (p >= var_orders.Size() || !var_orders[p]) { InitVarOrder(p); }
+      return var_orders[p]->TraceFiniteElementForGeometry(geom);
+   }
+
    /// Variable order version of DofForGeometry().
    /** The order parameter @a p represents the order of the highest-dimensional
        FiniteElement%s the fixed-order collection we want to query. In general,

fem/fespace.cpp

@@ -2454,7 +2454,7 @@ void FiniteElementSpace::Construct()
       else
       {
          // the simple case: all faces are of the same geometry and order
-         uni_fdof = fec->GetNumDof(mesh->GetFaceGeometry(0), order);
+         uni_fdof = fec->GetNumDof(mesh->GetTypicalFaceGeometry(), order);
          nfdofs = mesh->GetNFaces() * uni_fdof;
          var_face_dofs.Clear(); // ensure any old var_face_dof table is dumped.
       }
@@ -3134,7 +3134,7 @@ void FiniteElementSpace::GetFaceInteriorDofs(int i, Array<int> &dofs) const
    }
    else
    {
-      auto geom = mesh->GetFaceGeometry(0);
+      auto geom = mesh->GetTypicalFaceGeometry();
       nf = fec->GetNumDof(geom, fec->GetOrder());
       base = i*nf;
    }
@@ -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();
@@ -3276,7 +3286,12 @@ const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
 const FiniteElement *FiniteElementSpace::GetTraceElement(
    int i, Geometry::Type geom_type) const
 {
-   return fec->TraceFiniteElementForGeometry(geom_type);
+   return fec->GetTraceFE(geom_type, GetElementOrder(i));
+}
+
+const FiniteElement *FiniteElementSpace::GetTypicalTraceElement() const
+{
+   return fec->TraceFiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
 }
 
 FiniteElementSpace::~FiniteElementSpace()

fem/fespace.hpp

@@ -1168,6 +1168,13 @@ class FiniteElementSpace
         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;
@@ -1185,6 +1192,12 @@ class FiniteElementSpace
    /// Return the trace element from element 'i' to the given 'geom_type'
    const FiniteElement *GetTraceElement(int i, Geometry::Type geom_type) const;
 
+   /// @brief Return a "typical" trace element.
+   ///
+   /// This can be used in situations where the local mesh partition may be
+   /// empty.
+   const FiniteElement *GetTypicalTraceElement() const;
+
    /** @brief Mark degrees of freedom associated with boundary elements with
        the specified boundary attributes (marked in 'bdr_attr_is_ess').
        For spaces with 'vdim' > 1, the 'component' parameter can be used
@@ -1341,18 +1354,19 @@ class FiniteElementSpace
    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/gridfunc.cpp

@@ -322,19 +322,7 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
 
 int GridFunction::VectorDim() const
 {
-   const FiniteElement *fe;
-   if (!fes->GetNE())
-   {
-      const FiniteElementCollection *fe_coll = fes->FEColl();
-      static const Geometry::Type geoms[3] =
-      { Geometry::SEGMENT, Geometry::TRIANGLE, Geometry::TETRAHEDRON };
-      fe = fe_coll->
-           FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
-   }
-   else
-   {
-      fe = fes->GetFE(0);
-   }
+   const FiniteElement *fe = fes->GetTypicalFE();
    if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
    {
       return fes->GetVDim();
@@ -345,17 +333,7 @@ int GridFunction::VectorDim() const
 
 int GridFunction::CurlDim() const
 {
-   const FiniteElement *fe;
-   if (!fes->GetNE())
-   {
-      static const Geometry::Type geoms[3] =
-      { Geometry::SEGMENT, Geometry::TRIANGLE, Geometry::TETRAHEDRON };
-      fe = fec->FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
-   }
-   else
-   {
-      fe = fes->GetFE(0);
-   }
+   const FiniteElement *fe = fes->GetTypicalFE();
    if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
    {
       return 2 * fes->GetMesh()->SpaceDimension() - 3;
@@ -1790,8 +1768,8 @@ void GridFunction::ProjectGridFunction(const GridFunction &src)
    {
       // Assuming that the projection matrix is the same for all elements
       sameP = true;
-      fes->GetFE(0)->Project(*src.fes->GetFE(0),
-                             *mesh->GetElementTransformation(0), P);
+      fes->GetTypicalFE()->Project(*src.fes->GetTypicalFE(),
+                                   *mesh->GetTypicalElementTransformation(), P);
    }
    const int vdim = fes->GetVDim();
    MFEM_VERIFY(vdim == src.fes->GetVDim(), "incompatible vector dimensions!");

fem/gslib.cpp

@@ -541,12 +541,12 @@ void FindPointsGSLIB::SetupIntegrationRuleForSplitMesh(Mesh *meshin,
    meshin->SetNodalFESpace(&nodal_fes);
    const int NEsplit = meshin->GetNE();
 
-   const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
+   const int dof_cnt = nodal_fes.GetTypicalFE()->GetDof(),
              pts_cnt = NEsplit * dof_cnt;
    Vector irlist(dim * pts_cnt);
 
    const TensorBasisElement *tbe =
-      dynamic_cast<const TensorBasisElement *>(nodal_fes.GetFE(0));
+      dynamic_cast<const TensorBasisElement *>(nodal_fes.GetTypicalFE());
    MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
    const Array<int> &dof_map = tbe->GetDofMap();
 
@@ -1182,7 +1182,7 @@ void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
    crystal_init(cr, gsl_comm);
    mesh = &m;
    dim  = mesh->Dimension();
-   const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
+   const FiniteElement *fe = mesh->GetNodalFESpace()->GetTypicalFE();
    unsigned dof1D = fe->GetOrder() + 1;
 
    SetupSplitMeshes();

fem/integ/bilininteg_convection_mf.cpp

@@ -21,8 +21,8 @@ void ConvectionIntegrator::AssembleMF(const FiniteElementSpace &fes)
    // Assuming the same element type
    Mesh *mesh = fes.GetMesh();
    if (mesh->GetNE() == 0) { return; }
-   const FiniteElement &el = *fes.GetFE(0);
-   ElementTransformation &Trans = *fes.GetElementTransformation(0);
+   const FiniteElement &el = *fes.GetTypicalFE();
+   ElementTransformation &Trans = *mesh->GetTypicalElementTransformation();
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, Trans);
    if (DeviceCanUseCeed())
    {

fem/integ/bilininteg_convection_pa.cpp

@@ -138,8 +138,8 @@ void ConvectionIntegrator::AssemblePA(const FiniteElementSpace &fes)
                          Device::GetDeviceMemoryType() : pa_mt;
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
-   const FiniteElement &el = *fes.GetFE(0);
-   ElementTransformation &Trans = *fes.GetElementTransformation(0);
+   const FiniteElement &el = *fes.GetTypicalFE();
+   ElementTransformation &Trans = *mesh->GetTypicalElementTransformation();
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, Trans);
    if (DeviceCanUseCeed())
    {