diff --git a/fem/fe/fe_pyramid.cpp b/fem/fe/fe_pyramid.cpp
index 7e4c37f64b..9be70436c6 100644
--- a/fem/fe/fe_pyramid.cpp
+++ b/fem/fe/fe_pyramid.cpp
@@ -83,6 +83,38 @@ DenseMatrix FuentesPyramid::grad_lam45(real_t x, real_t y, real_t z)
    return dlam;
 }
 
+Vector FuentesPyramid::lam15_grad_lam15(real_t x, real_t y, real_t z)
+{
+   Vector lam = lam15(x, y, z);
+   Vector lamdlam(3);
+   add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam1(x, y, z), lamdlam);
+   return lamdlam;
+}
+
+Vector FuentesPyramid::lam25_grad_lam25(real_t x, real_t y, real_t z)
+{
+   Vector lam = lam25(x, y, z);
+   Vector lamdlam(3);
+   add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam2(x, y, z), lamdlam);
+   return lamdlam;
+}
+
+Vector FuentesPyramid::lam35_grad_lam35(real_t x, real_t y, real_t z)
+{
+   Vector lam = lam35(x, y, z);
+   Vector lamdlam(3);
+   add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam3(x, y, z), lamdlam);
+   return lamdlam;
+}
+
+Vector FuentesPyramid::lam45_grad_lam45(real_t x, real_t y, real_t z)
+{
+   Vector lam = lam45(x, y, z);
+   Vector lamdlam(3);
+   add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam4(x, y, z), lamdlam);
+   return lamdlam;
+}
+
 Vector FuentesPyramid::lam125_grad_lam125(real_t x, real_t y, real_t z)
 {
    Vector lgl({-x * z / (1.0 - z), y - 1.0, z});
@@ -223,6 +255,14 @@ Vector FuentesPyramid::nu01_grad_nu01(real_t z, Vector xy, unsigned int ab)
    return nudnu;
 }
 
+Vector FuentesPyramid::nu12_grad_nu12(real_t z, Vector xy, unsigned int ab)
+{
+   Vector nu = nu12(z, xy, ab);
+   Vector nudnu(3);
+   add(nu(0), grad_nu2(z, xy, ab), -nu(1), grad_nu1(z, xy, ab), nudnu);
+   return nudnu;
+}
+
 Vector FuentesPyramid::nu012_grad_nu012(real_t z, Vector xy, unsigned int ab)
 {
    Vector nu(nu012(z, xy, ab));
@@ -1039,6 +1079,40 @@ void FuentesPyramid::E_Q(int p, Vector s, const DenseMatrix &grad_s,
       }
 }
 
+void FuentesPyramid::E_T(int p, Vector s, Vector sds, DenseTensor &u) const
+{
+   MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
+   MFEM_ASSERT(s.Size() >= 3, "Size of s must be 3 or larger");
+   MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
+   MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
+   MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
+   MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
+
+#ifdef MFEM_THREAD_SAFE
+   Vector      E_T_vtmp;
+   DenseMatrix E_T_mtmp;
+#endif
+   Vector &L_j = E_T_vtmp;
+   DenseMatrix &E_E_i = E_T_mtmp;
+
+   E_E_i.SetSize(p, 3);
+   E_E(p, s, sds, E_E_i);
+
+   L_j.SetSize(p);
+   for (int i=0; i<p; i++)
+   {
+      const real_t alpha = 2.0 * i + 1.0;
+      CalcHomogenizedIntJacobi(p - 1, alpha, s[0] + s[1], s[2], L_j);
+
+      u(i, 0, 0) = 0.0; u(i, 0, 1) = 0.0; u(i, 0, 2) = 0.0;
+      for (int j=1; i+j<p; j++)
+         for (int k=0; k<3; k++)
+         {
+            u(i, j, k) = L_j(j) * E_E_i(i, k);
+         }
+   }
+}
+
 void FuentesPyramid::V_Q(int p, Vector s, Vector sds,
                          Vector t, Vector tdt, DenseTensor &u) const
 {
diff --git a/fem/fe/fe_pyramid.hpp b/fem/fe/fe_pyramid.hpp
index 9bddc1334c..5fc80975f3 100644
--- a/fem/fe/fe_pyramid.hpp
+++ b/fem/fe/fe_pyramid.hpp
@@ -36,6 +36,8 @@ private:
    mutable DenseMatrix E_Q_mtmp1;
    mutable DenseMatrix E_Q_mtmp2;
    mutable DenseMatrix E_Q_mtmp3;
+   mutable Vector      E_T_vtmp;
+   mutable DenseMatrix E_T_mtmp;
    mutable DenseMatrix V_Q_mtmp1;
    mutable DenseMatrix V_Q_mtmp2;
    mutable Vector      V_T_vtmp1;
@@ -91,6 +93,11 @@ public:
    static DenseMatrix grad_lam35(real_t x, real_t y, real_t z);
    static DenseMatrix grad_lam45(real_t x, real_t y, real_t z);
 
+   static Vector lam15_grad_lam15(real_t x, real_t y, real_t z);
+   static Vector lam25_grad_lam25(real_t x, real_t y, real_t z);
+   static Vector lam35_grad_lam35(real_t x, real_t y, real_t z);
+   static Vector lam45_grad_lam45(real_t x, real_t y, real_t z);
+
    /// Three component vectors associated with triangular faces
    static Vector lam125(real_t x, real_t y, real_t z)
    { return Vector({lam1(x, y, z), lam2(x, y, z), lam5(x, y, z)}); }
@@ -161,13 +168,18 @@ public:
 
    static Vector nu01(real_t z, Vector xy, unsigned int ab)
    { return Vector({nu0(z, xy, ab), nu1(z, xy, ab)}); }
+   static Vector nu12(real_t z, Vector xy, unsigned int ab)
+   { return Vector({nu1(z, xy, ab), nu2(z, xy, ab)}); }
    static Vector nu012(real_t z, Vector xy, unsigned int ab)
    { return Vector({nu0(z, xy, ab), nu1(z, xy, ab), nu2(z, xy, ab)}); }
+   static Vector nu120(real_t z, Vector xy, unsigned int ab)
+   { return Vector({nu1(z, xy, ab), nu2(z, xy, ab), nu0(z, xy, ab)}); }
 
    static DenseMatrix grad_nu01(real_t z, Vector xy, unsigned int ab);
    static DenseMatrix grad_nu012(real_t z, Vector xy, unsigned int ab);
 
    static Vector nu01_grad_nu01(real_t z, Vector xy, unsigned int ab);
+   static Vector nu12_grad_nu12(real_t z, Vector xy, unsigned int ab);
    static Vector nu012_grad_nu012(real_t z, Vector xy, unsigned int ab);
 
    /// Shifted and Scaled Legendre Polynomials
@@ -377,6 +389,8 @@ public:
             Vector t, const DenseMatrix &grad_t,
             DenseTensor &u, DenseTensor &curl_u) const;
 
+   void E_T(int p, Vector s, Vector sds, DenseTensor &u) const;
+
    /** This is a vector-valued function associated with the quadrilateral face
        of a pyramid