diff --git a/fem/fe/fe_base.cpp b/fem/fe/fe_base.cpp
index 521a8d9d0b..6d971225b2 100644
--- a/fem/fe/fe_base.cpp
+++ b/fem/fe/fe_base.cpp
@@ -2168,6 +2168,64 @@ void Poly_1D::CalcDBinomTerms(const int p, const real_t x, const real_t y,
    }
 }
 
+void Poly_1D::CalcDxBinomTerms(const int p, const real_t x, const real_t y,
+                               real_t *u)
+{
+   if (p == 0)
+   {
+      u[0] = 0.;
+   }
+   else
+   {
+      int i;
+      const int *b = Binom(p);
+      real_t z = 1.;
+
+      for (i = 1; i < p; i++)
+      {
+         u[i] = i * b[i]*z;
+         z *= x;
+      }
+      u[p] = i * z;
+      z = y;
+      for (i--; i > 0; i--)
+      {
+         u[i] *= z;
+         z *= y;
+      }
+      u[0] = 0;
+   }
+}
+
+void Poly_1D::CalcDyBinomTerms(const int p, const real_t x, const real_t y,
+                               real_t *u)
+{
+   if (p == 0)
+   {
+      u[0] = 0.;
+   }
+   else
+   {
+      int i;
+      const int *b = Binom(p);
+      real_t z = x;
+
+      for (i = 1; i < p; i++)
+      {
+         u[i] = b[i]*z;
+         z *= x;
+      }
+      u[p] = 0.;
+      z = 1.;
+      for (i--; i > 0; i--)
+      {
+         u[i] *= (p - i) * z;
+         z *= y;
+      }
+      u[0] = p * z;
+   }
+}
+
 void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u)
 {
    // use the recursive definition for [-1,1]:
diff --git a/fem/fe/fe_base.hpp b/fem/fe/fe_base.hpp
index 1ee7aacc3d..9174dff8db 100644
--- a/fem/fe/fe_base.hpp
+++ b/fem/fe/fe_base.hpp
@@ -1159,6 +1159,16 @@ public:
        in the already allocated @a d array.*/
    static void CalcDBinomTerms(const int p, const real_t x, const real_t y,
                                real_t *d);
+   /** @brief Compute the derivatives (w.r.t. x) of the terms in the expansion
+       of the binomial (x + y)^p.  Store the results in the already allocated
+       @a d array.*/
+   static void CalcDxBinomTerms(const int p, const real_t x, const real_t y,
+                                real_t *d);
+   /** @brief Compute the derivatives (w.r.t. y) of the terms in the expansion
+       of the binomial (x + y)^p.  Store the results in the already allocated
+       @a d array.*/
+   static void CalcDyBinomTerms(const int p, const real_t x, const real_t y,
+                                real_t *d);
 
    /** @brief Compute the values of the Bernstein basis functions of order
        @a p at coordinate @a x and store the results in the already allocated