add margin to WeightRevDiff test

minor compilation error fixes

fem/coefficient.cpp

@@ -14,6 +14,7 @@
 #include "fem.hpp"
 
 #include <cmath>
+#include <cstddef>
 #include <limits>
 
 namespace mfem
@@ -48,6 +49,15 @@ ElementTransformation *RefinedToCoarse(
    return coarse_T;
 }
 
+void Coefficient::EvalRevDiff(const double Q_bar,
+                              ElementTransformation &T,
+                              const IntegrationPoint &ip,
+                              DenseMatrix &PointMat_bar)
+{
+   MFEM_ABORT("Coefficient::EvalRevDiff\n"
+              "\tEvalRevDiff not implemented for this coefficient!\n");
+}
+
 double PWConstCoefficient::Eval(ElementTransformation & T,
                                 const IntegrationPoint & ip)
 {
@@ -119,6 +129,34 @@ double FunctionCoefficient::Eval(ElementTransformation & T,
    }
 }
 
+void FunctionCoefficient::EvalRevDiff(const double Q_bar,
+                                      ElementTransformation &T,
+                                      const IntegrationPoint &ip,
+                                      DenseMatrix &PointMat_bar)
+{
+   int space_dim = T.GetSpaceDim();
+   double x[3] = {};
+   Vector transip(x, space_dim);
+   T.Transform(ip, transip);
+
+   double x_bar[3] = {};
+   Vector transip_bar(x_bar, space_dim);
+   if (Function)
+   {
+      MFEM_ASSERT(FunctionRevDiff, "EvalRevDiff: reverse-mode differentiated "
+                  "version of Function must be provided");
+      FunctionRevDiff(transip, Q_bar, transip_bar);
+   }
+   else
+   {
+      MFEM_ASSERT(TDFunctionRevDiff, "EvalRevDiff: reverse-mode differentiated"
+                  " version of TDFunction must be provided");
+      TDFunctionRevDiff(transip, GetTime(), Q_bar, transip_bar);
+   }
+   static_cast<IsoparametricTransformation &>(T).TransformRevDiff(
+      ip, transip_bar, PointMat_bar);
+}
+
 double GridFunctionCoefficient::Eval (ElementTransformation &T,
                                       const IntegrationPoint &ip)
 {
@@ -189,6 +227,15 @@ void RestrictedCoefficient::SetTime(double t)
    this->Coefficient::SetTime(t);
 }
 
+void VectorCoefficient::EvalRevDiff(const Vector &V_bar,
+                                    ElementTransformation &T,
+                                    const IntegrationPoint &ip,
+                                    DenseMatrix &PointMat_bar)
+{
+   MFEM_ABORT("VectorCoefficient::EvalRevDiff\n"
+              "\tEvalRevDiff not implemented for this coefficient!\n");
+}
+
 void VectorCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                              const IntegrationRule &ir)
 {
@@ -283,6 +330,35 @@ void VectorFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
    }
 }
 
+void VectorFunctionCoefficient::EvalRevDiff(const Vector &V_bar,
+                                            ElementTransformation &T,
+                                            const IntegrationPoint &ip,
+                                            DenseMatrix &PointMat_bar)
+{
+   MFEM_ASSERT( Q == NULL, "EvalRevDiff: not implemented for use with Q.")
+
+   double x[3];
+   Vector transip(x, vdim);
+   double x_bar[3];
+   Vector transip_bar(x_bar, vdim);
+   T.Transform(ip, transip);
+   transip_bar = 0.0;
+   if (Function)
+   {
+      MFEM_ASSERT(FunctionRevDiff, "EvalRevDiff: reverse-mode differentiated "
+                  "version of Function must be provided");
+      FunctionRevDiff(transip, V_bar, transip_bar);
+   }
+   else
+   {
+      MFEM_ASSERT(TDFunctionRevDiff, "EvalRevDiff: reverse-mode differentiated"
+                  " version of TDFunction must be provided");
+      TDFunctionRevDiff(transip, GetTime(), V_bar, transip_bar);
+   }
+   static_cast<IsoparametricTransformation &>(T).TransformRevDiff(
+      ip, transip_bar, PointMat_bar);
+}
+
 VectorArrayCoefficient::VectorArrayCoefficient (int dim)
    : VectorCoefficient(dim), Coeff(dim), ownCoeff(dim)
 {
@@ -797,6 +873,25 @@ void ProductCoefficient::SetTime(double t)
    this->Coefficient::SetTime(t);
 }
 
+void ProductCoefficient::EvalRevDiff(const double Q_bar,
+                                     ElementTransformation &T,
+                                     const IntegrationPoint &ip,
+                                     DenseMatrix &PointMat_bar)
+{
+   if (a == nullptr)
+   {
+      b->EvalRevDiff(Q_bar * aConst, T, ip, PointMat_bar);
+   }
+   else
+   {
+      double a_val = a->Eval(T, ip);
+      double b_val = b->Eval(T, ip);
+
+      a->EvalRevDiff(Q_bar * b_val, T, ip, PointMat_bar);
+      b->EvalRevDiff(Q_bar * a_val, T, ip, PointMat_bar);
+   }
+}
+
 void RatioCoefficient::SetTime(double t)
 {
    if (a) { a->SetTime(t); }
@@ -967,6 +1062,35 @@ void ScalarVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
    V *= sa;
 }
 
+void ScalarVectorProductCoefficient::EvalRevDiff(
+   const Vector &V_bar,
+   ElementTransformation &T,
+   const IntegrationPoint &ip,
+   DenseMatrix &PointMat_bar)
+{
+#ifdef MFEM_THREAD_SAFE
+   Vector W(V_bar.Size());
+   Vector W_bar(V_bar.Size());
+#else
+   W.SetSize(V_bar.Size());
+   W_bar.SetSize(V_bar.Size());
+#endif
+
+   double sa = (a == nullptr) ? aConst : a->Eval(T, ip);
+   b->Eval(W, T, ip);
+   W *= sa;
+
+   /// reverse pass
+   W_bar = 0.0;
+   add(W_bar, sa, V_bar, W_bar);
+   b->EvalRevDiff(W_bar, T, ip, PointMat_bar);
+   if (a != nullptr)
+   {
+      const double sa_bar = V_bar * W;
+      a->EvalRevDiff(sa_bar, T, ip, PointMat_bar);
+   }
+}
+
 NormalizedVectorCoefficient::NormalizedVectorCoefficient(VectorCoefficient &A,
                                                          double tol_)
    : VectorCoefficient(A.GetVDim()), a(&A), tol(tol_)

fem/coefficient.hpp

@@ -58,6 +58,21 @@ public:
    virtual double Eval(ElementTransformation &T,
                        const IntegrationPoint &ip) = 0;
 
+   /** @brief Reverse-mode differentiation of Eval w.r.t. the mesh node
+       locations in the element described by @a T, accumulating the result in
+       @a PointMat_bar */
+   /** @param[in] Q_bar - derivative of some output w.r.t. result of Eval */
+   /** @param[in] T - an element transformation */
+   /** @param[in] ip - defines location in reference space */
+   /** @param[inout] PointMat_bar - derivative of output w.r.t. mesh nodes */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void EvalRevDiff(const double Q_bar,
+                            ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar);
+
    /** @brief Evaluate the coefficient in the element described by @a T at the
        point @a ip at time @a t. */
    /** @note When this method is called, the caller must make sure that the
@@ -211,7 +226,15 @@ class FunctionCoefficient : public Coefficient
 {
 protected:
    std::function<double(const Vector &)> Function;
+   std::function<void(const Vector &,
+                      const double,
+                      Vector &)> FunctionRevDiff;
+
    std::function<double(const Vector &, double)> TDFunction;
+   std::function<void(const Vector &,
+                      double,
+                      const double,
+                      Vector &)> TDFunctionRevDiff;
 
 public:
    /// Define a time-independent coefficient from a std function
@@ -226,6 +249,23 @@ public:
       : TDFunction(std::move(TDF))
    { }
 
+   /// Construct time-independent coefficient that can be differentiated
+   FunctionCoefficient(std::function<double(const Vector &)> F,
+                       std::function<void(const Vector &,
+                                          const double,
+                                          Vector &)> dF)
+      : Function(F), FunctionRevDiff(dF)
+   { }
+
+   /// Construct time-dependent coefficient that can be differentiated
+   FunctionCoefficient(std::function<double(const Vector &, double)> TDF,
+                       std::function<void(const Vector &,
+                                          double,
+                                          const double,
+                                          Vector &)> dTDF)
+      : TDFunction(TDF), TDFunctionRevDiff(dTDF)
+   { }
+
    /// (DEPRECATED) Define a time-independent coefficient from a C-function
    /** @deprecated Use the method where the C-function, @a f, uses a const
        Vector argument instead of Vector. */
@@ -247,6 +287,11 @@ public:
    /// Evaluate the coefficient at @a ip.
    virtual double Eval(ElementTransformation &T,
                        const IntegrationPoint &ip);
+
+   virtual void EvalRevDiff(const double Q_bar,
+                            ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar);
 };
 
 class GridFunction;
@@ -456,6 +501,20 @@ public:
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip) = 0;
 
+   /** @brief Reverse-mode differentiation of Eval w.r.t. the mesh node
+       locations in the element described by @a T, accumulating the result in
+       @a PointMat_bar */
+   /** @param[in] V_bar - derivative of some output with respect to `V` */
+   /** @param[in] T - an element transformation */
+   /** @param[in] ip - defines location in reference space */
+   /** @param[inout] PointMat_bar - derivative of output w.r.t. mesh nodes */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void EvalRevDiff(const Vector &V_bar, ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar);
+
    /** @brief Evaluate the vector coefficient in the element described by @a T
        at all points of @a ir, storing the result in @a M. */
    /** The dimensions of @a M are GetVDim() by ir.GetNPoints() and they must be
@@ -490,6 +549,10 @@ public:
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip) { V = vec; }
 
+   virtual void EvalRevDiff(const Vector &V_bar, ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar) { }
+
    /// Return a reference to the constant vector in this class.
    const Vector& GetVec() { return vec; }
 };
@@ -582,8 +645,17 @@ class VectorFunctionCoefficient : public VectorCoefficient
 {
 private:
    std::function<void(const Vector &, Vector &)> Function;
+   std::function<void(const Vector &,
+                      const Vector &,
+                      Vector &)> FunctionRevDiff;
+
    std::function<void(const Vector &, double, Vector &)> TDFunction;
+   std::function<void(const Vector &,
+                      double,
+                      const Vector &,
+                      Vector &)> TDFunctionRevDiff;
    Coefficient *Q;
+   // Coefficient *dQ;
 
 public:
    /// Define a time-independent vector coefficient from a std function
@@ -606,11 +678,37 @@ public:
       : VectorCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
    { }
 
+   /// Construct time-independent vector coefficient that can be differentiated
+   VectorFunctionCoefficient(int dim,
+                             std::function<void(const Vector &,
+                                                Vector &)> F,
+                             std::function<void(const Vector &,
+                                                const Vector &,
+                                                Vector &)> dF)
+      : VectorCoefficient(dim), Function(std::move(F)),
+        FunctionRevDiff(std::move(dF)), Q(NULL)
+   { }
+
+   /// Construct time-dependent vector coefficient that can be differentiated
+   VectorFunctionCoefficient(int dim,
+                             std::function<void(const Vector &,
+                                                double,
+                                                Vector &)> TDF,
+                             std::function<void(const Vector &,
+                                                double, const Vector &, Vector &)> dTDF)
+      : VectorCoefficient(dim), TDFunction(std::move(TDF)),
+        TDFunctionRevDiff(std::move(dTDF)), Q(NULL)
+   { }
+
    using VectorCoefficient::Eval;
    /// Evaluate the vector coefficient at @a ip.
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip);
 
+   virtual void EvalRevDiff(const Vector &V_bar, ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar);
+
    virtual ~VectorFunctionCoefficient() { }
 };
 
@@ -1388,6 +1486,11 @@ public:
    virtual double Eval(ElementTransformation &T,
                        const IntegrationPoint &ip)
    { return ((a == NULL ) ? aConst : a->Eval(T, ip) ) * b->Eval(T, ip); }
+
+   void EvalRevDiff(const double Q_bar,
+                    ElementTransformation &T,
+                    const IntegrationPoint &ip,
+                    DenseMatrix &PointMat_bar);
 };
 
 /** @brief Scalar coefficient defined as the ratio of two scalars where one or
@@ -1657,6 +1760,9 @@ private:
    double aConst;
    Coefficient * a;
    VectorCoefficient * b;
+#ifndef MFEM_THREAD_SAFE
+   Vector W, W_bar;
+#endif
 
 public:
    /// Constructor with constant and vector coefficient.  Result is A * B.
@@ -1687,6 +1793,11 @@ public:
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip);
    using VectorCoefficient::Eval;
+
+   virtual void EvalRevDiff(const Vector &V_bar, ElementTransformation &T,
+                            const IntegrationPoint &ip,
+                            DenseMatrix &PointMat_bar);
+
 };
 
 /// Vector coefficient defined as a normalized vector field (returns v/|v|)

fem/eltrans.cpp

@@ -532,6 +532,82 @@ void IsoparametricTransformation::Transform (const DenseMatrix &matrix,
    }
 }
 
+void IsoparametricTransformation::TransformRevDiff(const IntegrationPoint &ip,
+                                                   const Vector &x_bar,
+                                                   DenseMatrix &PointMat_bar)
+{
+   MFEM_ASSERT((PointMat_bar.Width() == PointMat.Width()) &&
+               (PointMat_bar.Height() == PointMat.Height()),
+               "PointMat_bar shape != PointMat shape");
+   shape.SetSize(FElem->GetDof());
+   FElem->CalcShape(ip, shape);
+   AddMultVWt(x_bar, shape, PointMat_bar);
+}
+
+void IsoparametricTransformation::JacobianRevDiff(const DenseMatrix &dFdx_bar,
+                                                  DenseMatrix &PointMat_bar)
+{
+   MFEM_ASSERT((PointMat_bar.Width() == PointMat.Width()) &&
+               (PointMat_bar.Height() == PointMat.Height()),
+               "PointMat_bar shape != PointMat shape");
+
+   dshape.SetSize(FElem->GetDof(), FElem->GetDim());
+   if (dshape.Width() > 0)
+   {
+      // The math here can be found in Giles' report "An extended collection of
+      // matrix derivative results for forward and reverse mode algorithmic
+      // differentiation"
+      FElem->CalcDShape(*IntPoint, dshape);
+      AddMultABt(dFdx_bar, dshape, PointMat_bar);
+   }
+}
+
+void IsoparametricTransformation::AdjugateJacobianRevDiff(
+   const DenseMatrix &adjJ_bar, DenseMatrix &PointMat_bar)
+{
+   Jacobian(); // Recompute the Jacobian, if necessary
+   double dFdx_bar_buffer[9];
+   DenseMatrix dFdx_bar(dFdx_bar_buffer, dFdx.Height(), dFdx.Width());
+   if (dFdx.Width() > 0)
+   {
+      CalcAdjugateRevDiff(dFdx, adjJ_bar, dFdx_bar);
+   }
+   JacobianRevDiff(dFdx_bar, PointMat_bar);
+}
+
+void IsoparametricTransformation::InverseJacobianRevDiff(
+   const DenseMatrix &invJ_bar, DenseMatrix &PointMat_bar)
+{
+   Jacobian(); // Recompute the Jacobian, if necessary
+   double dFdx_bar_buffer[9];
+   DenseMatrix dFdx_bar(dFdx_bar_buffer, dFdx.Height(), dFdx.Width());
+   if (dFdx.Width() > 0)
+   {
+      CalcInverseRevDiff(dFdx, invJ_bar, dFdx_bar);
+   }
+   JacobianRevDiff(dFdx_bar, PointMat_bar);
+}
+
+void IsoparametricTransformation::WeightRevDiff(DenseMatrix &PointMat_bar)
+{
+   Jacobian(); // Recompute the Jacobian, if necessary
+   double dFdx_bar_buffer[9];
+   DenseMatrix dFdx_bar(dFdx_bar_buffer, dFdx.Height(), dFdx.Width());
+   dFdx.WeightRevDiff(dFdx_bar);
+   JacobianRevDiff(dFdx_bar, PointMat_bar);
+}
+
+void IsoparametricTransformation::WeightRevDiff(double weight_bar,
+                                                DenseMatrix &PointMat_bar)
+{
+   Jacobian(); // Recompute the Jacobian, if necessary
+   double dFdx_bar_buffer[9];
+   DenseMatrix dFdx_bar(dFdx_bar_buffer, dFdx.Height(), dFdx.Width());
+   dFdx.WeightRevDiff(dFdx_bar);
+   dFdx_bar *= weight_bar;
+   JacobianRevDiff(dFdx_bar, PointMat_bar);
+}
+
 void IntegrationPointTransformation::Transform (const IntegrationPoint &ip1,
                                                 IntegrationPoint &ip2)
 {

fem/eltrans.hpp

@@ -446,6 +446,58 @@ public:
       return inv_tr.Transform(v, ip);
    }
 
+   /// @brief Reverse-mode differentiation of Transform() w.r.t. PointMat
+   /// @param[in] ip - specifies the location in reference space
+   /// @param[in] x_bar - derivative of some output w.r.t. x
+   /// @param[out] PointMat_bar - derivative of output w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void TransformRevDiff(const IntegrationPoint &ip, const Vector &x_bar,
+                         DenseMatrix &PointMat_bar);
+
+   /// @brief Reverse-mode differentiation of Jacobian() w.r.t. PointMat
+   /// @param[in] dFdx_bar - derivative of functional w.r.t. Jacobian
+   /// @param[out] PointMat_bar - derivative w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void JacobianRevDiff(const DenseMatrix &dFdx_bar,
+                        DenseMatrix &PointMat_bar);
+
+   /// @brief Reverse-mode differentiation of AdjugateJacobian() w.r.t. PointMat
+   /// @param[in] adjJ_bar - derivative of functional w.r.t. Adjugate
+   /// @param[out] PointMat_bar - derivative w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void AdjugateJacobianRevDiff(const DenseMatrix &adjJ_bar,
+                                DenseMatrix &PointMat_bar);
+
+   /// @brief Reverse-mode differentiation of InverseJacobian() w.r.t PointMat
+   /// @param[in] invJ_bar - derivative of functional w.r.t. Inverse
+   /// @param[out] PointMat_bar - derivative w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void InverseJacobianRevDiff(const DenseMatrix &adjJ_bar,
+                               DenseMatrix &PointMat_bar);
+
+   /// @brief Reverse-mode differentiation of Weight()
+   /// @param[out] PointMat_bar - derivative of functional w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void WeightRevDiff(DenseMatrix &PointMat_bar);
+
+   /// @brief Reverse-mode differentiation of Weight()
+   /// @param[in] weight_bar - derivative of functional w.r.t Weight
+   /// @param[out] PointMat_bar - derivative of functional w.r.t. PointMat
+   /// @note PointMat_bar must have the same shape as PointMat
+   /// @warning This routine does not initialize PointMat_bar, and instead
+   /// accumulates (with += or -=) contributions to its derivative.
+   void WeightRevDiff(double weight_bar, DenseMatrix &PointMat_bar);
+
    virtual ~IsoparametricTransformation() { }
 
    MFEM_DEPRECATED void FinalizeTransformation() {}

fem/fe/fe_base.cpp

@@ -30,6 +30,7 @@ FiniteElement::FiniteElement(int D, Geometry::Type G, int Do, int O, int F)
    deriv_map_type = VALUE;
    for (int i = 0; i < Geometry::MaxDim; i++) { orders[i] = -1; }
 #ifndef MFEM_THREAD_SAFE
+   shape.SetSize(dof);
    vshape.SetSize(dof, dim);
 #endif
 }
@@ -46,6 +47,13 @@ void FiniteElement::CalcVShape (
    MFEM_ABORT("method is not implemented for this class");
 }
 
+void FiniteElement::CalcVShapeRevDiff(ElementTransformation &Trans,
+                                      const DenseMatrix &shape_bar,
+                                      DenseMatrix &PointMat_bar) const
+{
+   MFEM_ABORT("method is not implemented for this class");
+}
+
 void FiniteElement::CalcDivShape (
    const IntegrationPoint &ip, Vector &divshape) const
 {
@@ -90,6 +98,60 @@ void FiniteElement::CalcPhysCurlShape(ElementTransformation &Trans,
    }
 }
 
+void FiniteElement::CalcPhysCurlShapeRevDiff(ElementTransformation &Trans,
+                                             const DenseMatrix &curlshape_bar,
+                                             DenseMatrix &PointMat_bar) const
+{
+   switch (dim)
+   {
+      case 3:
+      {
+#ifdef MFEM_THREAD_SAFE
+         DenseMatrix vshape(dof, dim);
+#endif
+         DenseMatrix vshapedxt(dof, dim);
+         DenseMatrix vshapedxt_bar(dof, dim);
+
+         CalcCurlShape(Trans.GetIntPoint(), vshape);
+         const auto &jac = Trans.Jacobian();
+         MultABt(vshape, jac, vshapedxt);
+
+         const double weight = Trans.Weight();
+         // curl_shape *= 1.0 / weight;
+
+         /// start reverse pass
+         auto &isotrans = dynamic_cast<IsoparametricTransformation&>(Trans);
+
+         /// curl_shape = vshapedxt / weight;
+         double weight_bar = 0.0;
+         for (int j = 0; j < curlshape_bar.Width(); ++j)
+         {
+            for (int i = 0; i < curlshape_bar.Height(); ++i)
+            {
+               weight_bar -= curlshape_bar(i,j) * vshapedxt(i,j) / pow(weight, 2);
+            }
+         }
+         vshapedxt_bar = curlshape_bar; vshapedxt_bar *= (1.0 / weight);
+
+         /// double weight = Trans.Weight();
+         isotrans.WeightRevDiff(weight_bar, PointMat_bar);
+
+         /// const auto &jac = Trans.Jacobian();
+         /// MultABt(vshape, jac, vshapedxt);
+         double jac_bar_buffer[9];
+         DenseMatrix jac_bar(jac_bar_buffer, jac.Width(), jac.Height());
+         MultAtB(vshapedxt_bar, vshape, jac_bar);
+         isotrans.JacobianRevDiff(jac_bar, PointMat_bar);
+         break;
+      }
+      case 2:
+         MFEM_ABORT("CalcPhysCurlShapeRevDiff not implemented!\n");
+         break;
+      default:
+         MFEM_ABORT("Invalid dimension, Dim = " << dim);
+   }
+}
+
 void FiniteElement::GetFaceDofs(int face, int **dofs, int *ndofs) const
 {
    MFEM_ABORT("method is not overloaded");
@@ -138,6 +200,16 @@ void FiniteElement::ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
    mfem_error ("FiniteElement::ProjectFromNodes() (vector) is not overloaded!");
 }
 
+void FiniteElement::ProjectRevDiff (
+   const Vector &P_bar,
+   VectorCoefficient &vc,
+   ElementTransformation &Trans,
+   DenseMatrix &PointMat_bar) const
+{
+   mfem_error ("FiniteElement::ProjectRevDiff (...) (vector) is not "
+               "overloaded !");
+}
+
 void FiniteElement::ProjectMatrixCoefficient(
    MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
 {
@@ -186,6 +258,26 @@ void FiniteElement::CalcPhysShape(ElementTransformation &Trans,
    }
 }
 
+void FiniteElement::CalcPhysShapeRevDiff(ElementTransformation &Trans,
+                                         const Vector &shape_bar,
+                                         DenseMatrix &PointMat_bar) const
+{
+   if (map_type == INTEGRAL)
+   {
+#ifdef MFEM_THREAD_SAFE
+      Vector shape(dof);
+#endif
+      CalcShape(Trans.GetIntPoint(), shape);
+      // shape /= Trans.Weight();
+      auto weight = Trans.Weight();
+      auto weight_bar = -(shape_bar * shape) / pow(weight, 2);
+
+      // cast the ElementTransformation
+      auto &isotrans = dynamic_cast<IsoparametricTransformation &>(Trans);
+      isotrans.WeightRevDiff(weight_bar, PointMat_bar);
+   }
+}
+
 void FiniteElement::CalcPhysDShape(ElementTransformation &Trans,
                                    DenseMatrix &dshape) const
 {
@@ -894,6 +986,46 @@ void VectorFiniteElement::CalcVShape_RT (
    shape *= (1.0 / Trans.Weight());
 }
 
+void VectorFiniteElement::CalcVShape_RTRevDiff(ElementTransformation &Trans,
+                                               const DenseMatrix &vshape_bar,
+                                               DenseMatrix &PointMat_bar) const
+{
+   MFEM_ASSERT(map_type == H_DIV, "");
+#ifdef MFEM_THREAD_SAFE
+   DenseMatrix vshape(dof, dim);
+   DenseMatrix vshapedxt(dof, dim);
+   DenseMatrix vshapedxt_bar(dof, dim);
+#else
+   vshapedxt.SetSize(dof, dim);
+   vshapedxt_bar.SetSize(dof, dim);
+#endif
+   CalcVShape(Trans.GetIntPoint(), vshape);
+   const auto &jac = Trans.Jacobian();
+   MultABt(vshape, jac, vshapedxt);
+
+   const double weight = Trans.Weight();
+   // shape *= (1.0 / weight);
+
+   /// start reverse pass
+   auto &isotrans = dynamic_cast<IsoparametricTransformation&>(Trans);
+
+   double weight_bar = 0.0;
+   for (int j = 0; j < vshape_bar.Width(); ++j)
+   {
+      for (int i = 0; i < vshape_bar.Height(); ++i)
+      {
+         weight_bar -= vshape_bar(i,j) * vshapedxt(i,j) / pow(weight,2);
+      }
+   }
+   isotrans.WeightRevDiff(weight_bar, PointMat_bar);
+
+   vshapedxt_bar = vshape_bar; vshapedxt_bar *= (1.0 / weight);
+   double jac_bar_buffer[9];
+   DenseMatrix jac_bar(jac_bar_buffer, jac.Width(), jac.Height());
+   MultAtB(vshapedxt_bar, vshape, jac_bar);
+   isotrans.JacobianRevDiff(jac_bar, PointMat_bar);
+}
+
 void VectorFiniteElement::CalcVShape_ND (
    ElementTransformation &Trans, DenseMatrix &shape) const
 {
@@ -905,6 +1037,46 @@ void VectorFiniteElement::CalcVShape_ND (
    Mult(vshape, Trans.InverseJacobian(), shape);
 }
 
+void VectorFiniteElement::CalcVShape_NDRevDiff(ElementTransformation &Trans,
+                                               const DenseMatrix &vshape_bar,
+                                               DenseMatrix &PointMat_bar) const
+{
+   MFEM_ASSERT(map_type == H_CURL, "");
+#ifdef MFEM_THREAD_SAFE
+   DenseMatrix vshape(dof, dim);
+   DenseMatrix vshapedxt(dof, dim);
+   DenseMatrix vshapedxt_bar(dof, dim);
+#else
+   vshapedxt.SetSize(dof, dim);
+   vshapedxt_bar.SetSize(dof, dim);
+#endif
+   CalcVShape(Trans.GetIntPoint(), vshape);
+   const auto &adjJ = Trans.AdjugateJacobian();
+   Mult(vshape, adjJ, vshapedxt);
+
+   const double weight = Trans.Weight();
+   // shape *= (1.0 / weight);
+
+   /// start reverse pass
+   auto &isotrans = dynamic_cast<IsoparametricTransformation&>(Trans);
+
+   double weight_bar = 0.0;
+   for (int j = 0; j < vshape_bar.Width(); ++j)
+   {
+      for (int i = 0; i < vshape_bar.Height(); ++i)
+      {
+         weight_bar -= vshape_bar(i,j) * vshapedxt(i,j) / pow(weight, 2);
+      }
+   }
+   isotrans.WeightRevDiff(weight_bar, PointMat_bar);
+
+   vshapedxt_bar = vshape_bar; vshapedxt_bar *= (1.0 / weight);
+   double adjJ_bar_buffer[9];
+   DenseMatrix adjJ_bar(adjJ_bar_buffer, adjJ.Width(), adjJ.Height());
+   MultAtB(vshape, vshapedxt_bar, adjJ_bar);
+   isotrans.AdjugateJacobianRevDiff(adjJ_bar, PointMat_bar);
+}
+
 void VectorFiniteElement::Project_RT(
    const double *nk, const Array<int> &d2n,
    VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
@@ -942,6 +1114,47 @@ void VectorFiniteElement::Project_RT(
    }
 }
 
+void VectorFiniteElement::Project_RTRevDiff(
+   const Vector &P_bar,
+   const double *nk, const Array<int> &d2n,
+   VectorCoefficient &vc, ElementTransformation &Trans,
+   DenseMatrix &PointMat_bar) const
+{
+   double vk[Geometry::MaxDim];
+   const int sdim = Trans.GetSpaceDim();
+   MFEM_ASSERT(vc.GetVDim() == sdim, "");
+   Vector xk(vk, sdim);
+   MFEM_ASSERT(dim == sdim, "VectorFiniteElement::Project_RTRevDiff\n"
+               "\tOnly implemented if space dim == reference dim!\n");
+
+   DenseMatrix temp_bar(PointMat_bar.Height(), PointMat_bar.Width());
+
+   IsoparametricTransformation &isotrans =
+      dynamic_cast<IsoparametricTransformation&>(Trans);
+
+   for (int k = 0; k < dof; k++)
+   {
+      temp_bar = 0.0;
+      isotrans.SetIntPoint(&Nodes.IntPoint(k));
+      vc.Eval(xk, isotrans, Nodes.IntPoint(k));
+      // dof_k = nk^t adj(J) xk
+      const Vector nk_vec(const_cast<double*>(nk + d2n[k]*dim), sdim);
+      double adjJ_bar_buffer[Geometry::MaxDim*Geometry::MaxDim];
+      DenseMatrix adjJ_bar(adjJ_bar_buffer, dim, dim);
+      MultVWt(nk_vec, xk, adjJ_bar);
+      isotrans.AdjugateJacobianRevDiff(adjJ_bar, temp_bar);
+
+      double V_bar_buffer[Geometry::MaxDim];
+      Vector V_bar(V_bar_buffer, sdim);
+      isotrans.AdjugateJacobian().MultTranspose(nk_vec, V_bar);
+      vc.EvalRevDiff(V_bar, isotrans,
+                     Nodes.IntPoint(k), temp_bar);
+
+      temp_bar *= P_bar(k);
+      PointMat_bar += temp_bar;
+   }
+}
+
 void VectorFiniteElement::ProjectMatrixCoefficient_RT(
    const double *nk, const Array<int> &d2n,
    MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
@@ -1159,6 +1372,48 @@ void VectorFiniteElement::Project_ND(
    }
 }
 
+void VectorFiniteElement::Project_NDRevDiff(
+   const Vector &P_bar,
+   const double *tk, const Array<int> &d2t,
+   VectorCoefficient &vc, ElementTransformation &Trans,
+   DenseMatrix &PointMat_bar) const
+{
+   double vk[Geometry::MaxDim];
+   const int sdim = Trans.GetSpaceDim();
+   MFEM_ASSERT(vc.GetVDim() == sdim, "");
+   Vector xk(vk, sdim);
+   MFEM_ASSERT(dim == sdim, "VectorFiniteElement::Project_NDRevDiff\n"
+               "\tOnly implemented if space dim == reference dim!\n");
+
+   DenseMatrix temp_bar(PointMat_bar.Height(), PointMat_bar.Width());
+
+   IsoparametricTransformation &isotrans =
+      dynamic_cast<IsoparametricTransformation&>(Trans);
+
+   for (int k = 0; k < dof; k++)
+   {
+      temp_bar = 0.0;
+      isotrans.SetIntPoint(&Nodes.IntPoint(k));
+      vc.Eval(xk, isotrans, Nodes.IntPoint(k));
+      // dof_k = nk^t J xk
+      const Vector tk_vec(const_cast<double*>(tk + d2t[k]*dim), sdim);
+
+      double J_bar_buffer[Geometry::MaxDim*Geometry::MaxDim];
+      DenseMatrix J_bar(J_bar_buffer, dim, dim);
+      MultVWt(xk, tk_vec, J_bar);
+      isotrans.JacobianRevDiff(J_bar, temp_bar);
+
+      double V_bar_buffer[Geometry::MaxDim];
+      Vector V_bar(V_bar_buffer, sdim);
+      isotrans.Jacobian().Mult(tk_vec, V_bar);
+      vc.EvalRevDiff(V_bar, isotrans,
+                     Nodes.IntPoint(k), temp_bar);
+
+      temp_bar *= P_bar(k);
+      PointMat_bar += temp_bar;
+   }
+}
+
 void VectorFiniteElement::ProjectMatrixCoefficient_ND(
    const double *tk, const Array<int> &d2t,
    MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const

fem/fe/fe_base.hpp

@@ -245,6 +245,7 @@ protected:
    mutable int orders[Geometry::MaxDim]; ///< Anisotropic orders
    IntegrationRule Nodes;
 #ifndef MFEM_THREAD_SAFE
+   mutable Vector shape;
    mutable DenseMatrix vshape; // Dof x Dim
 #endif
    /// Container for all DofToQuad objects created by the FiniteElement.
@@ -362,6 +363,10 @@ public:
    /** The size (#dof) of the result Vector @a shape must be set in advance. */
    void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const;
 
+   void CalcPhysShapeRevDiff(ElementTransformation &Trans,
+                             const Vector &shape_bar,
+                             DenseMatrix &PointMat_bar) const;
+
    /** @brief Evaluate the gradients of all shape functions of a scalar finite
        element in reference space at the given point @a ip. */
    /** Each row of the result DenseMatrix @a dshape contains the derivatives of
@@ -400,6 +405,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const;
 
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const;
+
    /// Equivalent to the CalcVShape() method with the same arguments.
    void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
    { CalcVShape(Trans, shape); }
@@ -435,6 +444,10 @@ public:
    void CalcPhysCurlShape(ElementTransformation &Trans,
                           DenseMatrix &curl_shape) const;
 
+   void CalcPhysCurlShapeRevDiff(ElementTransformation &Trans,
+                                 const DenseMatrix &curlshape_bar,
+                                 DenseMatrix &PointMat_bar) const;
+
    /** @brief Get the dofs associated with the given @a face.
        @a *dofs is set to an internal array of the local dofc on the
        face, while *ndofs is set to the number of dofs on that face.
@@ -519,6 +532,14 @@ public:
    virtual void Project(VectorCoefficient &vc,
                         ElementTransformation &Trans, Vector &dofs) const;
 
+   /** Given a vector coefficient and a transformation, compute the derivative of
+          its projection (approximation) in the local finite dimensional space
+          w.r.t. the mesh nodes (VectorFiniteElements) */
+   virtual void ProjectRevDiff(const Vector &P_bar,
+                               VectorCoefficient &vc,
+                               ElementTransformation &Trans,
+                               DenseMatrix &PointMat_bar) const;
+
    /** @brief Given a vector of values at the finite element nodes and a
        transformation, compute its projection (approximation) in the local
        finite dimensional space in terms of the degrees of freedom. Valid for
@@ -791,15 +812,24 @@ protected:
 #ifndef MFEM_THREAD_SAFE
    mutable DenseMatrix J, Jinv;
    mutable DenseMatrix curlshape, curlshape_J;
+   mutable DenseMatrix vshapedxt, vshapedxt_bar;
 #endif
    void SetDerivMembers();
 
    void CalcVShape_RT(ElementTransformation &Trans,
                       DenseMatrix &shape) const;
 
+   void CalcVShape_RTRevDiff(ElementTransformation &Trans,
+                             const DenseMatrix &shape_bar,
+                             DenseMatrix &PointMat_bar) const;
+
    void CalcVShape_ND(ElementTransformation &Trans,
                       DenseMatrix &shape) const;
 
+   void CalcVShape_NDRevDiff(ElementTransformation &Trans,
+                             const DenseMatrix &shape_bar,
+                             DenseMatrix &PointMat_bar) const;
+
    /** @brief Project a vector coefficient onto the RT basis functions
        @param nk    Face normal vectors for this element type
        @param d2n   Offset into nk for each degree of freedom
@@ -845,6 +875,23 @@ protected:
                    const FiniteElement &fe, ElementTransformation &Trans,
                    DenseMatrix &I) const;
 
+   /** Reverse-mode differentiation of Project_ND w.r.t. the mesh node
+       locations in the element described by @a T
+       @param[in] P_bar - derivative of function with respect to the projection
+       @param[in] nk - Face normal vectors for this element type
+       @param[in] d2n - Offset into nk for each degree of freedom
+       @param[in] vc - VectorCoefficient being projected
+       @param[in] Trans - an element transformation
+       @param[out] PointMat_bar - derivative of projected degrees of freedom w.r.t.
+                              mesh nodes
+       @warning - only implemented for the same space and reference dimension
+   */
+   void Project_RTRevDiff(const Vector &P_bar,
+                          const double *nk, const Array<int> &d2n,
+                          VectorCoefficient &vc,
+                          ElementTransformation &Trans,
+                          DenseMatrix &PointMat_bar) const;
+
    // rotated gradient in 2D
    void ProjectGrad_RT(const double *nk, const Array<int> &d2n,
                        const FiniteElement &fe, ElementTransformation &Trans,
@@ -883,6 +930,22 @@ protected:
                    Vector &vc, ElementTransformation &Trans,
                    Vector &dofs) const;
 
+   /** Reverse-mode differentiation of Project_ND w.r.t. the mesh node
+       locations in the element described by @a T
+       @param[in] P_bar - derivative of output with respect to the projection
+       @param[in] tk - Edge tangent vectors for this element type
+       @param[in] d2t - Offset into tk for each degree of freedom
+       @param[in] vc - Vector coefficient being projected
+       @param[in] Trans - Transformation from reference to physical coordinates
+       @param[out] PointMat_bar - derivative of some output w.r.t. mesh nodes
+       @warning - only implemented for the same space and reference dimension
+   */
+   void Project_NDRevDiff(const Vector &P_bar,
+                          const double *tk, const Array<int> &d2t,
+                          VectorCoefficient &vc,
+                          ElementTransformation &Trans,
+                          DenseMatrix &PointMat_bar) const;
+
    /// Project the rows of the matrix coefficient in an ND space
    void ProjectMatrixCoefficient_ND(
       const double *tk, const Array<int> &d2t,

fem/fe/fe_nd.hpp

@@ -42,6 +42,11 @@ public:
                            DenseMatrix &shape) const
    { CalcVShape_ND(Trans, shape); }
 
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_NDRevDiff(Trans, shape_bar, PointMat_bar); }
+
    virtual void CalcCurlShape(const IntegrationPoint &ip,
                               DenseMatrix &curl_shape) const;
 
@@ -120,6 +125,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_ND(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_NDRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcCurlShape(const IntegrationPoint &ip,
                               DenseMatrix &curl_shape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -139,6 +148,11 @@ public:
       if (obasis1d.IsIntegratedType()) { ProjectIntegrated(vc, Trans, dofs); }
       else { Project_ND(tk, dof2tk, vc, Trans, dofs); }
    }
+   virtual void ProjectRevDiff(const Vector &P_bar,
+                               VectorCoefficient &vc,
+                               ElementTransformation &Trans,
+                               DenseMatrix &dofs_bar) const
+   { Project_NDRevDiff(P_bar, tk, dof2tk, vc, Trans, dofs_bar); }
    virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
                                  Vector &dofs) const
    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
@@ -182,6 +196,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_ND(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_NDRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcCurlShape(const IntegrationPoint &ip,
                               DenseMatrix &curl_shape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -198,6 +216,11 @@ public:
    virtual void Project(VectorCoefficient &vc,
                         ElementTransformation &Trans, Vector &dofs) const
    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
+   virtual void ProjectRevDiff(const Vector &P_bar,
+                               VectorCoefficient &vc,
+                               ElementTransformation &Trans,
+                               DenseMatrix &dofs_bar) const
+   { Project_NDRevDiff(P_bar, tk, dof2tk, vc, Trans, dofs_bar); }
    virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
                                  Vector &dofs) const
    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
@@ -241,6 +264,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_ND(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_NDRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcCurlShape(const IntegrationPoint &ip,
                               DenseMatrix &curl_shape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -292,6 +319,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_ND(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_NDRevDiff(Trans, shape_bar, PointMat_bar); }
    // virtual void CalcCurlShape(const IntegrationPoint &ip,
    //                            DenseMatrix &curl_shape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,

fem/fe/fe_rt.hpp

@@ -41,6 +41,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_RT(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_RTRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -60,6 +64,11 @@ public:
       if (obasis1d.IsIntegratedType()) { ProjectIntegrated(vc, Trans, dofs); }
       else { Project_RT(nk, dof2nk, vc, Trans, dofs); }
    }
+   virtual void ProjectRevDiff(const Vector &P_bar,
+                               VectorCoefficient &vc,
+                               ElementTransformation &Trans,
+                               DenseMatrix &dofs_bar) const
+   { Project_RTRevDiff(P_bar, nk, dof2nk, vc, Trans, dofs_bar); }
    virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
                                  Vector &dofs) const
    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
@@ -110,6 +119,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_RT(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_RTRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -172,6 +185,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_RT(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_RTRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -232,6 +249,10 @@ public:
    virtual void CalcVShape(ElementTransformation &Trans,
                            DenseMatrix &shape) const
    { CalcVShape_RT(Trans, shape); }
+   virtual void CalcVShapeRevDiff(ElementTransformation &Trans,
+                                  const DenseMatrix &shape_bar,
+                                  DenseMatrix &PointMat_bar) const
+   { CalcVShape_RTRevDiff(Trans, shape_bar, PointMat_bar); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
@@ -248,6 +269,11 @@ public:
    virtual void Project(VectorCoefficient &vc,
                         ElementTransformation &Trans, Vector &dofs) const
    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
+   virtual void ProjectRevDiff(const Vector &P_bar,
+                               VectorCoefficient &vc,
+                               ElementTransformation &Trans,
+                               DenseMatrix &dofs_bar) const
+   { Project_RTRevDiff(P_bar, nk, dof2nk, vc, Trans, dofs_bar); }
    virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
                                  Vector &dofs) const
    { Project_RT(nk, dof2nk, vc, Trans, dofs); }

linalg/densemat.cpp

@@ -513,6 +513,119 @@ double DenseMatrix::Weight() const
    return 0.0;
 }
 
+void DenseMatrix::DetRevDiff(DenseMatrix &A_bar) const
+{
+   MFEM_ASSERT(Height() == Width() && Height() > 0,
+               "The matrix must be square and "
+               << "sized larger than zero to compute the determinant."
+               << "  Height() = " << Height()
+               << ", Width() = " << Width());
+
+   switch (Height())
+   {
+      case 1:
+         // return data[0];
+         A_bar(0,0) = 1.0;
+         return;
+
+      case 2:
+         // return data[0] * data[3] - data[1] * data[2];
+         A_bar(0,0) = data[3]; // data[0]
+         A_bar(1,1) = data[0]; // data[3]
+         A_bar(1,0) = -data[2]; // data[1]
+         A_bar(0,1) = -data[1]; // data[2]
+         return;
+
+      case 3:
+      {
+         const double *d = data;
+         // return
+         //   d[0] * (d[4] * d[8] - d[5] * d[7]) +
+         //   d[3] * (d[2] * d[7] - d[1] * d[8]) +
+         //   d[6] * (d[1] * d[5] - d[2] * d[4]);
+         A_bar(0,0) = d[4]*d[8] - d[5]*d[7]; // d[0]
+         A_bar(1,0) = d[6]*d[5] - d[3]*d[8]; // d[1]
+         A_bar(2,0) = d[3]*d[7] - d[6]*d[4]; // d[2]
+         A_bar(0,1) = d[2]*d[7] - d[1]*d[8]; // d[3]
+         A_bar(1,1) = d[0]*d[8] - d[6]*d[2]; // d[4]
+         A_bar(2,1) = d[6]*d[1] - d[0]*d[7]; // d[5]
+         A_bar(0,2) = d[1]*d[5] - d[2]*d[4]; // d[6]
+         A_bar(1,2) = d[3]*d[2] - d[0]*d[5]; // d[7]
+         A_bar(2,2) = d[0]*d[4] - d[3]*d[1]; // d[8]
+         return;
+
+      }
+      default:
+      {
+         // In the general case we compute the gradient of the determinant
+         // using the relation from Mike Giles document:
+         // "An extended collection of matrix derivative results for forward
+         // and reverse mode algorithmic differentiation"
+         DenseMatrixInverse lu_factors(*this);
+         lu_factors.GetInverseMatrix(A_bar);
+         A_bar.Transpose();
+         A_bar *= lu_factors.Det();
+         return;
+
+      }
+   }
+   // not reachable
+}
+
+void DenseMatrix::WeightRevDiff(DenseMatrix &A_bar) const
+{
+#ifdef MFEM_DEBUG
+   if (Height() != A_bar.Height() || Width() != A_bar.Width())
+   {
+      mfem_error("DenseMatrix::WeightRevDiff()");
+   }
+#endif
+   if (Height() == Width())
+   {
+      // return Det();
+      DetRevDiff(A_bar);
+      return;
+   }
+   else if ((Height() == 2) && (Width() == 1))
+   {
+      // return sqrt(data[0] * data[0] + data[1] * data[1]);
+      double wgt = sqrt(data[0] * data[0] + data[1] * data[1]);
+      A_bar(0,0) = data[0]/wgt;
+      A_bar(1,0) = data[1]/wgt;
+      return;
+   }
+   else if ((Height() == 3) && (Width() == 1))
+   {
+      // return sqrt(data[0] * data[0] + data[1] * data[1] + data[2] * data[2]);
+      double wgt = sqrt(data[0] * data[0] + data[1] * data[1] + data[2] * data[2]);
+      for (int i = 0; i < 3; ++i)
+      {
+         A_bar(i,0) = data[i]/wgt;
+      }
+      return;
+   }
+   else if ((Height() == 3) && (Width() == 2))
+   {
+      const double *d = data;
+      double E = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
+      double G = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
+      double F = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
+      double wgt = sqrt(E * G - F * F);
+      // start reverse sweep
+      double E_bar = 0.5*G/wgt;
+      double G_bar = 0.5*E/wgt;
+      double F_bar = -F/wgt;
+      A_bar(0,0) = F_bar*d[3] + 2.0*E_bar*d[0]; // d[0]
+      A_bar(1,0) = F_bar*d[4] + 2.0*E_bar*d[1]; // d[1]
+      A_bar(2,0) = F_bar*d[5] + 2.0*E_bar*d[2]; // d[2]
+      A_bar(0,1) = F_bar*d[0] + 2.0*G_bar*d[3]; // d[3]
+      A_bar(1,1) = F_bar*d[1] + 2.0*G_bar*d[4]; // d[4]
+      A_bar(2,1) = F_bar*d[2] + 2.0*G_bar*d[5]; // d[5]
+      return;
+   }
+   mfem_error("DenseMatrix::WeightRevDiff()");
+}
+
 void DenseMatrix::Set(double alpha, const double *A)
 {
    const int s = Width()*Height();
@@ -2175,6 +2288,143 @@ void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
    }
 }
 
+void CalcAdjugateRevDiff(const DenseMatrix &a, const DenseMatrix &adja_bar,
+                         DenseMatrix &a_bar)
+{
+#ifdef MFEM_DEBUG
+   if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
+   {
+      mfem_error("CalcAdjugateRevDiff(...)");
+   }
+   if (a.Width() != a_bar.Width() ||
+       a.Height() != a_bar.Height() ||
+       a_bar.Width() != adja_bar.Height() ||
+       a_bar.Height() != adja_bar.Width())
+   {
+      mfem_error("CalcAdjugateRefDiff(...)");
+   }
+#endif
+
+   if (a.Width() < a.Height())
+   {
+      const double *d = a.Data();
+      const double *ad_bar = adja_bar.Data();
+      double *d_bar = a_bar.Data();
+      if (a.Width() == 1)
+      {
+         // N x 1, N = 2,3
+         // ad[0] = d[0];
+         d_bar[0] = ad_bar[0];
+         // ad[1] = d[1];
+         d_bar[1] = ad_bar[1];
+         if (a.Height() == 3)
+         {
+            // ad[2] = d[2];
+            d_bar[2] = ad_bar[2];
+         }
+      }
+      else
+      {
+         // 3 x 2
+         // e, g, and f are needed during the reverse sweep
+         double e, g, f;
+         e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+         g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+         f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+
+         // start reverse sweep
+         a_bar = 0.0; // this zeros out d_bar[]
+         double e_bar = 0.0;
+         double g_bar = 0.0;
+         double f_bar = 0.0;
+         // ad[0] = d[0]*g - d[3]*f;
+         d_bar[0] += g*ad_bar[0];
+         d_bar[3] -= f*ad_bar[0];
+         g_bar += d[0]*ad_bar[0];
+         f_bar -= d[3]*ad_bar[0];
+         // ad[1] = d[3]*e - d[0]*f;
+         d_bar[3] += e*ad_bar[1];
+         d_bar[0] -= f*ad_bar[1];
+         e_bar += d[3]*ad_bar[1];
+         f_bar -= d[0]*ad_bar[1];
+         // ad[2] = d[1]*g - d[4]*f;
+         d_bar[1] += g*ad_bar[2];
+         d_bar[4] -= f*ad_bar[2];
+         g_bar += d[1]*ad_bar[2];
+         f_bar -= d[4]*ad_bar[2];
+         // ad[3] = d[4]*e - d[1]*f;
+         d_bar[4] += e*ad_bar[3];
+         d_bar[1] -= f*ad_bar[3];
+         e_bar += d[4]*ad_bar[3];
+         f_bar -= d[1]*ad_bar[3];
+         // ad[4] = d[2]*g - d[5]*f;
+         d_bar[2] += g*ad_bar[4];
+         d_bar[5] -= f*ad_bar[4];
+         g_bar += d[2]*ad_bar[4];
+         f_bar -= d[5]*ad_bar[4];
+         // ad[5] = d[5]*e - d[2]*f;
+         d_bar[5] += e*ad_bar[5];
+         d_bar[2] -= f*ad_bar[5];
+         e_bar += d[5]*ad_bar[5];
+         f_bar -= d[2]*ad_bar[5];
+
+         // e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+         d_bar[0] += 2.0*d[0]*e_bar;
+         d_bar[1] += 2.0*d[1]*e_bar;
+         d_bar[2] += 2.0*d[2]*e_bar;
+         // g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+         d_bar[3] += 2.0*d[3]*g_bar;
+         d_bar[4] += 2.0*d[4]*g_bar;
+         d_bar[5] += 2.0*d[5]*g_bar;
+         // f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+         d_bar[0] += d[3]*f_bar;
+         d_bar[3] += d[0]*f_bar;
+         d_bar[1] += d[4]*f_bar;
+         d_bar[4] += d[1]*f_bar;
+         d_bar[2] += d[5]*f_bar;
+         d_bar[5] += d[2]*f_bar;
+
+      }
+      return;
+   }
+
+   if (a.Width() == 1)
+   {
+      // adja(0,0) = 1.0;
+      a_bar(0,0) = 0.0;
+   }
+   else if (a.Width() == 2)
+   {
+      // adja(0,0) =  a(1,1);
+      a_bar(1,1) = adja_bar(0,0);
+      // adja(0,1) = -a(0,1);
+      a_bar(0,1) = -adja_bar(0,1);
+      // adja(1,0) = -a(1,0);
+      a_bar(1,0) = -adja_bar(1,0);
+      // adja(1,1) =  a(0,0);
+      a_bar(0,0) = adja_bar(1,1);
+   }
+   else
+   {
+      a_bar = 0.0;
+      for (int di1 = 0; di1 < 3; ++di1)
+      {
+         int it11 = (di1 + 1) % 3;
+         int it12 = (di1 + 2) % 3;
+         for (int di2 = 0; di2 < 3; ++di2)
+         {
+            int it21 = (di2 + 1) % 3;
+            int it22 = (di2 + 2) % 3;
+            // adja(di2,di1) = a(it11,it21)*a(it12,it22) - a(it11,it22)*a(it12,it21);
+            a_bar(it11,it21) += a(it12,it22)*adja_bar(di2,di1);
+            a_bar(it12,it22) += a(it11,it21)*adja_bar(di2,di1);
+            a_bar(it11,it22) -= a(it12,it21)*adja_bar(di2,di1);
+            a_bar(it12,it21) -= a(it11,it22)*adja_bar(di2,di1);
+         }
+      }
+   }
+}
+
 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
 {
    MFEM_ASSERT(a.Width() <= a.Height() && a.Width() >= 1 && a.Height() <= 3, "");
@@ -2209,7 +2459,9 @@ void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
             g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
             f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
             t = 1.0 / (e*g - f*f);
-            e *= t; g *= t; f *= t;
+            e *= t;
+            g *= t;
+            f *= t;
 
             id[0] = d[0]*g - d[3]*f;
             id[1] = d[3]*e - d[0]*f;
@@ -2281,6 +2533,194 @@ void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
    }
 }
 
+void CalcInverseRevDiff(const DenseMatrix &a, const DenseMatrix &inva_bar,
+                        DenseMatrix &a_bar)
+{
+#ifdef MFEM_DEBUG
+   if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
+   {
+      mfem_error("CalcInverseRevDiff(...)");
+   }
+   if (a.Width() != a_bar.Width() ||
+       a.Height() != a_bar.Height() ||
+       a_bar.Width() != inva_bar.Height() ||
+       a_bar.Height() != inva_bar.Width())
+   {
+      mfem_error("CalcInverseRevDiff(...)");
+   }
+#endif
+
+   if (a.Width() < a.Height())
+   {
+      const double *d = a.Data();
+      const double *id_bar = inva_bar.Data();
+      double *d_bar = a_bar.Data();
+      if (a.Height() == 2)
+      {
+         double t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
+
+         /// id[0] = d[0] * t;
+         d_bar[0] += id_bar[0] * t;
+         double t_bar = id_bar[0] * d[0];
+
+         /// id[1] = d[1] * t;
+         d_bar[1] += id_bar[1] * t;
+         t_bar += id_bar[1] * d[1];
+
+         /// t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
+         d_bar[0] -= t_bar * 2 * d[0] / pow(d[0]*d[0] + d[1]*d[1], 2);
+         d_bar[1] -= t_bar * 2 * d[1] / pow(d[0]*d[0] + d[1]*d[1], 2);
+      }
+      else
+      {
+         if (a.Width() == 1)
+         {
+            double t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
+
+            /// id[0] = d[0] * t;
+            d_bar[0] += id_bar[0] * t;
+            double t_bar = id_bar[0] * d[0];
+
+            /// id[1] = d[1] * t;
+            d_bar[1] += id_bar[1] * t;
+            t_bar += id_bar[1] * d[1];
+
+            /// id[2] = d[2] * t;
+            d_bar[2] += id_bar[2] * t;
+            t_bar += id_bar[2] * d[2];
+
+            /// t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
+            d_bar[0] -= t_bar * 2 * d[0] / pow(d[0]*d[0] + d[1]*d[1] + d[2]*d[2], 2);
+            d_bar[1] -= t_bar * 2 * d[1] / pow(d[0]*d[0] + d[1]*d[1] + d[2]*d[2], 2);
+            d_bar[2] -= t_bar * 2 * d[2] / pow(d[0]*d[0] + d[1]*d[1] + d[2]*d[2], 2);
+         }
+         else
+         {
+            double e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+            double g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+            double f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+            double t = 1.0 / (e*g - f*f);
+
+            double ee = e * t;
+            double gg = g * t;
+            double ff = f * t;
+
+            /// id[0] = d[0]*g - d[3]*f;
+            d_bar[0] += id_bar[0]*gg;
+            double gg_bar = id_bar[0] * d[0];
+            d_bar[3] += -id_bar[0] * ff;
+            double ff_bar = -id_bar[0] * d[3];
+
+            /// id[1] = d[3]*e - d[0]*f;
+            d_bar[3] += id_bar[1] * ee;
+            double ee_bar = id_bar[1] * d[3];
+            d_bar[0] += -id_bar[1] * ff;
+            ff_bar -= id_bar[1] * d[0];
+
+            /// id[2] = d[1]*g - d[4]*f;
+            d_bar[1] += id_bar[2] * gg;
+            gg_bar += id_bar[2] * d[1];
+            d_bar[4] += -id_bar[2] * ff;
+            ff_bar -= id_bar[2] * d[4];
+
+            /// id[3] = d[4]*e - d[1]*f;
+            d_bar[4] += id_bar[3] * ee;
+            ee_bar += id_bar[3] * d[4];
+            d_bar[1] += -id_bar[3] * ff;
+            ff_bar -= id_bar[3] * d[1];
+
+            /// id[4] = d[2]*g - d[5]*f;
+            d_bar[2] += id_bar[4] * gg;
+            gg_bar += id_bar[4] * d[2];
+            d_bar[5] += -id_bar[4] * ff;
+            ff_bar -= id_bar[4] * d[5];
+
+            /// id[5] = d[5]*e - d[2]*f;
+            d_bar[5] += id_bar[5] * ee;
+            ee_bar += id_bar[5] * d[5];
+            d_bar[2] += -id_bar[5] * ff;
+            ff_bar -= id_bar[5] * d[2];
+
+
+            /// double ff = f * t;
+            double t_bar = ff_bar * f;
+            double f_bar = ff_bar * t;
+
+            /// double gg = g * t;
+            t_bar += gg_bar * g;
+            double g_bar = gg_bar * t;
+
+            /// double ee = e * t;
+            t_bar += ee_bar * e;
+            double e_bar = ee_bar * t;
+
+            // /// f *= t;
+            // double t_bar = f_bar * f / t;
+            // f_bar *= t;
+
+            // /// g *= t;
+            // t_bar += g_bar * g / t;
+            // g_bar *= t;
+
+            // /// e *= t;
+            // t_bar += e_bar * e / t;
+            // e_bar *= t;
+
+            /// double t = 1.0 / (e*g - f*f);
+            e_bar -= t_bar * g / pow(e*g - f*f, 2);
+            g_bar -= t_bar * e / pow(e*g - f*f, 2);
+            f_bar += t_bar * 2*f / pow(e*g - f*f, 2);
+
+            /// double f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+            d_bar[0] += f_bar * d[3];
+            d_bar[3] += f_bar * d[0];
+            d_bar[1] += f_bar * d[4];
+            d_bar[4] += f_bar * d[1];
+            d_bar[2] += f_bar * d[5];
+            d_bar[5] += f_bar * d[2];
+
+            /// double g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+            d_bar[3] += g_bar * 2 * d[3];
+            d_bar[4] += g_bar * 2 * d[4];
+            d_bar[5] += g_bar * 2 * d[5];
+
+            /// double e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+            d_bar[0] += e_bar * 2 * d[0];
+            d_bar[1] += e_bar * 2 * d[1];
+            d_bar[2] += e_bar * 2 * d[2];
+         }
+      }
+      return;
+   }
+
+#ifdef MFEM_DEBUG
+   double t = a.Det();
+   MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(a.FNorm()/a.Width(), a.Width()),
+               "singular matrix!");
+#endif
+
+   double inva_buffer[9] = {};
+   DenseMatrix inva(inva_buffer, a.Height(), a.Width());
+
+   switch (a.Height())
+   {
+      case 1:
+         inva(0,0) = 1.0 / a.Det();
+         break;
+      case 2:
+         kernels::CalcInverse<2>(a.Data(), inva.Data());
+         break;
+      case 3:
+         kernels::CalcInverse<3>(a.Data(), inva.Data());
+         break;
+   }
+
+   double tmp_buffer[9] = {};
+   DenseMatrix tmp(tmp_buffer, a.Height(), a.Width());
+   MultAtB(inva, inva_bar, tmp);
+   AddMult_a_ABt(-1, tmp, inva, a_bar);
+}
+
 void CalcOrtho(const DenseMatrix &J, Vector &n)
 {
    MFEM_ASSERT( ((J.Height() == 2 && J.Width() == 1)
@@ -2307,6 +2747,50 @@ void CalcOrtho(const DenseMatrix &J, Vector &n)
    }
 }
 
+void CalcOrthoRevDiff(const DenseMatrix &J, const Vector &n_bar,
+                      DenseMatrix &J_bar)
+{
+   MFEM_ASSERT(((J.Height() == 2 && J.Width() == 1) ||
+                (J.Height() == 3 && J.Width() == 2)) &&
+               (J.Height() == n_bar.Size()),
+               "Matrix must be 3x2 or 2x1, "
+               << "and the Vector must be sized with the rows. "
+               << " J.Height() = " << J.Height()
+               << ", J.Width() = " << J.Width()
+               << ", n_bar.Size() = " << n_bar.Size());
+   MFEM_ASSERT((J.Height() == J_bar.Height() && J.Width() == J_bar.Width()),
+               "Input matrix and derivative matrix must be the same size.");
+
+   const double *d = J.Data();
+   double *d_bar = J_bar.Data();
+   if (J.Height() == 2)
+   {
+      // n(0) =  d[1];
+      d_bar[1] = n_bar(0);
+      // n(1) = -d[0];
+      d_bar[0] = -n_bar(1);
+   }
+   else
+   {
+      J_bar = 0.0;
+      // n(0) = d[1]*d[5] - d[2]*d[4];
+      d_bar[1] += d[5]*n_bar(0);
+      d_bar[5] += d[1]*n_bar(0);
+      d_bar[2] -= d[4]*n_bar(0);
+      d_bar[4] -= d[2]*n_bar(0);
+      // n(1) = d[2]*d[3] - d[0]*d[5];
+      d_bar[2] += d[3]*n_bar(1);
+      d_bar[3] += d[2]*n_bar(1);
+      d_bar[0] -= d[5]*n_bar(1);
+      d_bar[5] -= d[0]*n_bar(1);
+      // n(2) = d[0]*d[4] - d[1]*d[3];
+      d_bar[0] += d[4]*n_bar(2);
+      d_bar[4] += d[0]*n_bar(2);
+      d_bar[1] -= d[3]*n_bar(2);
+      d_bar[3] -= d[1]*n_bar(2);
+   }
+}
+
 void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
 {
    const int height = a.Height();

linalg/densemat.hpp

@@ -197,6 +197,12 @@ public:
 
    double Weight() const;
 
+   /// Evaluate the derivative of Det() w.r.t. the matrix entries
+   void DetRevDiff(DenseMatrix &A_bar) const;
+
+   /// Evaluate the derivative of Weight() w.r.t. the matrix entries
+   void WeightRevDiff(DenseMatrix &A_bar) const;
+
    /** @brief Set the matrix to alpha * A, assuming that A has the same
        dimensions as the matrix and uses column-major layout. */
    void Set(double alpha, const double *A);
@@ -457,6 +463,10 @@ void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja);
 /// Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
 void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat);
 
+/// Reverse-mode sensitivities of adj(A) with respect to the entries in A
+void CalcAdjugateRevDiff(const DenseMatrix &a, const DenseMatrix &adja_bar,
+                         DenseMatrix &a_bar);
+
 /** Calculate the inverse of a matrix (for NxN matrices, N=1,2,3) or the
     left inverse (A^t.A)^{-1}.A^t (for 2x1, 3x1, or 3x2 matrices) */
 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
@@ -464,11 +474,19 @@ void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
 /// Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
 
+/// Reverse-mode sensitivities of inv(A) with respect to the entries in A
+void CalcInverseRevDiff(const DenseMatrix &a, const DenseMatrix &inva_bar,
+                        DenseMatrix &a_bar);
+
 /** For a given Nx(N-1) (N=2,3) matrix J, compute a vector n such that
     n_k = (-1)^{k+1} det(J_k), k=1,..,N, where J_k is the matrix J with the
     k-th row removed. Note: J^t.n = 0, det([n|J])=|n|^2=det(J^t.J). */
 void CalcOrtho(const DenseMatrix &J, Vector &n);
 
+/// The reverse-mode differentiation of CalcOrtho()
+void CalcOrthoRevDiff(const DenseMatrix &J, const Vector &n_bar,
+                      DenseMatrix &J_bar);
+
 /// Calculate the matrix A.At
 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
 

tests/unit/CMakeLists.txt

@@ -56,17 +56,20 @@ set(UNIT_TESTS_SRCS
   fem/test_blocknonlinearform.cpp
   fem/test_calcshape.cpp
   fem/test_coefficient.cpp
+  fem/test_coeff_revdiff.cpp
   fem/test_datacollection.cpp
   fem/test_derefine.cpp
   fem/test_estimator.cpp
   fem/test_face_elem_trans.cpp
   fem/test_face_permutation.cpp
   fem/test_fe.cpp
+  fem/test_fe_revdiff.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp
+  fem/test_eltrans_revdiff.cpp
   fem/test_lexicographic_ordering.cpp
   fem/test_lin_interp.cpp
   fem/test_linear_fes.cpp

tests/unit/fem/test_coeff_revdiff.cpp

@@ -0,0 +1,325 @@
+// Copyright (c) 2010-2021, 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 "catch.hpp"
+
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <random>
+
+using namespace mfem;
+
+namespace
+{
+// String used to define a single element mesh, with a c-shaped quad element
+std::string mesh_str =
+   "MFEM mesh v1.0"                    "\n\n"
+   "dimension"                           "\n"
+   "2"                                 "\n\n"
+   "elements"                            "\n"
+   "1"                                   "\n"
+   "1 3 0 1 2 3"                       "\n\n"
+   "boundary"                            "\n"
+   "0"                                 "\n\n"
+   "vertices"                            "\n"
+   "4"                                 "\n\n"
+   "nodes"                               "\n"
+   "FiniteElementSpace"                  "\n"
+   "FiniteElementCollection: Quadratic"  "\n"
+   "VDim: 2"                             "\n"
+   "Ordering: 1"                         "\n"
+   "0 0"                                 "\n"
+   "0 2"                                 "\n"
+   "0 6"                                 "\n"
+   "0 8"                                 "\n"
+   "0 1"                                 "\n"
+   "-6 4"                                "\n"
+   "0 7"                                 "\n"
+   "-8 4"                                "\n"
+   "-7 4"                                "\n";
+
+double scalar_func(const Vector &x)
+{
+   double q = 0;
+   for (int i = 0; i < x.Size(); ++i)
+   {
+      q += pow(x(i), 2);
+   }
+   return q;
+}
+
+void scalar_funcRevDiff(const mfem::Vector &x, const double q_bar,
+                        mfem::Vector &x_bar)
+{
+   for (int i = 0; i < x.Size(); ++i)
+   {
+      x_bar(i) += q_bar * 2 * x(i);
+   }
+}
+
+double scalar_func2(const Vector &x)
+{
+   double q = 0;
+   for (int i = 0; i < x.Size(); ++i)
+   {
+      q += x(i);
+   }
+   return q;
+}
+
+void scalar_func2RevDiff(const mfem::Vector &x, const double q_bar,
+                         mfem::Vector &x_bar)
+{
+   for (int i = 0; i < x.Size(); ++i)
+   {
+      x_bar(i) += q_bar;
+   }
+}
+
+void func2D(const Vector &x, Vector &y)
+{
+   y.SetSize(2);
+   y(0) = x(0)*x(0) - x(1);
+   y(1) = x(0) * exp(x(1));
+}
+
+void func2DRevDiff(const Vector &x, const Vector &v_bar, Vector &x_bar)
+{
+   x_bar(0) = v_bar(0) * 2*x(0) + v_bar(1) * exp(x(1));
+   x_bar(1) = -v_bar(0) + v_bar(1) * x(0) * exp(x(1));
+}
+
+void func3D(const Vector &x, Vector &y)
+{
+   y.SetSize(3);
+   y(0) = x(0)*x(0) - x(1);
+   y(1) = x(0) * exp(x(1));
+   y(2) = x(2)*x(0) - x(1);
+}
+
+void func3DRevDiff(const Vector &x, const Vector &v_bar, Vector &x_bar)
+{
+   x_bar(0) = v_bar(0) * 2*x(0) + v_bar(1) * exp(x(1)) + v_bar(2)*x(2);
+   x_bar(1) = -v_bar(0) + v_bar(1) * x(0) * exp(x(1)) - v_bar(2);
+   x_bar(2) = v_bar(2) * x(0);
+}
+
+void runScalarTest(Mesh &mesh, Coefficient &q)
+{
+   constexpr double eps_fd = 1e-5;
+   std::default_random_engine generator;
+   std::uniform_real_distribution<double> distribution(-1.0,1.0);
+
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      H1_FECollection fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      DenseMatrix &coords = trans.GetPointMat();
+      DenseMatrix coords_bar(coords.Height(), coords.Width());
+
+      double Q_bar = distribution(generator);
+
+      int order = trans.OrderW() + 2 * el.GetOrder();
+      const IntegrationRule *ir = &IntRules.Get(el.GetGeomType(), order);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint (&ip);
+
+         // reverse-mode differentiation of Eval
+         coords_bar = 0.0;
+         q.EvalRevDiff(Q_bar, trans, ip, coords_bar);
+
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               double Q_fd = q.Eval(trans, ip);
+               coords(di, n) -= 2.0*eps_fd;
+               Q_fd -= q.Eval(trans, ip);
+               Q_fd /= (2.0 * eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = Q_bar * Q_fd;
+
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+}
+
+void runVectorTest(Mesh &mesh, VectorCoefficient &vc)
+{
+   constexpr double eps_fd = 1e-5;
+   std::default_random_engine generator;
+   std::uniform_real_distribution<double> distribution(-1.0,1.0);
+
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      ND_FECollection fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      DenseMatrix &coords = trans.GetPointMat();
+      DenseMatrix coords_bar(coords.Height(), coords.Width());
+
+      Vector V_bar(dim), V_fd(dim), V_pert(dim);
+      for (int i = 0; i < V_bar.Size(); ++i)
+      {
+         V_bar(i) = distribution(generator);
+      }
+
+      int order = trans.OrderW() + 2 * el.GetOrder();
+      const IntegrationRule *ir = &IntRules.Get(el.GetGeomType(), order);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint (&ip);
+
+         // reverse-mode differentiation of Eval
+         coords_bar = 0.0;
+         vc.EvalRevDiff(V_bar, trans, ip, coords_bar);
+
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               vc.Eval(V_fd, trans, ip);
+               coords(di, n) -= 2.0*eps_fd;
+               vc.Eval(V_pert, trans, ip);
+               V_fd -= V_pert;
+               V_fd /= (2.0 * eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = V_bar * V_fd;
+
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+}
+
+} // anonymous namespace
+
+namespace coeff_revdiff
+{
+
+TEST_CASE("CoeffRevDiff::FunctionCoefficient::EvalRevDiff_2D")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+
+   REQUIRE( mesh2D.GetNE() == 1 );
+   REQUIRE( mesh2D.GetNodes() != NULL );
+
+   FunctionCoefficient c1(scalar_func, scalar_funcRevDiff);
+   runScalarTest(mesh2D, c1);
+
+   FunctionCoefficient c2(scalar_func2, scalar_func2RevDiff);
+   runScalarTest(mesh2D, c2);
+}
+
+TEST_CASE("CoeffRevDiff::ProductCoefficient::EvalRevDiff_2D 1 Coeff")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+
+   REQUIRE( mesh2D.GetNE() == 1 );
+   REQUIRE( mesh2D.GetNodes() != NULL );
+
+   FunctionCoefficient c1(scalar_func, scalar_funcRevDiff);
+   ProductCoefficient prod(2.0, c1);
+
+   runScalarTest(mesh2D, prod);
+}
+
+TEST_CASE("CoeffRevDiff::ProductCoefficient::EvalRevDiff_2D 2 Coeffs")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+
+   REQUIRE( mesh2D.GetNE() == 1 );
+   REQUIRE( mesh2D.GetNodes() != NULL );
+
+   FunctionCoefficient c1(scalar_func, scalar_funcRevDiff);
+   FunctionCoefficient c2(scalar_func2, scalar_func2RevDiff);
+
+   ProductCoefficient prod(c1, c2);
+
+   runScalarTest(mesh2D, prod);
+}
+
+
+TEST_CASE("CoeffRevDiff::VectorFunctionCoefficient::EvalRevDiff_2D")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+
+   REQUIRE( mesh2D.GetNE() == 1 );
+   REQUIRE( mesh2D.GetNodes() != NULL );
+
+   VectorFunctionCoefficient vc2D(2, func2D, func2DRevDiff);
+
+   runVectorTest(mesh2D, vc2D);
+}
+
+TEST_CASE("CoeffRevDiff::VectorFunctionCoefficient::EvalRevDiff_3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2,
+                                       Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   mesh3D.EnsureNodes();
+
+   VectorFunctionCoefficient vc3D(3, func3D, func3DRevDiff);
+
+   runVectorTest(mesh3D, vc3D);
+}
+
+TEST_CASE("CoeffRevDiff::ScalarVectorProductCoefficient::EvalRevDiff_3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2,
+                                       Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   mesh3D.EnsureNodes();
+
+   VectorFunctionCoefficient vfc(3, func3D, func3DRevDiff);
+
+   ScalarVectorProductCoefficient vc(2.0, vfc);
+
+   runVectorTest(mesh3D, vc);
+}
+
+} // namespace coeff_revdiff

tests/unit/fem/test_eltrans_revdiff.cpp

@@ -0,0 +1,268 @@
+// Copyright (c) 2010-2021, 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 "catch.hpp"
+
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <fstream>
+
+using namespace mfem;
+
+// String used to define a single element mesh, with a c-shaped quad element
+std::string mesh_str =
+   "MFEM mesh v1.0"                    "\n\n"
+   "dimension"                           "\n"
+   "2"                                 "\n\n"
+   "elements"                            "\n"
+   "1"                                   "\n"
+   "1 3 0 1 2 3"                       "\n\n"
+   "boundary"                            "\n"
+   "0"                                 "\n\n"
+   "vertices"                            "\n"
+   "4"                                 "\n\n"
+   "nodes"                               "\n"
+   "FiniteElementSpace"                  "\n"
+   "FiniteElementCollection: Quadratic"  "\n"
+   "VDim: 2"                             "\n"
+   "Ordering: 1"                         "\n"
+   "0 0"                                 "\n"
+   "0 2"                                 "\n"
+   "0 6"                                 "\n"
+   "0 8"                                 "\n"
+   "0 1"                                 "\n"
+   "-6 4"                                "\n"
+   "0 7"                                 "\n"
+   "-8 4"                                "\n"
+   "-7 4"                                "\n";
+
+TEST_CASE("IsoparametricTransformation reverse-mode differentiation",
+          "[IsoparametricTransformation]")
+{
+   constexpr double eps_fd = 1e-5;
+
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh(meshStr);
+
+   REQUIRE( mesh.GetNE() == 1 );
+   REQUIRE( mesh.GetNodes() != NULL );
+
+   bool dumpMesh = false;
+   if (dumpMesh)
+   {
+      std::ofstream mesh_ostream("isoparametric-revdiff-mesh.vtk");
+      mesh_ostream.precision(14);
+      int refine = 10;
+      mesh.PrintVTK(mesh_ostream, refine);
+   }
+
+   // Create the transformation and get integration rule
+   IsoparametricTransformation trans;
+   mesh.GetElementTransformation(0, &trans);
+   const int intorder = 5;
+   const IntegrationRule *ir = &IntRules.Get(mesh.GetElementBaseGeometry(0),
+                                             intorder);
+   DenseMatrix &coords = trans.GetPointMat();
+   DenseMatrix coords_bar(coords.Height(), coords.Width());
+
+   SECTION("TransformRevDiff")
+   {
+      // x_bar(i) is the weight on the (i)th entry of the coordinate x;
+      // the values are not important for this test.
+      double x_bar_data[4] = {2.5, -3.2};
+      Vector x_bar(x_bar_data, 2);
+      Vector x_fd(2), x_pert(2);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         // reverse-mode differentiation of coordinate transformation
+         coords_bar = 0.0;
+         trans.TransformRevDiff(ip, x_bar, coords_bar);
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.Transform(ip, x_fd);
+               coords(di, n) -= 2.0*eps_fd;
+               trans.Transform(ip, x_pert);
+               x_fd -= x_pert;
+               x_fd *= 1.0/(2.0*eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = 0.0;
+               for (int j = 0; j < x_bar.Size(); ++j)
+               {
+                  x_bar_fd += x_bar(j)*x_fd(j);
+               }
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+
+   SECTION("JacobianRevDiff")
+   {
+      // dFdx_bar(i,j) is the weight on the (i,j)th entry of the Jacobian;
+      // the values are not important for this test.
+      double dFdx_bar_data[4] = {2.0, -3.0, 4.0, -1.0};
+      DenseMatrix dFdx_bar(dFdx_bar_data, 2, 2);
+      DenseMatrix dFdx_fd(2,2);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         // reverse-mode differentiation of Jacobian of mapping
+         coords_bar = 0.0;
+         trans.JacobianRevDiff(dFdx_bar, coords_bar);
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Jacobian
+               dFdx_fd = trans.Jacobian();
+               coords(di, n) -= 2.0*eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Jacobian
+               dFdx_fd -= trans.Jacobian();
+               dFdx_fd *= 1.0/(2.0*eps_fd);
+               coords(di, n) += eps_fd;
+               double dFdx_bar_fd = 0.0;
+               for (int j = 0; j < dFdx_bar.Height(); ++j)
+               {
+                  for (int k = 0; k < dFdx_bar.Width(); ++k)
+                  {
+                     dFdx_bar_fd += dFdx_bar(j,k)*dFdx_fd(j,k);
+                  }
+               }
+               REQUIRE(coords_bar(di, n) == Approx(dFdx_bar_fd));
+            }
+         }
+      }
+   }
+
+   SECTION("AdjugateJacobianRevDiff")
+   {
+      // adjJ_bar(i,j) is the weight on the (i,j)th entry of the Adjugate;
+      // the values are not important for this test.
+      double adjJ_bar_data[4] = {2.0, -3.0, 4.0, -1.0};
+      DenseMatrix adjJ_bar(adjJ_bar_data, 2, 2);
+      DenseMatrix adjJ_fd(2,2);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         // reverse-mode differentiation of Adjugate of mapping
+         coords_bar = 0.0;
+         trans.AdjugateJacobianRevDiff(adjJ_bar, coords_bar);
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Adjugate
+               adjJ_fd = trans.AdjugateJacobian();
+               coords(di, n) -= 2.0*eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Adjugate
+               adjJ_fd -= trans.AdjugateJacobian();
+               adjJ_fd *= 1.0/(2.0*eps_fd);
+               coords(di, n) += eps_fd;
+               double adjJ_bar_fd = 0.0;
+               for (int j = 0; j < adjJ_bar.Height(); ++j)
+               {
+                  for (int k = 0; k < adjJ_bar.Width(); ++k)
+                  {
+                     adjJ_bar_fd += adjJ_bar(j,k)*adjJ_fd(j,k);
+                  }
+               }
+               REQUIRE(coords_bar(di, n) == Approx(adjJ_bar_fd));
+            }
+         }
+      }
+   }
+
+   SECTION("InverseJacobianRevDiff")
+   {
+      // invJ_bar(i,j) is the weight on the (i,j)th entry of the Inverse;
+      // the values are not important for this test.
+      double invJ_bar_data[4] = {2.0, -3.0, 4.0, -1.0};
+      DenseMatrix invJ_bar(invJ_bar_data, 2, 2);
+      DenseMatrix invJ_fd(2,2);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         // reverse-mode differentiation of Inverse of mapping
+         coords_bar = 0.0;
+         trans.InverseJacobianRevDiff(invJ_bar, coords_bar);
+         // get the weighted derivatives using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Inverse
+               invJ_fd = trans.InverseJacobian();
+               coords(di, n) -= 2.0*eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Inverse
+               invJ_fd -= trans.InverseJacobian();
+               invJ_fd *= 1.0/(2.0*eps_fd);
+               coords(di, n) += eps_fd;
+               double invJ_bar_fd = 0.0;
+               for (int j = 0; j < invJ_bar.Height(); ++j)
+               {
+                  for (int k = 0; k < invJ_bar.Width(); ++k)
+                  {
+                     invJ_bar_fd += invJ_bar(j,k)*invJ_fd(j,k);
+                  }
+               }
+               REQUIRE(coords_bar(di, n) == Approx(invJ_bar_fd));
+            }
+         }
+      }
+   }
+
+   SECTION("WeightRevDiff")
+   {
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         // get the gradient of the Weight() using reverse mode
+         coords_bar = 0.0;
+         trans.WeightRevDiff(coords_bar);
+         // get the gradient of the Weight() using finite difference method
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Weight
+               double dWeight_fd = trans.Weight();
+               coords(di, n) -= 2.0*eps_fd;
+               trans.SetIntPoint(&ip); // force re-evaluation of Weight
+               dWeight_fd -= trans.Weight();
+               dWeight_fd /= (2.0*eps_fd);
+               coords(di, n) += eps_fd;
+               REQUIRE(coords_bar(di, n) == Approx(dWeight_fd).margin(1e-10));
+            }
+         }
+      }
+   }
+
+}

tests/unit/fem/test_fe_revdiff.cpp

@@ -0,0 +1,419 @@
+// Copyright (c) 2010-2021, 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 "catch.hpp"
+
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <random>
+
+using namespace mfem;
+
+namespace
+{
+// String used to define a single element mesh, with a c-shaped quad element
+std::string mesh_str =
+   "MFEM mesh v1.0"                    "\n\n"
+   "dimension"                           "\n"
+   "2"                                 "\n\n"
+   "elements"                            "\n"
+   "1"                                   "\n"
+   "1 3 0 1 2 3"                       "\n\n"
+   "boundary"                            "\n"
+   "0"                                 "\n\n"
+   "vertices"                            "\n"
+   "4"                                 "\n\n"
+   "nodes"                               "\n"
+   "FiniteElementSpace"                  "\n"
+   "FiniteElementCollection: Quadratic"  "\n"
+   "VDim: 2"                             "\n"
+   "Ordering: 1"                         "\n"
+   "0 0"                                 "\n"
+   "0 2"                                 "\n"
+   "0 6"                                 "\n"
+   "0 8"                                 "\n"
+   "0 1"                                 "\n"
+   "-6 4"                                "\n"
+   "0 7"                                 "\n"
+   "-8 4"                                "\n"
+   "-7 4"                                "\n";
+
+void func2D(const Vector &x, Vector &y)
+{
+   y.SetSize(2);
+   y(0) = x(0)*x(0) - x(1);
+   y(1) = x(0) * exp(x(1));
+}
+
+void func2DRevDiff(const Vector &x, const Vector &v_bar, Vector &x_bar)
+{
+   x_bar(0) = v_bar(0) * 2*x(0) + v_bar(1) * exp(x(1));
+   x_bar(1) = -v_bar(0) + v_bar(1) * x(0) * exp(x(1));
+}
+
+void func3D(const Vector &x, Vector &y)
+{
+   y.SetSize(3);
+   y(0) = x(0)*x(0) - x(1);
+   y(1) = x(0) * exp(x(1));
+   y(2) = x(2)*x(0) - x(1);
+}
+
+void func3DRevDiff(const Vector &x, const Vector &v_bar, Vector &x_bar)
+{
+   x_bar(0) = v_bar(0) * 2*x(0) + v_bar(1) * exp(x(1)) + v_bar(2)*x(2);
+   x_bar(1) = -v_bar(0) + v_bar(1) * x(0) * exp(x(1)) - v_bar(2);
+   x_bar(2) = v_bar(2) * x(0);
+}
+
+} // anonymous namespace
+
+namespace fe_revdiff
+{
+
+template<typename T>
+void runProjectRevDiffTest(Mesh &mesh, VectorCoefficient &vc);
+
+template<typename T>
+void runCalcPhysShapeRevDiffTest(Mesh &mesh);
+
+template<typename T>
+void runCalcVShapeRevDiffTest(Mesh &mesh);
+
+void runCalcPhysCurlShapeRevDiffTest(Mesh &mesh);
+
+constexpr double eps_fd = 1e-5;
+std::default_random_engine generator;
+std::uniform_real_distribution<double> distribution(-1.0,1.0);
+
+TEST_CASE("VectorFiniteElement::ProjectRevDiff - 2D")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+
+   REQUIRE(mesh2D.GetNE() == 1);
+   REQUIRE(mesh2D.GetNodes() != nullptr);
+
+   VectorFunctionCoefficient vc2D(2, func2D, func2DRevDiff);
+
+   runProjectRevDiffTest<RT_FECollection>(mesh2D, vc2D);
+   runProjectRevDiffTest<ND_FECollection>(mesh2D, vc2D);
+}
+
+TEST_CASE("VectorFiniteElement::ProjectRevDiff - 3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2, Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   // mesh3D.ReorientTetMesh();
+   mesh3D.EnsureNodes();
+
+   VectorFunctionCoefficient vc3D(3, func3D, func3DRevDiff);
+
+   runProjectRevDiffTest<RT_FECollection>(mesh3D, vc3D);
+   runProjectRevDiffTest<ND_FECollection>(mesh3D, vc3D);
+}
+
+TEST_CASE("FiniteElement::CalcPhysShapeRevDiff")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+   REQUIRE(mesh2D.GetNE() == 1);
+   REQUIRE(mesh2D.GetNodes() != nullptr);
+   runCalcPhysShapeRevDiffTest<H1_FECollection>(mesh2D);
+   runCalcPhysShapeRevDiffTest<L2_FECollection>(mesh2D);
+}
+
+TEST_CASE("VectorFiniteElement::CalcVShape_RTRevDiff - 2D")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+   REQUIRE(mesh2D.GetNE() == 1);
+   REQUIRE(mesh2D.GetNodes() != nullptr);
+   runCalcVShapeRevDiffTest<RT_FECollection>(mesh2D);
+}
+
+TEST_CASE("VectorFiniteElement::CalcVShape_NDRevDiff - 2D")
+{
+   // Create quadratic mesh with single C-shaped quadrilateral
+   std::stringstream meshStr;
+   meshStr << mesh_str;
+   Mesh mesh2D(meshStr);
+   REQUIRE(mesh2D.GetNE() == 1);
+   REQUIRE(mesh2D.GetNodes() != nullptr);
+   runCalcVShapeRevDiffTest<ND_FECollection>(mesh2D);
+}
+
+TEST_CASE("VectorFiniteElement::CalcVShape_RTRevDiff - 3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2, Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   // mesh3D.ReorientTetMesh();
+   mesh3D.EnsureNodes();
+   runCalcVShapeRevDiffTest<RT_FECollection>(mesh3D);
+}
+
+TEST_CASE("VectorFiniteElement::CalcVShape_NDRevDiff - 3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2, Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   // mesh3D.ReorientTetMesh();
+   mesh3D.EnsureNodes();
+   runCalcVShapeRevDiffTest<ND_FECollection>(mesh3D);
+}
+
+TEST_CASE("FiniteElement::CalcPhysCurlShapeRevDiff - 3D")
+{
+   auto mesh3D = Mesh::MakeCartesian3D(2, 2, 2, Element::TETRAHEDRON,
+                                       2.0, 1.0, 3.0, true);
+   // mesh3D.ReorientTetMesh();
+   mesh3D.EnsureNodes();
+   runCalcPhysCurlShapeRevDiffTest(mesh3D);
+}
+
+template<typename T>
+void runProjectRevDiffTest(Mesh &mesh, VectorCoefficient &vc)
+{
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      T fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      // P_bar is the vector contracted with the derivative of the projection
+      // the values are not important for this test
+      const int dof = el.GetDof();
+      Vector P_bar(dof);
+      for (int i = 0; i < P_bar.Size(); ++i)
+      {
+         P_bar(i) = distribution(generator);
+      }
+
+      // reverse-mode differentiation of projection
+      DenseMatrix &coords = trans.GetPointMat();
+      DenseMatrix coords_bar(coords.Height(), coords.Width());
+      coords_bar = 0.0;
+      el.ProjectRevDiff(P_bar, vc, trans, coords_bar);
+
+      // get the weighted derivatives using finite difference method
+      Vector dofs_fd(dof), dofs_pert(dof);
+      for (int n = 0; n < coords.Width(); ++n)
+      {
+         for (int di = 0; di < coords.Height(); ++di)
+         {
+            coords(di, n) += eps_fd;
+            trans.Reset();
+            el.Project(vc, trans, dofs_fd);
+            coords(di, n) -= 2.0 * eps_fd;
+            trans.Reset();
+            el.Project(vc, trans, dofs_pert);
+            dofs_fd -= dofs_pert;
+            dofs_fd *= 1.0 / (2.0 * eps_fd);
+            coords(di, n) += eps_fd;
+            double x_bar_fd = P_bar * dofs_fd;
+
+            REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+         }
+      }
+   }
+}
+
+template<typename T>
+void runCalcPhysShapeRevDiffTest(Mesh &mesh)
+{
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      T fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      int order = trans.OrderW() + 2 * el.GetOrder();
+      const IntegrationRule *ir = &IntRules.Get(el.GetGeomType(), order);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+
+         const int dof = el.GetDof();
+         Vector shape_bar(dof);
+         for (int k = 0; k < shape_bar.Size(); ++k)
+         {
+            shape_bar(k) = distribution(generator);
+         }
+
+         // reverse-mode differentiation CalcVShape
+         DenseMatrix &coords = trans.GetPointMat();
+         DenseMatrix coords_bar(coords.Height(), coords.Width());
+         coords_bar = 0.0;
+         el.CalcPhysShapeRevDiff(trans, shape_bar, coords_bar);
+
+         // get the weighted derivatives using finite difference method
+         Vector shape_fd(dof), shape_pert(dof);
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.Reset();
+               el.CalcPhysShape(trans, shape_fd);
+               coords(di, n) -= 2.0 * eps_fd;
+               trans.Reset();
+               el.CalcPhysShape(trans, shape_pert);
+               shape_fd -= shape_pert;
+               shape_fd *= 1.0 / (2.0 * eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = shape_bar * shape_fd;
+
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+}
+
+template<typename T>
+void runCalcVShapeRevDiffTest(Mesh &mesh)
+{
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      T fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      int order = trans.OrderW() + 2 * el.GetOrder();
+      const IntegrationRule *ir = &IntRules.Get(el.GetGeomType(), order);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+
+         const int dof = el.GetDof();
+         const int el_dim = el.GetDim();
+         DenseMatrix vshape_bar(dof, el_dim);
+         for (int k = 0; k < vshape_bar.Width(); ++k)
+         {
+            for (int j = 0; j < vshape_bar.Height(); ++j)
+            {
+               vshape_bar(j, k) = distribution(generator);
+            }
+         }
+
+         // reverse-mode differentiation CalcVShape
+         DenseMatrix &coords = trans.GetPointMat();
+         DenseMatrix coords_bar(coords.Height(), coords.Width());
+         coords_bar = 0.0;
+         el.CalcVShapeRevDiff(trans, vshape_bar, coords_bar);
+
+         // get the weighted derivatives using finite difference method
+         DenseMatrix vshape_fd(dof, el_dim), vshape_pert(dof, el_dim);
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.Reset();
+               el.CalcVShape(trans, vshape_fd);
+               coords(di, n) -= 2.0 * eps_fd;
+               trans.Reset();
+               el.CalcVShape(trans, vshape_pert);
+               vshape_fd -= vshape_pert;
+               vshape_fd *= 1.0 / (2.0 * eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = vshape_bar * vshape_fd;
+
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+}
+
+void runCalcPhysCurlShapeRevDiffTest(Mesh &mesh)
+{
+   for (int p = 1; p <= 4; ++p)
+   {
+      const int dim = mesh.Dimension();
+      ND_FECollection fec(p, dim);
+      FiniteElementSpace fes(&mesh, &fec);
+
+      const FiniteElement &el = *fes.GetFE(0);
+      IsoparametricTransformation trans;
+      mesh.GetElementTransformation(0, &trans);
+
+      int order = trans.OrderW() + 2 * el.GetOrder();
+      const IntegrationRule *ir = &IntRules.Get(el.GetGeomType(), order);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+
+         const int dof = el.GetDof();
+         const int el_dim = el.GetDim();
+         DenseMatrix curlshape_bar(dof, el_dim);
+         for (int k = 0; k < curlshape_bar.Width(); ++k)
+         {
+            for (int j = 0; j < curlshape_bar.Height(); ++j)
+            {
+               curlshape_bar(j, k) = distribution(generator);
+            }
+         }
+
+         // reverse-mode differentiation CalcPhysCurlShape
+         DenseMatrix &coords = trans.GetPointMat();
+         DenseMatrix coords_bar(coords.Height(), coords.Width());
+         coords_bar = 0.0;
+         el.CalcPhysCurlShapeRevDiff(trans, curlshape_bar, coords_bar);
+
+         // get the weighted derivatives using finite difference method
+         DenseMatrix curlshape_fd(dof, el_dim), curlshape_pert(dof, el_dim);
+         for (int n = 0; n < coords.Width(); ++n)
+         {
+            for (int di = 0; di < coords.Height(); ++di)
+            {
+               coords(di, n) += eps_fd;
+               trans.Reset();
+               el.CalcPhysCurlShape(trans, curlshape_fd);
+               coords(di, n) -= 2.0 * eps_fd;
+               trans.Reset();
+               el.CalcPhysCurlShape(trans, curlshape_pert);
+               curlshape_fd -= curlshape_pert;
+               curlshape_fd *= 1.0 / (2.0 * eps_fd);
+               coords(di, n) += eps_fd;
+               double x_bar_fd = curlshape_bar * curlshape_fd;
+
+               REQUIRE(coords_bar(di, n) == Approx(x_bar_fd));
+            }
+         }
+      }
+   }
+}
+
+} // namespace fe_revdiff

tests/unit/linalg/test_matrix_dense.cpp

@@ -341,3 +341,517 @@ TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
       REQUIRE(t3.Data()[i] == t1.Data()[i]);
    }
 }
+
+TEST_CASE("DenseMatrix CalcAdjugateRevDiff", "[DenseMatrix]")
+{
+   constexpr double eps_fd = 1e-5; // 2nd-order finite-difference step size
+
+   SECTION("1x1 matrix")
+   {
+      double A_data[1] = { 3.1415926};
+      double adjA_bar_data[1] = {-2.0};
+
+      DenseMatrix A(A_data, 1, 1);
+      DenseMatrix adjA_bar(adjA_bar_data, 1, 1);
+      DenseMatrix A_bar(1, 1), adjA_fd(1, 1);
+      DenseMatrix A_pert(1, 1), adjA_pert(1, 1);
+
+      // Compute the derivative using reverse mode
+      CalcAdjugateRevDiff(A, adjA_bar, A_bar);
+
+      // Compute the derivative using central finite-difference approximation
+      A_pert = A;
+      A_pert(0, 0) += eps_fd;
+      CalcAdjugate(A_pert, adjA_fd);
+      A_pert(0, 0) -= 2.0 * eps_fd;
+      CalcAdjugate(A_pert, adjA_pert);
+      adjA_fd -= adjA_pert;
+      adjA_fd *= 1/(2.0 * eps_fd);
+      // sum up derivative with weights
+      double A_bar_fd = adjA_fd(0, 0) * adjA_bar(0, 0);
+      REQUIRE(A_bar(0, 0) == Approx(A_bar_fd));
+   }
+
+   SECTION("2x1 matrix")
+   {
+      double A_data[2] = {2.0, -3.0};
+      double adjA_bar_data[2] = {-1.5, 4.0};
+
+      DenseMatrix A(A_data, 2, 1);
+      DenseMatrix adjA_bar(adjA_bar_data, 1, 2);
+      DenseMatrix A_bar(2,1), adjA_fd(1,2);
+      DenseMatrix A_pert(2,1), adjA_pert(1,2);
+
+      // Compute the derivatives using reverse mode
+      CalcAdjugateRevDiff(A, adjA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 2; ++i)
+      {
+         // Pertrub A(i,0) and evaluate derivative of adjugate
+         A_pert = A;
+         A_pert(i, 0) += eps_fd;
+         CalcAdjugate(A_pert, adjA_fd);
+         A_pert(i, 0) -= 2.0 * eps_fd;
+         CalcAdjugate(A_pert, adjA_pert);
+         adjA_fd -= adjA_pert;
+         adjA_fd *= 1 / (2.0 * eps_fd);
+         // sum up derivative with weights
+         double A_bar_fd = 0.0;
+         for (int k = 0; k < 2; ++k)
+         {
+            A_bar_fd += adjA_fd(0, k) * adjA_bar(0, k);
+         }
+         REQUIRE(A_bar(i, 0) == Approx(A_bar_fd));
+      }
+   }
+
+   SECTION("2x2 matrix")
+   {
+      double A_data[4] = {2.0, -3.0, 4.0, -1.0};
+      double adjA_bar_data[4] = {1.0, 4.0, 2.0, -3.0};
+
+      DenseMatrix A(A_data, 2, 2);
+      DenseMatrix adjA_bar(adjA_bar_data, 2, 2);
+      DenseMatrix A_bar(2,2), adjA_fd(2,2);
+      DenseMatrix A_pert(2,2), adjA_pert(2,2);
+
+      // Compute the derivatives using reverse mode
+      CalcAdjugateRevDiff(A, adjA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 2; ++i)
+      {
+         for (int j = 0; j < 2; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of adjugate
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcAdjugate(A_pert, adjA_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcAdjugate(A_pert, adjA_pert);
+            adjA_fd -= adjA_pert;
+            adjA_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 2; ++k)
+            {
+               for (int l = 0; l < 2; ++l)
+               {
+                  A_bar_fd += adjA_fd(k,l)*adjA_bar(k,l);
+               }
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+
+   SECTION("3x1 matrix")
+   {
+      double A_data[3] = {2.0, -3.0, 3.1415926};
+      double adjA_bar_data[3] = {-1.5, 4.0, 2.71828};
+
+      DenseMatrix A(A_data, 3, 1);
+      DenseMatrix adjA_bar(adjA_bar_data, 1, 3);
+      DenseMatrix A_bar(3, 1), adjA_fd(1, 3);
+      DenseMatrix A_pert(3, 1), adjA_pert(1, 3);
+
+      // Compute the derivatives using reverse mode
+      CalcAdjugateRevDiff(A, adjA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         // Pertrub A(i,0) and evaluate derivative of adjugate
+         A_pert = A;
+         A_pert(i, 0) += eps_fd;
+         CalcAdjugate(A_pert, adjA_fd);
+         A_pert(i, 0) -= 2.0 * eps_fd;
+         CalcAdjugate(A_pert, adjA_pert);
+         adjA_fd -= adjA_pert;
+         adjA_fd *= 1 / (2.0 * eps_fd);
+         // sum up derivative with weights
+         double A_bar_fd = 0.0;
+         for (int k = 0; k < 3; ++k)
+         {
+            A_bar_fd += adjA_fd(0, k) * adjA_bar(0, k);
+         }
+         REQUIRE(A_bar(i, 0) == Approx(A_bar_fd));
+      }
+   }
+
+   SECTION("3x3 matrix")
+   {
+      double A_data[9] = {1.0, 5.0, 3.0, -2.0, 6.0, -9.0, 4.0, -7.0, 8.0};
+      double adjA_bar_data[9] = {3.0, 6.0, -8.0, 1.0, -7.0, 5.0, 2.0, 4.0, -9.0};
+
+      DenseMatrix A(A_data, 3, 3);
+      DenseMatrix adjA_bar(adjA_bar_data, 3, 3);
+      DenseMatrix A_bar(3,3), adjA_fd(3,3);
+      DenseMatrix A_pert(3,3), adjA_pert(3,3);
+
+      // Compute the derivatives using reverse mode
+      CalcAdjugateRevDiff(A, adjA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         for (int j = 0; j < 3; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of adjugate
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcAdjugate(A_pert, adjA_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcAdjugate(A_pert, adjA_pert);
+            adjA_fd -= adjA_pert;
+            adjA_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 3; ++k)
+            {
+               for (int l = 0; l < 3; ++l)
+               {
+                  A_bar_fd += adjA_fd(k,l)*adjA_bar(k,l);
+               }
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+}
+
+TEST_CASE("DenseMatrix CalcInverseRevDiff", "[DenseMatrix]")
+{
+   constexpr double eps_fd = 1e-5; // 2nd-order finite-difference step size
+
+   SECTION("1x1 matrix")
+   {
+      double A_data[1] = { 3.1415926};
+      double invA_bar_data[1] = {-2.0};
+
+      DenseMatrix A(A_data, 1, 1);
+      DenseMatrix invA_bar(invA_bar_data, 1, 1);
+      DenseMatrix A_bar(1, 1), invA_fd(1, 1);
+      DenseMatrix A_pert(1, 1), invA_pert(1, 1);
+
+      // Compute the derivative using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivative using central finite-difference approximation
+      A_pert = A;
+      A_pert(0, 0) += eps_fd;
+      CalcInverse(A_pert, invA_fd);
+      A_pert(0, 0) -= 2.0 * eps_fd;
+      CalcInverse(A_pert, invA_pert);
+      invA_fd -= invA_pert;
+      invA_fd *= 1/(2.0 * eps_fd);
+      // sum up derivative with weights
+      double A_bar_fd = invA_fd(0, 0) * invA_bar(0, 0);
+      REQUIRE(A_bar(0, 0) == Approx(A_bar_fd));
+   }
+
+   SECTION("2x1 matrix")
+   {
+      double A_data[2] = {2.0, -3.0};
+      double invA_bar_data[2] = {-1.5, 4.0};
+
+      DenseMatrix A(A_data, 2, 1);
+      DenseMatrix invA_bar(invA_bar_data, 1, 2);
+      DenseMatrix A_bar(2,1), invA_fd(1,2);
+      DenseMatrix A_pert(2,1), invA_pert(1,2);
+
+      // Compute the derivatives using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 2; ++i)
+      {
+         // Pertrub A(i,0) and evaluate derivative of inverse
+         A_pert = A;
+         A_pert(i, 0) += eps_fd;
+         CalcInverse(A_pert, invA_fd);
+         A_pert(i, 0) -= 2.0 * eps_fd;
+         CalcInverse(A_pert, invA_pert);
+         invA_fd -= invA_pert;
+         invA_fd *= 1 / (2.0 * eps_fd);
+         // sum up derivative with weights
+         double A_bar_fd = 0.0;
+         for (int k = 0; k < 2; ++k)
+         {
+            A_bar_fd += invA_fd(0, k) * invA_bar(0, k);
+         }
+         REQUIRE(A_bar(i, 0) == Approx(A_bar_fd));
+      }
+   }
+
+   SECTION("2x2 matrix")
+   {
+      double A_data[4] = {2.0, -3.0, 4.0, -1.0};
+      double invA_bar_data[4] = {1.0, 4.0, 2.0, -3.0};
+
+      DenseMatrix A(A_data, 2, 2);
+      DenseMatrix invA_bar(invA_bar_data, 2, 2);
+      DenseMatrix A_bar(2,2), invA_fd(2,2);
+      DenseMatrix A_pert(2,2), invA_pert(2,2);
+
+      // Compute the derivatives using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 2; ++i)
+      {
+         for (int j = 0; j < 2; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of inverse
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcInverse(A_pert, invA_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcInverse(A_pert, invA_pert);
+            invA_fd -= invA_pert;
+            invA_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 2; ++k)
+            {
+               for (int l = 0; l < 2; ++l)
+               {
+                  A_bar_fd += invA_fd(k,l)*invA_bar(k,l);
+               }
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+
+   SECTION("3x1 matrix")
+   {
+      double A_data[3] = {2.0, -3.0, 3.1415926};
+      double invA_bar_data[3] = {-1.5, 4.0, 2.71828};
+
+      DenseMatrix A(A_data, 3, 1);
+      DenseMatrix invA_bar(invA_bar_data, 1, 3);
+      DenseMatrix A_bar(3, 1), invA_fd(1, 3);
+      DenseMatrix A_pert(3, 1), invA_pert(1, 3);
+
+      // Compute the derivatives using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         // Pertrub A(i,0) and evaluate derivative of inverse
+         A_pert = A;
+         A_pert(i, 0) += eps_fd;
+         CalcInverse(A_pert, invA_fd);
+         A_pert(i, 0) -= 2.0 * eps_fd;
+         CalcInverse(A_pert, invA_pert);
+         invA_fd -= invA_pert;
+         invA_fd *= 1 / (2.0 * eps_fd);
+         // sum up derivative with weights
+         double A_bar_fd = 0.0;
+         for (int k = 0; k < 3; ++k)
+         {
+            A_bar_fd += invA_fd(0, k) * invA_bar(0, k);
+         }
+         REQUIRE(A_bar(i, 0) == Approx(A_bar_fd));
+      }
+   }
+
+   SECTION("3x2 matrix")
+   {
+      double A_data[6] = {1.0, 5.0, 3.0, -2.0, 6.0, -9.0};
+      double invA_bar_data[6] = {3.0, 6.0, -8.0, 1.0, -7.0, 5.0};
+      // double invA_bar_data[6] = {1.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+      DenseMatrix invA_bar(invA_bar_data, 2, 3);
+      DenseMatrix invA_pert(2, 3), invA_fd(2, 3);
+      DenseMatrix A(A_data, 3, 2);
+      DenseMatrix A_bar(3, 2), A_pert(3, 2);
+
+      // Compute the derivatives using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         for (int j = 0; j < 2; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of inverse
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcInverse(A_pert, invA_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcInverse(A_pert, invA_pert);
+            invA_fd -= invA_pert;
+            invA_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 3; ++k)
+            {
+               for (int l = 0; l < 2; ++l)
+               {
+                  A_bar_fd += invA_fd(l, k)*invA_bar(l, k);
+               }
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+
+   SECTION("3x3 matrix")
+   {
+      double A_data[9] = {1.0, 5.0, 3.0, -2.0, 6.0, -9.0, 4.0, -7.0, 8.0};
+      double invA_bar_data[9] = {3.0, 6.0, -8.0, 1.0, -7.0, 5.0, 2.0, 4.0, -9.0};
+
+      DenseMatrix A(A_data, 3, 3);
+      DenseMatrix invA_bar(invA_bar_data, 3, 3);
+      DenseMatrix A_bar(3,3), invA_fd(3,3);
+      DenseMatrix A_pert(3,3), invA_pert(3,3);
+
+      // Compute the derivatives using reverse mode
+      A_bar = 0.0;
+      CalcInverseRevDiff(A, invA_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         for (int j = 0; j < 3; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of inverse
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcInverse(A_pert, invA_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcInverse(A_pert, invA_pert);
+            invA_fd -= invA_pert;
+            invA_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 3; ++k)
+            {
+               for (int l = 0; l < 3; ++l)
+               {
+                  A_bar_fd += invA_fd(k,l)*invA_bar(k,l);
+               }
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+}
+
+TEST_CASE("DenseMatrix WeightRevDiff", "[DenseMatrix]")
+{
+   // This also tests DenseMatrix::DetRevDiff indirectly
+   constexpr double eps_fd = 1e-5; // 2nd-order finite-difference step size
+   double A_data[9] = {1.0, 5.0, 3.0, -2.0, 6.0, -9.0, 4.0, -7.0, 8.0};
+
+   for (int height = 1; height <= 3; ++height)
+   {
+      for (int width = 1; width <= height; ++width)
+      {
+         DenseMatrix A(A_data, height, width);
+         DenseMatrix weight_bar(height, width);
+         DenseMatrix A_pert(height, width);
+
+         // Compute the gradient of A.Weight() using reverse mode AD
+         A.WeightRevDiff(weight_bar);
+
+         // Compute the gradient of A.Weight using 2nd order finite-difference
+         for (int i = 0; i < height; ++i)
+         {
+            for (int j = 0; j < width; ++j)
+            {
+               // Perturb A(i,j) in + and - directions and evaluate Weight()
+               A_pert = A;
+               A_pert(i,j) += eps_fd;
+               double dweight = A_pert.Weight();
+               A_pert(i,j) -= 2.0*eps_fd;
+               dweight -= A_pert.Weight();
+               dweight /= (2.0*eps_fd);
+               REQUIRE(weight_bar(i,j) == Approx(dweight));
+            }
+         }
+      }
+   }
+}
+
+TEST_CASE("DenseMatrix CalcOrthoRevDiff", "[DenseMatrix]")
+{
+   constexpr double eps_fd = 1e-5; // 2nd-order finite-difference step size
+
+   SECTION("2x1 matrix")
+   {
+      double A_data[2] = {2.0, -3.0};
+      double n_bar_data[2] = {-1.5, 4.0};
+      Vector n_bar(n_bar_data, 2);
+      Vector n_pert(2), n_fd(2);
+      DenseMatrix A(A_data, 2, 1);
+      DenseMatrix A_bar(2,1), A_pert(2,1);
+
+      // Compute the derivatives using reverse mode
+      CalcOrthoRevDiff(A, n_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 2; ++i)
+      {
+         // Pertrub A(i,0) and evaluate derivative of adjugate
+         A_pert = A;
+         A_pert(i, 0) += eps_fd;
+         CalcOrtho(A_pert, n_fd);
+         A_pert(i, 0) -= 2.0 * eps_fd;
+         CalcOrtho(A_pert, n_pert);
+         n_fd -= n_pert;
+         n_fd *= 1 / (2.0 * eps_fd);
+         // sum up derivative with weights
+         double A_bar_fd = 0.0;
+         for (int k = 0; k < 2; ++k)
+         {
+            A_bar_fd += n_fd(k) * n_bar(k);
+         }
+         REQUIRE(A_bar(i, 0) == Approx(A_bar_fd));
+      }
+   }
+
+   SECTION("3x2 matrix")
+   {
+      double A_data[6] = {1.0, 5.0, 3.0, -2.0, 6.0, -9.0};
+      double n_bar_data[3] = {1.0, 4.0, -3.0};
+      Vector n_bar(n_bar_data, 3);
+      Vector n_pert(3), n_fd(3);
+      DenseMatrix A(A_data, 3, 2);
+      DenseMatrix A_bar(3,2), A_pert(3,2);
+
+      // Compute the derivatives using reverse mode
+      CalcOrthoRevDiff(A, n_bar, A_bar);
+
+      // Compute the derivatives using central finite-difference approximation
+      for (int i = 0; i < 3; ++i)
+      {
+         for (int j = 0; j < 2; ++j)
+         {
+            // Pertrub A(i,j) and evaluate derivative of adjugate
+            A_pert = A;
+            A_pert(i,j) += eps_fd;
+            CalcOrtho(A_pert, n_fd);
+            A_pert(i,j) -= 2.0*eps_fd;
+            CalcOrtho(A_pert, n_pert);
+            n_fd -= n_pert;
+            n_fd *= 1/(2.0*eps_fd);
+            // sum up derivative with weights
+            double A_bar_fd = 0.0;
+            for (int k = 0; k < 3; ++k)
+            {
+               A_bar_fd += n_fd(k)*n_bar(k);
+            }
+            REQUIRE(A_bar(i,j) == Approx(A_bar_fd));
+         }
+      }
+   }
+}