diff --git a/fem/fe.cpp b/fem/fe.cpp
index d9ac775ae3..7dcc8dfed5 100644
--- a/fem/fe.cpp
+++ b/fem/fe.cpp
@@ -10984,11 +10984,13 @@ RT_HexahedronElement::RT_HexahedronElement(const int p,
                                            const int ob_type)
    : VectorTensorFiniteElement(3, 3*(p + 1)*(p + 1)*(p + 2), p + 1, cb_type,
                                ob_type, H_DIV, DofMapType::L2_DOF_MAP),
-     dof2nk(dof)
+     dof2nk(dof),
+     cp(poly1d.ClosedPoints(p + 1, cb_type))
 {
+   if (obasis1d.IsIntegratedType()) { is_nodal = false; }
+
    dof_map.SetSize(dof);
 
-   const double *cp = poly1d.ClosedPoints(p + 1, cb_type);
    const double *op = poly1d.OpenPoints(p, ob_type);
    const int dof3 = dof/3;
 
@@ -11307,6 +11309,56 @@ void RT_HexahedronElement::CalcDivShape(const IntegrationPoint &ip,
          }
 }
 
+void RT_HexahedronElement::ProjectIntegrated(VectorCoefficient &vc,
+                                             ElementTransformation &Trans,
+                                             Vector &dofs) const
+{
+   MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
+   double vq[Geometry::MaxDim];
+   Vector xq(vq, vc.GetVDim());
+
+   const IntegrationRule &ir2d = IntRules.Get(Geometry::SQUARE, order);
+   const int nqpt = ir2d.GetNPoints();
+
+   IntegrationPoint ip3d;
+
+   int o = 0;
+   for (int c = 0; c < 3; c++)
+   {
+      int im = (c == 0) ? order + 1 : order;
+      int jm = (c == 1) ? order + 1 : order;
+      int km = (c == 2) ? order + 1 : order;
+      for (int k = 0; k < km; k++)
+         for (int j = 0; j < jm; j++)
+            for (int i = 0; i < im; i++)
+            {
+               int idx = dof_map[o++];
+               if (idx < 0) { idx = -1 - idx; }
+               int ic1, ic2;
+               if (c == 0) { ic1 = j; ic2 = k; }
+               else if (c == 1) { ic1 = i; ic2 = k; }
+               else { ic1 = i; ic2 = j; }
+               const double h1 = cp[ic1+1] - cp[ic1];
+               const double h2 = cp[ic2+1] - cp[ic2];
+               double val = 0.0;
+               for (int q = 0; q < nqpt; q++)
+               {
+                  const IntegrationPoint &ip2d = ir2d.IntPoint(q);
+                  if (c == 0) { ip3d.Set3(cp[i], cp[j] + h1*ip2d.x, cp[j] + h2*ip2d.y); }
+                  else if (c == 1) { ip3d.Set3(cp[i] + h1*ip2d.x, cp[j], cp[k] + h2*ip2d.y); }
+                  else { ip3d.Set3(cp[i] + h1*ip2d.x, cp[j] + h2*ip2d.y, cp[k]); }
+                  Trans.SetIntPoint(&ip3d);
+                  vc.Eval(xq, Trans, ip3d);
+                  // nk^t adj(J) xq
+                  const double ipval
+                     = Trans.AdjugateJacobian().InnerProduct(vq, nk + dof2nk[idx]*dim);
+                  val += ip2d.weight*ipval;
+               }
+               dofs(idx) = val*h1*h2;
+            }
+   }
+}
+
 
 const double RT_TriangleElement::nk[6] =
 { 0., -1., 1., 1., -1., 0. };
diff --git a/fem/fe.hpp b/fem/fe.hpp
index 2a31a1090a..88c5ce5b4f 100644
--- a/fem/fe.hpp
+++ b/fem/fe.hpp
@@ -2858,6 +2858,7 @@ class RT_HexahedronElement : public VectorTensorFiniteElement
    mutable Vector dshape_cx, dshape_cy, dshape_cz;
 #endif
    Array<int> dof2nk;
+   const double *cp;
 
 public:
    /** @brief Construct the RT_HexahedronElement of order @a p and closed and
@@ -2886,7 +2887,10 @@ public:
    using FiniteElement::Project;
    virtual void Project(VectorCoefficient &vc,
                         ElementTransformation &Trans, Vector &dofs) const
-   { Project_RT(nk, dof2nk, vc, Trans, dofs); }
+   {
+      if (obasis1d.IsIntegratedType()) { ProjectIntegrated(vc, Trans, dofs); }
+      else { Project_RT(nk, dof2nk, vc, Trans, dofs); }
+   }
    virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
                                  Vector &dofs) const
    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
@@ -2900,6 +2904,11 @@ public:
                             ElementTransformation &Trans,
                             DenseMatrix &curl) const
    { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+
+protected:
+   void ProjectIntegrated(VectorCoefficient &vc,
+                          ElementTransformation &Trans,
+                          Vector &dofs) const;
 };