Move face map logic from restriction.cpp into FiniteElement derived classes

Use inheritance instead of switch statement to dispatch

CHANGELOG

@@ -18,7 +18,10 @@ Meshing improvements
 
 Discretization improvements
 ---------------------------
-- TBD
+- Face restriction operators for Nedelec and Raviart-Thomas finite element
+  spaces are now supported through the H1_ND_RT_FaceRestriction class.
+
+- VectorFEBoundaryFluxLFIntegrator is now supported on device/GPU.
 
 Linear and nonlinear solvers
 ----------------------------

fem/CMakeLists.txt

@@ -44,6 +44,7 @@ set(SRCS
   eltrans.cpp
   estimators.cpp
   fe.cpp
+  fe/face_map_utils.cpp
   fe/fe_base.cpp
   fe/fe_fixed_order.cpp
   fe/fe_h1.cpp
@@ -74,6 +75,7 @@ set(SRCS
   linearform_ext.cpp
   lininteg.cpp
   lininteg_boundary.cpp
+  lininteg_boundary_flux.cpp
   lininteg_domain.cpp
   lininteg_domain_grad.cpp
   lor/lor.cpp
@@ -151,6 +153,7 @@ set(HDRS
   eltrans.hpp
   estimators.hpp
   fe.hpp
+  fe/face_map_utils.hpp
   fe/fe_base.hpp
   fe/fe_fixed_order.hpp
   fe/fe_h1.hpp

fem/fe/face_map_utils.cpp

@@ -0,0 +1,60 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// Finite Element Base classes
+
+#include "face_map_utils.hpp"
+
+namespace mfem
+{
+
+std::pair<int,int> GetFaceNormal3D(const int face_id)
+{
+   switch (face_id)
+   {
+      case 0: return std::make_pair(2, 0); // z = 0
+      case 1: return std::make_pair(1, 0); // y = 0
+      case 2: return std::make_pair(0, 1); // x = 1
+      case 3: return std::make_pair(1, 1); // y = 1
+      case 4: return std::make_pair(0, 0); // x = 0
+      case 5: return std::make_pair(2, 1); // z = 1
+      default: MFEM_ABORT("Invalid face ID.")
+   }
+   return std::make_pair(-1, -1); // invalid
+}
+
+void FillFaceMap(const int n_face_dofs_per_component,
+                 const std::vector<int> offsets,
+                 const std::vector<int> &strides,
+                 const std::vector<int> &n_dofs_per_dim,
+                 Array<int> &face_map)
+{
+   const int n_components = offsets.size();
+   const int face_dim = strides.size() / n_components;
+   for (int comp = 0; comp < n_components; ++comp)
+   {
+      const int offset = offsets[comp];
+      for (int i = 0; i < n_face_dofs_per_component; ++i)
+      {
+         int idx = offset;
+         int j = i;
+         for (int d = 0; d < face_dim; ++d)
+         {
+            const int dof1d = n_dofs_per_dim[comp*(face_dim) + d];
+            idx += strides[comp*(face_dim) + d]*(j % dof1d);
+            j /= dof1d;
+         }
+         face_map[comp*n_face_dofs_per_component + i] = idx;
+      }
+   }
+}
+
+}

fem/fe/face_map_utils.hpp

@@ -0,0 +1,49 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_FACE_MAP_UTILS_HPP
+#define MFEM_FACE_MAP_UTILS_HPP
+
+#include "../../general/array.hpp"
+#include <utility> // std::pair
+#include <vector>
+
+namespace mfem
+{
+
+/// Each face of a hexahedron is given by a level set x_i = l, where x_i is one
+/// of x, y, or z (corresponding to i = 0, i=1, i = 3), and l is either 0 or 1.
+/// Returns i and level.
+std::pair<int,int> GetFaceNormal3D(const int face_id);
+
+/// @brief Fills in the entries of the lexicographic face_map.
+///
+/// For use in FiniteElement::GetFaceMap.
+///
+/// n_face_dofs_per_component is the number of DOFs for each vector component
+/// on the face (there is only one vector component in all cases except for 3D
+/// Nedelec elements, where the face DOFs have two components to span the
+/// tangent space).
+///
+/// The DOFs for the i-th vector component begin at offsets[i] (i.e. the number
+/// of vector components is given by offsets.size()).
+///
+/// The DOFs for each vector component are arranged in a Cartesian grid defined
+/// by strides and n_dofs_per_dim.
+void FillFaceMap(const int n_face_dofs_per_component,
+                 const std::vector<int> offsets,
+                 const std::vector<int> &strides,
+                 const std::vector<int> &n_dofs_per_dim,
+                 Array<int> &face_map);
+
+} // namespace mfem
+
+#endif

fem/fe/fe_base.cpp

@@ -12,6 +12,7 @@
 // Finite Element Base classes
 
 #include "fe_base.hpp"
+#include "face_map_utils.hpp"
 #include "../coefficient.hpp"
 
 namespace mfem
@@ -370,6 +371,12 @@ const DofToQuad &FiniteElement::GetDofToQuad(const IntegrationRule &,
    return *dof2quad_array[0]; // suppress a warning
 }
 
+void FiniteElement::GetFaceMap(const int face_id,
+                               Array<int> &face_map) const
+{
+   MFEM_ABORT("method is not implemented for this element");
+}
+
 FiniteElement::~FiniteElement()
 {
    for (int i = 0; i < dof2quad_array.Size(); i++)
@@ -2427,6 +2434,55 @@ void NodalTensorFiniteElement::SetMapType(const int map_type)
    }
 }
 
+void NodalTensorFiniteElement::GetFaceMap(const int face_id,
+                                          Array<int> &face_map) const
+{
+   const int dof1d = order + 1;
+   int n_face_dofs = pow(dof1d, dim - 1);
+   std::vector<int> offsets, strides;
+   switch (dim)
+   {
+      case 1:
+         offsets = {(face_id == 0) ? 0 : dof1d - 1};
+         break;
+      case 2:
+         strides = {(face_id == 0 || face_id == 2) ? 1 : dof1d};
+         switch (face_id)
+         {
+            case 0: offsets = {0}; break; // y = 0
+            case 1: offsets = {dof1d - 1}; break; // x = 1
+            case 2: offsets = {(dof1d-1)*dof1d}; break; // y = 1
+            case 3: offsets = {0}; break; // x = 0
+         }
+         break;
+      case 3:
+      {
+         const auto f = GetFaceNormal3D(face_id);
+         const int face_normal = f.first, level = f.second;
+         if (face_normal == 0) // x-normal
+         {
+            offsets = {level ? dof1d-1 : 0};
+            strides = {dof1d, dof1d*dof1d};
+         }
+         else if (face_normal == 1) // y-normal
+         {
+            offsets = {level ? (dof1d-1)*dof1d : 0};
+            strides = {1, dof1d*dof1d};
+         }
+         else if (face_normal == 2) // z-normal
+         {
+            offsets = {level ? (dof1d-1)*dof1d*dof1d : 0};
+            strides = {1, dof1d};
+         }
+         break;
+      }
+   }
+
+   // same number of DOFs in each dimension, repeat dof1d (dim - 1) times
+   std::vector<int> n_dofs(dim - 1, dof1d);
+   FillFaceMap(n_face_dofs, offsets, strides, n_dofs, face_map);
+}
+
 VectorTensorFiniteElement::VectorTensorFiniteElement(const int dims,
                                                      const int d,
                                                      const int p,

fem/fe/fe_base.hpp

@@ -579,6 +579,10 @@ public:
    /** See the documentation for DofToQuad for more details. */
    virtual const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
                                          DofToQuad::Mode mode) const;
+
+   virtual void GetFaceMap(const int face_id,
+                           Array<int> &face_map) const;
+
    /// Deconstruct the FiniteElement
    virtual ~FiniteElement();
 
@@ -1267,6 +1271,8 @@ public:
          NodalFiniteElement::GetTransferMatrix(fe, Trans, I);
       }
    }
+
+   virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
 };
 
 class VectorTensorFiniteElement : public VectorFiniteElement,

fem/fe/fe_nd.cpp

@@ -12,6 +12,7 @@
 // Nedelec Finite Element classes
 
 #include "fe_nd.hpp"
+#include "face_map_utils.hpp"
 #include "../coefficient.hpp"
 
 namespace mfem
@@ -481,6 +482,50 @@ void ND_HexahedronElement::CalcCurlShape(const IntegrationPoint &ip,
          }
 }
 
+void ND_HexahedronElement::GetFaceMap(const int face_id,
+                                      Array<int> &face_map) const
+{
+   const int p = order;
+   const int pp1 = p + 1;
+   const int n_face_dofs_per_component = p*pp1;
+   const int n_dof_per_dim = p*pp1*pp1;
+
+   std::vector<int> n_dofs = {p, pp1, pp1, p};
+   std::vector<int> offsets, strides;
+
+   const auto f = GetFaceNormal3D(face_id);
+   const int face_normal = f.first, level = f.second;
+   if (face_normal == 0) // x-normal
+   {
+      offsets =
+      {
+         n_dof_per_dim + (level ? pp1 - 1 : 0),
+         2*n_dof_per_dim + (level ? pp1 - 1 : 0)
+      };
+      strides = {pp1, p*pp1, pp1, pp1*pp1};
+   }
+   else if (face_normal == 1) // y-normal
+   {
+      offsets =
+      {
+         level ? p*(pp1 - 1) : 0,
+         2*n_dof_per_dim + (level ? pp1*(pp1 - 1) : 0)
+      };
+      strides = {1, p*pp1, 1, pp1*pp1};
+   }
+   else if (face_normal == 2) // z-normal
+   {
+      offsets =
+      {
+         level ? p*pp1*(pp1 - 1) : 0,
+         n_dof_per_dim + (level ? p*pp1*(pp1 - 1) : 0)
+      };
+      strides = {1, p, 1, pp1};
+   }
+
+   FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs, face_map);
+}
+
 const double ND_QuadrilateralElement::tk[8] =
 { 1.,0.,  0.,1., -1.,0., 0.,-1. };
 
@@ -771,6 +816,25 @@ void ND_QuadrilateralElement::CalcCurlShape(const IntegrationPoint &ip,
       }
 }
 
+void ND_QuadrilateralElement::GetFaceMap(const int face_id,
+                                         Array<int> &face_map) const
+{
+   const int p = order;
+   const int pp1 = order + 1;
+   const int n_face_dofs_per_component = p;
+   std::vector<int> strides = {(face_id == 0 || face_id == 2) ? 1 : pp1};
+   std::vector<int> n_dofs = {p};
+   std::vector<int> offsets;
+   switch (face_id)
+   {
+      case 0: offsets = {0}; break; // y = 0
+      case 1: offsets = {p*pp1 + pp1 - 1}; break; // x = 1
+      case 2: offsets = {p*(pp1 - 1)}; break; // y = 1
+      case 3: offsets = {p*pp1}; break; // x = 0
+   }
+   FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs, face_map);
+}
+
 
 const double ND_TetrahedronElement::tk[18] =
 { 1.,0.,0.,  0.,1.,0.,  0.,0.,1.,  -1.,1.,0.,  -1.,0.,1.,  0.,-1.,1. };

fem/fe/fe_nd.hpp

@@ -91,6 +91,8 @@ public:
                             DenseMatrix &curl) const
    { ProjectCurl_ND(tk, dof2tk, fe, Trans, curl); }
 
+   virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
+
 protected:
    void ProjectIntegrated(VectorCoefficient &vc,
                           ElementTransformation &Trans,
@@ -155,6 +157,8 @@ public:
                             DenseMatrix &grad) const
    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
 
+   virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
+
 protected:
    void ProjectIntegrated(VectorCoefficient &vc,
                           ElementTransformation &Trans,

fem/fe/fe_rt.cpp

@@ -12,6 +12,7 @@
 // Raviart-Thomas Finite Element classes
 
 #include "fe_rt.hpp"
+#include "face_map_utils.hpp"
 #include "../coefficient.hpp"
 
 namespace mfem
@@ -297,6 +298,27 @@ void RT_QuadrilateralElement::ProjectIntegrated(VectorCoefficient &vc,
    }
 }
 
+void RT_QuadrilateralElement::GetFaceMap(const int face_id,
+                                         Array<int> &face_map) const
+{
+   const int p = order;
+   const int pp1 = p + 1;
+   const int n_face_dofs = p;
+
+   std::vector<int> offsets;
+   std::vector<int> strides = {(face_id == 0 || face_id == 2) ? 1 : pp1};
+   switch (face_id)
+   {
+      case 0: offsets = {p*pp1}; break; // y = 0
+      case 1: offsets = {pp1 - 1}; break; // x = 1
+      case 2: offsets = {p*pp1 + p*(pp1 - 1)}; break; // y = 1
+      case 3: offsets = {0}; break; // x = 0
+   }
+
+   std::vector<int> n_dofs(dim - 1, p);
+   FillFaceMap(n_face_dofs, offsets, strides, n_dofs, face_map);
+}
+
 
 const double RT_HexahedronElement::nk[18] =
 { 0.,0.,-1.,  0.,-1.,0.,  1.,0.,0.,  0.,1.,0.,  -1.,0.,0.,  0.,0.,1. };
@@ -686,6 +708,35 @@ void RT_HexahedronElement::ProjectIntegrated(VectorCoefficient &vc,
    }
 }
 
+void RT_HexahedronElement::GetFaceMap(const int face_id,
+                                      Array<int> &face_map) const
+{
+   const int p = order;
+   const int pp1 = p + 1;
+   int n_face_dofs = p*p;
+   std::vector<int> strides, offsets;
+   const int n_dof_per_dim = p*p*pp1;
+   const auto f = GetFaceNormal3D(face_id);
+   const int face_normal = f.first, level = f.second;
+   if (face_normal == 0) // x-normal
+   {
+      offsets = {level ? pp1 - 1 : 0};
+      strides = {pp1, p*pp1};
+   }
+   else if (face_normal == 1) // y-normal
+   {
+      offsets = {n_dof_per_dim + (level ? p*(pp1 - 1) : 0)};
+      strides = {1, p*pp1};
+   }
+   else if (face_normal == 2) // z-normal
+   {
+      offsets = {2*n_dof_per_dim + (level ? p*p*(pp1 - 1) : 0)};
+      strides = {1, p};
+   }
+   std::vector<int> n_dofs = {p, p};
+   FillFaceMap(n_face_dofs, offsets, strides, n_dofs, face_map);
+}
+
 
 const double RT_TriangleElement::nk[6] =
 { 0., -1., 1., 1., -1., 0. };

fem/fe/fe_rt.hpp

@@ -82,6 +82,8 @@ public:
                             DenseMatrix &curl) const
    { ProjectGrad_RT(nk, dof2nk, fe, Trans, curl); }
 
+   virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
+
 protected:
    void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans,
                           Vector &dofs) const;
@@ -145,6 +147,10 @@ public:
                             DenseMatrix &curl) const
    { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
 
+   /// @brief Return the mapping from lexicographically ordered DOFs to face
+   /// DOFs corresponding to local face @a face_id.
+   virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
+
 protected:
    void ProjectIntegrated(VectorCoefficient &vc,
                           ElementTransformation &Trans,

fem/fespace.cpp

@@ -1317,7 +1317,7 @@ const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
       }
       else
       {
-         res = new H1FaceRestriction(*this, e_ordering, type);
+         res = new H1_ND_RT_FaceRestriction(*this, e_ordering, type);
       }
       L2F[key] = res;
       return res;

fem/linearform.cpp

@@ -101,7 +101,7 @@ void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
    interior_face_integs.Append(lfi);
 }
 
-bool LinearForm::SupportsDevice()
+bool LinearForm::SupportsDevice() const
 {
    // return false for NURBS meshes, so we don’t convert it to non-NURBS
    // through Assemble, AssembleDevice, GetGeometricFactors and EnsureNodes

fem/linearform.hpp

@@ -203,7 +203,7 @@ public:
    void Assemble();
 
    /// Return true if assembly on device is supported, false otherwise.
-   virtual bool SupportsDevice();
+   virtual bool SupportsDevice() const;
 
    /// Assembles delta functions of the linear form
    void AssembleDelta();

fem/lininteg.hpp

@@ -31,7 +31,7 @@ protected:
 public:
 
    /// Method probing for assembly on device
-   virtual bool SupportsDevice() { return false; }
+   virtual bool SupportsDevice() const { return false; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -121,7 +121,7 @@ public:
    DomainLFIntegrator(Coefficient &QF, const IntegrationRule *ir)
       : DeltaLFIntegrator(QF, ir), Q(QF), oa(1), ob(1) { }
 
-   virtual bool SupportsDevice() { return true; }
+   virtual bool SupportsDevice() const { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -154,7 +154,7 @@ public:
    DomainLFGradIntegrator(VectorCoefficient &QF)
       : DeltaLFIntegrator(QF), Q(QF) { }
 
-   virtual bool SupportsDevice() { return true; }
+   virtual bool SupportsDevice() const { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -187,7 +187,7 @@ public:
    BoundaryLFIntegrator(Coefficient &QG, int a = 1, int b = 1)
       : Q(QG), oa(a), ob(b) { }
 
-   virtual bool SupportsDevice() { return true; }
+   virtual bool SupportsDevice() const { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -217,7 +217,7 @@ public:
    BoundaryNormalLFIntegrator(VectorCoefficient &QG, int a = 1, int b = 1)
       : Q(QG), oa(a), ob(b) { }
 
-   virtual bool SupportsDevice() { return true; }
+   virtual bool SupportsDevice() const { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -262,7 +262,7 @@ public:
    VectorDomainLFIntegrator(VectorCoefficient &QF)
       : DeltaLFIntegrator(QF), Q(QF) { }
 
-   virtual bool SupportsDevice() { return true; }
+   virtual bool SupportsDevice() const { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -296,7 +296,7 @@ public:
    VectorDomainLFGradIntegrator(VectorCoefficient &QF)
       : DeltaLFIntegrator(QF), Q(QF) { }
 
-   virtual bool SupportsDevice() override { return true; }
+   virtual bool SupportsDevice() const override { return true; }
 
    /// Method defining assembly on device
    virtual void AssembleDevice(const FiniteElementSpace &fes,
@@ -456,6 +456,12 @@ public:
                                        Vector &elvect);
 
    using LinearFormIntegrator::AssembleRHSElementVect;
+
+   virtual bool SupportsDevice() const { return true; }
+
+   virtual void AssembleDevice(const FiniteElementSpace &fes,
+                               const Array<int> &markers,
+                               Vector &b);
 };
 
 /// Class for boundary integration \f$ L(v) = (n \times f, v) \f$

fem/lininteg_boundary_flux.cpp

@@ -0,0 +1,180 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "fem.hpp"
+#include "../fem/kernels.hpp"
+#include "../general/forall.hpp"
+
+namespace mfem
+{
+
+template<int T_D1D = 0, int T_Q1D = 0> static
+void BFLFEvalAssemble2D(const int nbe, const int d, const int q,
+                        const int *markers, const double *b,
+                        const double *weights, const Vector &coeff, double *y)
+{
+   const auto F = coeff.Read();
+   const auto M = Reshape(markers, nbe);
+   const auto B = Reshape(b, q, d);
+   const auto W = Reshape(weights, q);
+   const bool const_coeff = coeff.Size() == 1;
+   const auto C = const_coeff ? Reshape(F,1,1) : Reshape(F,q,nbe);
+   auto Y = Reshape(y, d, nbe);
+
+   MFEM_FORALL(e, nbe,
+   {
+      if (M(e) == 0) { return; } // ignore
+
+      constexpr int Q = T_Q1D ? T_Q1D : MAX_Q1D;
+      double QQ[Q];
+
+      for (int qx = 0; qx < q; ++qx)
+      {
+         const double coeff_val = const_coeff ? C(0,0) : C(qx,e);
+         QQ[qx] = W(qx) * coeff_val;
+      }
+      for (int dx = 0; dx < d; ++dx)
+      {
+         double u = 0;
+         for (int qx = 0; qx < q; ++qx) { u += QQ[qx] * B(qx,dx); }
+         Y(dx,e) += u;
+      }
+   });
+}
+
+template<int T_D1D = 0, int T_Q1D = 0> static
+void BFLFEvalAssemble3D(const int nbe, const int d, const int q,
+                        const int *markers, const double *b,
+                        const double *weights, const Vector &coeff, double *y)
+{
+   const auto F = coeff.Read();
+   const auto M = Reshape(markers, nbe);
+   const auto B = Reshape(b, q, d);
+   const auto W = Reshape(weights, q, q);
+   const bool const_coeff = coeff.Size() == 1;
+   const auto C = const_coeff ? Reshape(F,1,1,1) : Reshape(F,q,q,nbe);
+   auto Y = Reshape(y, d, d, nbe);
+
+   MFEM_FORALL_2D(e, nbe, q, q, 1,
+   {
+      if (M(e) == 0) { return; } // ignore
+
+      constexpr int Q = T_Q1D ? T_Q1D : MAX_Q1D;
+      constexpr int D = T_D1D ? T_D1D : MAX_D1D;
+
+      MFEM_SHARED double sBt[Q*D];
+      MFEM_SHARED double sQQ[Q*Q];
+      MFEM_SHARED double sQD[Q*D];
+
+      const DeviceMatrix Bt(sBt, d, q);
+      kernels::internal::LoadB<D,Q>(d, q, B, sBt);
+
+      const DeviceMatrix QQ(sQQ, q, q);
+      const DeviceMatrix QD(sQD, q, d);
+
+      MFEM_FOREACH_THREAD(x,x,q)
+      {
+         MFEM_FOREACH_THREAD(y,y,q)
+         {
+            const double coeff_val = const_coeff ? C(0,0,0) : C(x,y,e);
+            QQ(y,x) = W(x,y) * coeff_val;
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(qy,y,q)
+      {
+         MFEM_FOREACH_THREAD(dx,x,d)
+         {
+            double u = 0.0;
+            for (int qx = 0; qx < q; ++qx) { u += QQ(qy,qx) * Bt(dx,qx); }
+            QD(qy,dx) = u;
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(dy,y,d)
+      {
+         MFEM_FOREACH_THREAD(dx,x,d)
+         {
+            double u = 0.0;
+            for (int qy = 0; qy < q; ++qy) { u += QD(qy,dx) * Bt(dy,qy); }
+            Y(dx,dy,e) += u;
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+static void BFLFEvalAssemble(const FiniteElementSpace &fes,
+                             const IntegrationRule &ir,
+                             const Array<int> &markers,
+                             const Vector &coeff,
+                             Vector &y)
+{
+   Mesh &mesh = *fes.GetMesh();
+   const int dim = mesh.Dimension();
+   const FiniteElement &el = *fes.GetBE(0);
+   const DofToQuad &maps = el.GetDofToQuad(ir, DofToQuad::TENSOR);
+   const int d = maps.ndof, q = maps.nqpt;
+   auto ker = (dim == 2) ? BFLFEvalAssemble2D<> : BFLFEvalAssemble3D<>;
+
+   if (dim==2)
+   {
+      if (d==1 && q==1) { ker=BFLFEvalAssemble2D<1,1>; }
+      if (d==2 && q==2) { ker=BFLFEvalAssemble2D<2,2>; }
+      if (d==3 && q==3) { ker=BFLFEvalAssemble2D<3,3>; }
+      if (d==4 && q==4) { ker=BFLFEvalAssemble2D<4,4>; }
+      if (d==5 && q==5) { ker=BFLFEvalAssemble2D<5,5>; }
+      if (d==2 && q==3) { ker=BFLFEvalAssemble2D<2,3>; }
+      if (d==3 && q==4) { ker=BFLFEvalAssemble2D<3,4>; }
+      if (d==4 && q==5) { ker=BFLFEvalAssemble2D<4,5>; }
+      if (d==5 && q==6) { ker=BFLFEvalAssemble2D<5,6>; }
+   }
+
+   if (dim==3)
+   {
+      if (d==1 && q==1) { ker=BFLFEvalAssemble3D<1,1>; }
+      if (d==2 && q==2) { ker=BFLFEvalAssemble3D<2,2>; }
+      if (d==3 && q==3) { ker=BFLFEvalAssemble3D<3,3>; }
+      if (d==4 && q==4) { ker=BFLFEvalAssemble3D<4,4>; }
+      if (d==5 && q==5) { ker=BFLFEvalAssemble3D<5,5>; }
+      if (d==2 && q==3) { ker=BFLFEvalAssemble3D<2,3>; }
+      if (d==3 && q==4) { ker=BFLFEvalAssemble3D<3,4>; }
+      if (d==4 && q==5) { ker=BFLFEvalAssemble3D<4,5>; }
+      if (d==5 && q==6) { ker=BFLFEvalAssemble3D<5,6>; }
+   }
+
+   MFEM_VERIFY(ker, "No kernel ndof " << d << " nqpt " << q);
+
+   const int nbe = fes.GetMesh()->GetNFbyType(FaceType::Boundary);
+   const int *M = markers.Read();
+   const double *B = maps.B.Read();
+   const double *W = ir.GetWeights().Read();
+   double *Y = y.ReadWrite();
+   ker(nbe, d, q, M, B, W, coeff, Y);
+}
+
+void VectorFEBoundaryFluxLFIntegrator::AssembleDevice(
+   const FiniteElementSpace &fes,
+   const Array<int> &markers,
+   Vector &b)
+{
+   const FiniteElement &fe = *fes.GetBE(0);
+   const int qorder = oa * fe.GetOrder() + ob;
+   const Geometry::Type gtype = fe.GetGeomType();
+   const IntegrationRule &ir = IntRule ? *IntRule : IntRules.Get(gtype, qorder);
+   Mesh &mesh = *fes.GetMesh();
+
+   FaceQuadratureSpace qs(mesh, ir, FaceType::Boundary);
+   CoefficientVector coeff(F, qs, CoefficientStorage::COMPRESSED);
+   BFLFEvalAssemble(fes, ir, markers, coeff, b);
+}
+
+} // namespace mfem

fem/pfespace.cpp

@@ -569,7 +569,7 @@ const FaceRestriction *ParFiniteElementSpace::GetFaceRestriction(
       {
          if (Conforming())
          {
-            res = new H1FaceRestriction(*this, e_ordering, type);
+            res = new H1_ND_RT_FaceRestriction(*this, e_ordering, type);
          }
          else
          {

fem/plinearform.cpp

@@ -54,7 +54,7 @@ void ParLinearForm::Assemble()
    }
 }
 
-bool ParLinearForm::SupportsDevice()
+bool ParLinearForm::SupportsDevice() const
 {
    bool parallel;
    bool local = LinearForm::SupportsDevice();

fem/plinearform.hpp

@@ -120,7 +120,7 @@ public:
    void Assemble();
 
    /// Return true if assembly on device is supported, false otherwise.
-   virtual bool SupportsDevice();
+   virtual bool SupportsDevice() const;
 
    void AssembleSharedFaces();
 

fem/restriction.cpp

@@ -594,119 +594,11 @@ void L2ElementRestriction::FillJAndData(const Vector &ea_data,
    });
 }
 
-/** Return the face degrees of freedom returned in Lexicographic order.
-    Note: Only for quad and hex */
-void GetFaceDofs(const int dim, const int face_id,
-                 const int dof1d, Array<int> &face_map)
-{
-   switch (dim)
-   {
-      case 1:
-         switch (face_id)
-         {
-            case 0: // WEST
-               face_map[0] = 0;
-               break;
-            case 1: // EAST
-               face_map[0] = dof1d-1;
-               break;
-         }
-         break;
-      case 2:
-         switch (face_id)
-         {
-            case 0: // SOUTH
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  face_map[i] = i;
-               }
-               break;
-            case 1: // EAST
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  face_map[i] = dof1d-1 + i*dof1d;
-               }
-               break;
-            case 2: // NORTH
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  face_map[i] = (dof1d-1)*dof1d + i;
-               }
-               break;
-            case 3: // WEST
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  face_map[i] = i*dof1d;
-               }
-               break;
-         }
-         break;
-      case 3:
-         switch (face_id)
-         {
-            case 0: // BOTTOM
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = i + j*dof1d;
-                  }
-               }
-               break;
-            case 1: // SOUTH
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = i + j*dof1d*dof1d;
-                  }
-               }
-               break;
-            case 2: // EAST
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = dof1d-1 + i*dof1d + j*dof1d*dof1d;
-                  }
-               }
-               break;
-            case 3: // NORTH
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = (dof1d-1)*dof1d + i + j*dof1d*dof1d;
-                  }
-               }
-               break;
-            case 4: // WEST
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = i*dof1d + j*dof1d*dof1d;
-                  }
-               }
-               break;
-            case 5: // TOP
-               for (int i = 0; i < dof1d; ++i)
-               {
-                  for (int j = 0; j < dof1d; ++j)
-                  {
-                     face_map[i+j*dof1d] = (dof1d-1)*dof1d*dof1d + i + j*dof1d;
-                  }
-               }
-               break;
-         }
-         break;
-   }
-}
-
-H1FaceRestriction::H1FaceRestriction(const FiniteElementSpace &fes,
-                                     const ElementDofOrdering e_ordering,
-                                     const FaceType type,
-                                     bool build)
+H1_ND_RT_FaceRestriction::H1_ND_RT_FaceRestriction(
+   const FiniteElementSpace &fes,
+   const ElementDofOrdering e_ordering,
+   const FaceType type,
+   bool build)
    : fes(fes),
      nf(fes.GetNFbyType(type)),
      vdim(fes.GetVDim()),
@@ -727,18 +619,36 @@ H1FaceRestriction::H1FaceRestriction(const FiniteElementSpace &fes,
 
    CheckFESpace(e_ordering);
 
-   ComputeScatterIndicesAndOffsets(e_ordering, type);
+   // Get the mapping from native DOF ordering to lexicographic ordering.
+   const FiniteElement *fe = fes.GetFE(0);
+   const TensorBasisElement* el =
+      dynamic_cast<const TensorBasisElement*>(fe);
+   const Array<int> &dof_map_ = el->GetDofMap();
+   if (dof_map_.Size() > 0)
+   {
+      dof_map.MakeRef(dof_map_);
+   }
+   else
+   {
+      // For certain types of elements dof_map_ is empty, in this case that
+      // means the element is already ordered lexicographically, so the
+      // permutation is the identity.
+      dof_map.SetSize(elem_dofs);
+      for (int i = 0; i < elem_dofs; ++i) { dof_map[i] = i; }
+   }
 
+   ComputeScatterIndicesAndOffsets(e_ordering, type);
    ComputeGatherIndices(e_ordering,type);
 }
 
-H1FaceRestriction::H1FaceRestriction(const FiniteElementSpace &fes,
-                                     const ElementDofOrdering e_ordering,
-                                     const FaceType type)
-   : H1FaceRestriction(fes, e_ordering, type, true)
+H1_ND_RT_FaceRestriction::H1_ND_RT_FaceRestriction(
+   const FiniteElementSpace &fes,
+   const ElementDofOrdering e_ordering,
+   const FaceType type)
+   : H1_ND_RT_FaceRestriction(fes, e_ordering, type, true)
 { }
 
-void H1FaceRestriction::Mult(const Vector& x, Vector& y) const
+void H1_ND_RT_FaceRestriction::Mult(const Vector& x, Vector& y) const
 {
    if (nf==0) { return; }
    // Assumes all elements have the same number of dofs
@@ -750,17 +660,20 @@ void H1FaceRestriction::Mult(const Vector& x, Vector& y) const
    auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
    MFEM_FORALL(i, nfdofs,
    {
-      const int idx = d_indices[i];
+      const int s_idx = d_indices[i];
+      const int sgn = (s_idx >= 0) ? 1 : -1;
+      const int idx = (s_idx >= 0) ? s_idx : -1 - s_idx;
       const int dof = i % nface_dofs;
       const int face = i / nface_dofs;
       for (int c = 0; c < vd; ++c)
       {
-         d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
+         d_y(dof, c, face) = sgn*d_x(t?c:idx, t?idx:c);
       }
    });
 }
 
-void H1FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
+void H1_ND_RT_FaceRestriction::AddMultTranspose(
+   const Vector& x, Vector& y) const
 {
    if (nf==0) { return; }
    // Assumes all elements have the same number of dofs
@@ -780,15 +693,17 @@ void H1FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
          double dof_value = 0;
          for (int j = offset; j < next_offset; ++j)
          {
-            const int idx_j = d_indices[j];
-            dof_value +=  d_x(idx_j % nface_dofs, c, idx_j / nface_dofs);
+            const int s_idx_j = d_indices[j];
+            const int sgn = (s_idx_j >= 0) ? 1 : -1;
+            const int idx_j = (s_idx_j >= 0) ? s_idx_j : -1 - s_idx_j;
+            dof_value += sgn*d_x(idx_j % nface_dofs, c, idx_j / nface_dofs);
          }
          d_y(t?c:i,t?i:c) += dof_value;
       }
    });
 }
 
-void H1FaceRestriction::CheckFESpace(const ElementDofOrdering e_ordering)
+void H1_ND_RT_FaceRestriction::CheckFESpace(const ElementDofOrdering e_ordering)
 {
 #ifdef MFEM_USE_MPI
 
@@ -810,7 +725,7 @@ void H1FaceRestriction::CheckFESpace(const ElementDofOrdering e_ordering)
                (tfe->GetBasisType()==BasisType::GaussLobatto ||
                 tfe->GetBasisType()==BasisType::Positive),
                "Only Gauss-Lobatto and Bernstein basis are supported in "
-               "H1FaceRestriction.");
+               "H1_ND_RT_FaceRestriction.");
 
    // Assuming all finite elements are using Gauss-Lobatto.
    const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
@@ -824,16 +739,11 @@ void H1FaceRestriction::CheckFESpace(const ElementDofOrdering e_ordering)
          if (el) { continue; }
          MFEM_ABORT("Finite element not suitable for lexicographic ordering");
       }
-      const FiniteElement *fe = fes.GetFaceElement(0);
-      const TensorBasisElement* el =
-         dynamic_cast<const TensorBasisElement*>(fe);
-      const Array<int> &fe_dof_map = el->GetDofMap();
-      MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
    }
 #endif
 }
 
-void H1FaceRestriction::ComputeScatterIndicesAndOffsets(
+void H1_ND_RT_FaceRestriction::ComputeScatterIndicesAndOffsets(
    const ElementDofOrdering ordering,
    const FaceType type)
 {
@@ -871,7 +781,7 @@ void H1FaceRestriction::ComputeScatterIndicesAndOffsets(
    }
 }
 
-void H1FaceRestriction::ComputeGatherIndices(
+void H1_ND_RT_FaceRestriction::ComputeGatherIndices(
    const ElementDofOrdering ordering,
    const FaceType type)
 {
@@ -904,7 +814,9 @@ void H1FaceRestriction::ComputeGatherIndices(
    gather_offsets[0] = 0;
 }
 
-void H1FaceRestriction::SetFaceDofsScatterIndices(
+static inline int absdof(int i) { return i < 0 ? -1-i : i; }
+
+void H1_ND_RT_FaceRestriction::SetFaceDofsScatterIndices(
    const Mesh::FaceInformation &face,
    const int face_index,
    const ElementDofOrdering ordering)
@@ -914,32 +826,29 @@ void H1FaceRestriction::SetFaceDofsScatterIndices(
    MFEM_ASSERT(face.element[0].orientation==0,
                "FaceRestriction used on degenerated mesh.");
 
-   const TensorBasisElement* el =
-      dynamic_cast<const TensorBasisElement*>(fes.GetFE(0));
-   const int *dof_map = el->GetDofMap().GetData();
+   fes.GetFE(0)->GetFaceMap(face.element[0].local_face_id, face_map);
+
    const Table& e2dTable = fes.GetElementToDofTable();
    const int* elem_map = e2dTable.GetJ();
-   const int face_id = face.element[0].local_face_id;
-   const int dim = fes.GetMesh()->Dimension();
-   const int dof1d = fes.GetFE(0)->GetOrder()+1;
    const int elem_index = face.element[0].index;
    const bool dof_reorder = (ordering == ElementDofOrdering::LEXICOGRAPHIC);
-   GetFaceDofs(dim, face_id, dof1d, face_map); // Only for quad and hex
 
    for (int face_dof = 0; face_dof < face_dofs; ++face_dof)
    {
       const int nat_volume_dof = face_map[face_dof];
-      const int volume_dof = (!dof_reorder)?
-                             nat_volume_dof:
-                             dof_map[nat_volume_dof];
-      const int global_dof = elem_map[elem_index*elem_dofs + volume_dof];
+      const int s_volume_dof = (!dof_reorder) ?
+                               nat_volume_dof :
+                               dof_map[nat_volume_dof]; // signed
+      const int volume_dof = absdof(s_volume_dof);
+      const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
+      const int global_dof = absdof(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
-      scatter_indices[restriction_dof] = global_dof;
+      scatter_indices[restriction_dof] = s_global_dof;
       ++gather_offsets[global_dof + 1];
    }
 }
 
-void H1FaceRestriction::SetFaceDofsGatherIndices(
+void H1_ND_RT_FaceRestriction::SetFaceDofsGatherIndices(
    const Mesh::FaceInformation &face,
    const int face_index,
    const ElementDofOrdering ordering)
@@ -947,25 +856,25 @@ void H1FaceRestriction::SetFaceDofsGatherIndices(
    MFEM_ASSERT(!(face.IsNonconformingCoarse()),
                "This method should not be used on nonconforming coarse faces.");
 
-   const TensorBasisElement* el =
-      dynamic_cast<const TensorBasisElement*>(fes.GetFE(0));
-   const int *dof_map = el->GetDofMap().GetData();
+   fes.GetFE(0)->GetFaceMap(face.element[0].local_face_id, face_map);
+
    const Table& e2dTable = fes.GetElementToDofTable();
    const int* elem_map = e2dTable.GetJ();
-   const int face_id = face.element[0].local_face_id;
-   const int dim = fes.GetMesh()->Dimension();
-   const int dof1d = fes.GetFE(0)->GetOrder()+1;
    const int elem_index = face.element[0].index;
    const bool dof_reorder = (ordering == ElementDofOrdering::LEXICOGRAPHIC);
-   GetFaceDofs(dim, face_id, dof1d, face_map); // Only for quad and hex
 
    for (int face_dof = 0; face_dof < face_dofs; ++face_dof)
    {
       const int nat_volume_dof = face_map[face_dof];
-      const int volume_dof = (!dof_reorder)?nat_volume_dof:dof_map[nat_volume_dof];
-      const int global_dof = elem_map[elem_index*elem_dofs + volume_dof];
+      const int s_volume_dof = (!dof_reorder)?nat_volume_dof:dof_map[nat_volume_dof];
+      const int volume_dof = absdof(s_volume_dof);
+      const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
+      const int sgn = (s_global_dof >= 0) ? 1 : -1;
+      const int global_dof = absdof(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
-      gather_indices[gather_offsets[global_dof]++] = restriction_dof;
+      const int s_restriction_dof = (sgn >= 0) ? restriction_dof : -1 -
+                                    restriction_dof;
+      gather_indices[gather_offsets[global_dof]++] = s_restriction_dof;
    }
 }
 
@@ -1517,10 +1426,8 @@ void L2FaceRestriction::SetFaceDofsScatterIndices1(
    const Table& e2dTable = fes.GetElementToDofTable();
    const int* elem_map = e2dTable.GetJ();
    const int face_id1 = face.element[0].local_face_id;
-   const int dim = fes.GetMesh()->Dimension();
-   const int dof1d = fes.GetFE(0)->GetOrder()+1;
    const int elem_index = face.element[0].index;
-   GetFaceDofs(dim, face_id1, dof1d, face_map); // Only for quad and hex
+   fes.GetFE(0)->GetFaceMap(face_id1, face_map);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
@@ -1546,7 +1453,7 @@ void L2FaceRestriction::PermuteAndSetFaceDofsScatterIndices2(
    const int orientation = face.element[1].orientation;
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetFE(0)->GetOrder()+1;
-   GetFaceDofs(dim, face_id2, dof1d, face_map); // Only for quad and hex
+   fes.GetFE(0)->GetFaceMap(face_id2, face_map);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
@@ -1574,7 +1481,7 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
    const int orientation = face.element[1].orientation;
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetFE(0)->GetOrder()+1;
-   GetFaceDofs(dim, face_id2, dof1d, face_map); // Only for quad and hex
+   fes.GetFE(0)->GetFaceMap(face_id2, face_map);
    Array<int> face_nbr_dofs;
    const ParFiniteElementSpace &pfes =
       static_cast<const ParFiniteElementSpace&>(this->fes);
@@ -1616,10 +1523,8 @@ void L2FaceRestriction::SetFaceDofsGatherIndices1(
    const Table& e2dTable = fes.GetElementToDofTable();
    const int* elem_map = e2dTable.GetJ();
    const int face_id1 = face.element[0].local_face_id;
-   const int dim = fes.GetMesh()->Dimension();
-   const int dof1d = fes.GetFE(0)->GetOrder()+1;
    const int elem_index = face.element[0].index;
-   GetFaceDofs(dim, face_id1, dof1d, face_map); // Only for quad and hex
+   fes.GetFE(0)->GetFaceMap(face_id1, face_map);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
@@ -1645,7 +1550,7 @@ void L2FaceRestriction::PermuteAndSetFaceDofsGatherIndices2(
    const int orientation = face.element[1].orientation;
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetFE(0)->GetOrder()+1;
-   GetFaceDofs(dim, face_id2, dof1d, face_map); // Only for quad and hex
+   fes.GetFE(0)->GetFaceMap(face_id2, face_map);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {

fem/restriction.hpp

@@ -216,7 +216,7 @@ public:
 /// Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
 /** Objects of this type are typically created and owned by FiniteElementSpace
     objects, see FiniteElementSpace::GetFaceRestriction(). */
-class H1FaceRestriction : public FaceRestriction
+class H1_ND_RT_FaceRestriction : public FaceRestriction
 {
 protected:
    const FiniteElementSpace &fes;
@@ -230,29 +230,30 @@ protected:
    Array<int> scatter_indices; // Scattering indices for element 1 on each face
    Array<int> gather_offsets; // offsets for the gathering indices of each dof
    Array<int> gather_indices; // gathering indices for each dof
+   Array<int> dof_map; // mapping to lexicographic ordering
 
-   /** @brief Construct an H1FaceRestriction.
+   /** @brief Construct an H1_ND_RT_FaceRestriction.
 
        @param[in] fes      The FiniteElementSpace on which this operates
        @param[in] ordering Request a specific element ordering
        @param[in] type     Request internal or boundary faces dofs
        @param[in] build    Request the NCL2FaceRestriction to compute the
                            scatter/gather indices. False should only be used
-                           when inheriting from H1FaceRestriction.
+                           when inheriting from H1_ND_RT_FaceRestriction.
    */
-   H1FaceRestriction(const FiniteElementSpace& fes,
-                     const ElementDofOrdering ordering,
-                     const FaceType type,
-                     bool build);
+   H1_ND_RT_FaceRestriction(const FiniteElementSpace& fes,
+                            const ElementDofOrdering ordering,
+                            const FaceType type,
+                            bool build);
 public:
-   /** @brief Construct an H1FaceRestriction.
+   /** @brief Construct an H1_ND_RT_FaceRestriction.
 
        @param[in] fes      The FiniteElementSpace on which this operates
        @param[in] ordering Request a specific element ordering
        @param[in] type     Request internal or boundary faces dofs */
-   H1FaceRestriction(const FiniteElementSpace& fes,
-                     const ElementDofOrdering ordering,
-                     const FaceType type);
+   H1_ND_RT_FaceRestriction(const FiniteElementSpace& fes,
+                            const ElementDofOrdering ordering,
+                            const FaceType type);
 
    /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
        face E-Vector.
@@ -303,7 +304,7 @@ private:
 protected:
    mutable Array<int> face_map; // Used in the computation of GetFaceDofs
 
-   /** @brief Verify that H1FaceRestriction is build from an H1 FESpace.
+   /** @brief Verify that H1_ND_RT_FaceRestriction is build from an H1 FESpace.
 
        @param[in] ordering The FESpace element ordering.
    */
@@ -332,6 +333,10 @@ protected:
                                  const ElementDofOrdering ordering);
 };
 
+/// @brief Alias for H1_ND_RT_FaceRestriction, for backwards compatibility and
+/// as base class for ParNCH1FaceRestriction.
+using H1FaceRestriction = H1_ND_RT_FaceRestriction;
+
 /// Operator that extracts Face degrees of freedom for L2 spaces.
 /** Objects of this type are typically created and owned by FiniteElementSpace
     objects, see FiniteElementSpace::GetFaceRestriction(). */

mesh/mesh.cpp

@@ -1319,6 +1319,7 @@ std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
          os << "NA";
          break;
    }
+   os << '\n';
    os << "element[0].location=";
    switch (info.element[0].location)
    {
@@ -1332,7 +1333,7 @@ std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
          os << "NA";
          break;
    }
-   os << std::endl;
+   os << '\n';
    os << "element[1].location=";
    switch (info.element[1].location)
    {
@@ -1346,7 +1347,7 @@ std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
          os << "NA";
          break;
    }
-   os << std::endl;
+   os << '\n';
    os << "element[0].conformity=";
    switch (info.element[0].conformity)
    {
@@ -1363,7 +1364,7 @@ std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
          os << "NA";
          break;
    }
-   os << std::endl;
+   os << '\n';
    os << "element[1].conformity=";
    switch (info.element[1].conformity)
    {
@@ -1380,13 +1381,13 @@ std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
          os << "NA";
          break;
    }
-   os << std::endl;
-   os << "element[0].index=" << info.element[0].index << std::endl
-      << "element[1].index=" << info.element[1].index << std::endl
-      << "element[0].local_face_id=" << info.element[0].local_face_id << std::endl
-      << "element[1].local_face_id=" << info.element[1].local_face_id << std::endl
-      << "element[0].orientation=" << info.element[0].orientation << std::endl
-      << "element[1].orientation=" << info.element[1].orientation << std::endl
+   os << '\n';
+   os << "element[0].index=" << info.element[0].index << '\n'
+      << "element[1].index=" << info.element[1].index << '\n'
+      << "element[0].local_face_id=" << info.element[0].local_face_id << '\n'
+      << "element[1].local_face_id=" << info.element[1].local_face_id << '\n'
+      << "element[0].orientation=" << info.element[0].orientation << '\n'
+      << "element[1].orientation=" << info.element[1].orientation << '\n'
       << "ncface=" << info.ncface << std::endl;
    return os;
 }

tests/unit/CMakeLists.txt

@@ -68,6 +68,7 @@ set(UNIT_TESTS_SRCS
   fem/test_estimator.cpp
   fem/test_face_elem_trans.cpp
   fem/test_face_permutation.cpp
+  fem/test_face_restriction.cpp
   fem/test_fe.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp

tests/unit/fem/test_face_restriction.cpp

@@ -0,0 +1,113 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+Mesh MakeCartesianMesh(int nx, int dim)
+{
+   if (dim == 2)
+   {
+      return Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL, true);
+   }
+   else
+   {
+      return Mesh::MakeCartesian3D(nx, nx, nx, Element::HEXAHEDRON);
+   }
+}
+
+TEST_CASE("Vector FE Face Restriction", "[FaceRestriction]")
+{
+   enum class SpaceType {RT, ND};
+   const auto space_type = GENERATE(SpaceType::RT, SpaceType::ND);
+   const int dim = GENERATE(2, 3);
+   const int nx = 3;
+   const int order = 4;
+
+   CAPTURE(dim);
+
+   Mesh mesh = MakeCartesianMesh(nx, dim);
+
+   int ndof_per_face;
+   std::unique_ptr<FiniteElementCollection> fec;
+   if (space_type == SpaceType::RT)
+   {
+      fec.reset(new RT_FECollection(order-1, dim));
+      ndof_per_face = int(pow(order, dim-1));
+   }
+   else
+   {
+      fec.reset(new ND_FECollection(order, dim));
+      ndof_per_face = (dim - 1)*order*int(pow(order + 1, dim - 2));
+   }
+
+   FiniteElementSpace fes(&mesh, fec.get());
+
+   auto ordering = ElementDofOrdering::LEXICOGRAPHIC;
+   auto ftype = FaceType::Boundary;
+   const int nfaces = fes.GetNFbyType(FaceType::Boundary);
+   const FaceRestriction *face_restr =
+      fes.GetFaceRestriction(ordering, ftype);
+
+   REQUIRE(face_restr != nullptr);
+
+   Array<int> bdr_dofs;
+   fes.GetBoundaryTrueDofs(bdr_dofs);
+
+   // Set gf to have random values on the boundary, zero on the interior
+   GridFunction gf(&fes);
+   gf.Randomize(0);
+   gf.SetSubVectorComplement(bdr_dofs, 0.0);
+
+   // Mapping to face E-vector and back to L-vector should give back the
+   // original grid function.
+   Vector face_vec(face_restr->Height());
+   REQUIRE(face_vec.Size() == nfaces*ndof_per_face);
+   face_restr->Mult(gf, face_vec);
+
+   if (space_type == SpaceType::ND && dim == 3)
+   {
+      // Adjust for multiplicity. In all other cases, each boundary DOF is
+      // unique (not shared between faces). In the case of 3D ND elements, some
+      // boundary DOFs are shared between two faces (i.e. those that lie on
+      // element edges).
+      //
+      // This adjustment will ensure that the original vector is recovered after
+      // multiplying by the transpose of the face restriction operator.
+      const int n = order*(order+1);
+      for (int f = 0; f < fes.GetNFbyType(ftype); ++f)
+      {
+         for (int d = 0; d < 2; ++d)
+         {
+            const int nx = (d == 0) ? order : order + 1;
+            const int ny = (d == 0) ? order + 1 : order;
+            for (int i = 0; i < n; ++i)
+            {
+               const int ix = i % nx;
+               const int iy = i / nx;
+               if ((d == 0 && (iy == 0 || iy == ny - 1)) ||
+                   (d == 1 && (ix == 0 || ix == nx - 1)))
+               {
+                  face_vec[f*ndof_per_face + d*n + i] *= 0.5;
+               }
+            }
+         }
+      }
+   }
+
+   GridFunction gf2(&fes);
+   face_restr->MultTranspose(face_vec, gf2);
+
+   gf2 -= gf;
+   REQUIRE(gf2.Normlinf() == MFEM_Approx(0.0));
+}

tests/unit/fem/test_linearform_ext.cpp

@@ -282,4 +282,33 @@ TEST_CASE("Linear Form Extension", "[LinearFormExtension], [CUDA]")
 
       REQUIRE(d1.Norml2() == MFEM_Approx(0.0));
    }
+
+   SECTION("VectorFE")
+   {
+      Mesh mesh(mesh_file);
+      const int dim = mesh.Dimension();
+
+      CAPTURE(mesh_file, dim, p);
+
+      RT_FECollection fec(p-1, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      FunctionCoefficient coeff(f);
+
+      LinearForm d1(&fes);
+      d1.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(coeff));
+      d1.UseFastAssembly(true);
+      d1.Assemble();
+
+      LinearForm d2(&fes);
+      d2.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(coeff));
+      d2.UseFastAssembly(false);
+      d2.Assemble();
+
+      CAPTURE(d1.Norml2(), d2.Norml2());
+
+      d1 -= d2;
+
+      REQUIRE(d1.Norml2() == MFEM_Approx(0.0));
+   }
 }