From 00bf53ed901020ea8b0950524e6a29bdb7f965a1 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 4 Apr 2024 15:48:32 -0400
Subject: [PATCH 01/17] introduce the ExponentialMatrixCoefficient class

---
 fem/bilininteg.hpp                      |  6 +--
 fem/coefficient.cpp                     | 24 ++++++++++
 fem/coefficient.hpp                     | 27 +++++++++++-
 linalg/densemat.cpp                     | 50 +++++++++++++++++++++
 linalg/densemat.hpp                     |  4 ++
 tests/unit/linalg/test_matrix_dense.cpp | 58 +++++++++++++++++++++++++
 6 files changed, 164 insertions(+), 5 deletions(-)

diff --git a/fem/bilininteg.hpp b/fem/bilininteg.hpp
index 8dff5e2f04..04868cc08d 100644
--- a/fem/bilininteg.hpp
+++ b/fem/bilininteg.hpp
@@ -1708,7 +1708,7 @@ public:
    { vector_fe.CalcPhysDShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot  \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot \nabla u, \nabla \cdot v)$ in 2D
     or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ and $v$ is in $H(div)$. */
 class MixedGradDivIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1747,7 +1747,7 @@ public:
    { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla v)$ in 2D
     or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H(div)$ and $v$ is in $H^1$. */
 class MixedDivGradIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1787,7 +1787,7 @@ public:
    { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla \cdot v)$ in 2D or 3D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla v)$ in 2D or 3D
     and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ or $L_2$ and $v$ is in $H^1$. */
 class MixedScalarWeakDivergenceIntegrator : public MixedScalarVectorIntegrator
 {
diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index a4d7ad2c77..a3878fa139 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -1326,6 +1326,30 @@ void InverseMatrixCoefficient::Eval(DenseMatrix &M,
    M.Invert();
 }
 
+ExponentialMatrixCoefficient::ExponentialMatrixCoefficient(MatrixCoefficient &A)
+   : MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth() && A.GetHeight() == 2,
+               "ExponentialMatrixCoefficient:  "
+               << "Argument must be a square 2x2 matrix."
+               << "  Height = " << A.GetHeight()
+               << ", Width = " << A.GetWidth());
+}
+
+void ExponentialMatrixCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void ExponentialMatrixCoefficient::Eval(DenseMatrix &M,
+                                    ElementTransformation &T,
+                                    const IntegrationPoint &ip)
+{
+   a->Eval(M, T, ip);
+   M.Exponential();
+}
+
 OuterProductCoefficient::OuterProductCoefficient(VectorCoefficient &A,
                                                  VectorCoefficient &B)
    : MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),
diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index ed899d545b..2529ba7589 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -2112,7 +2112,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the transpose a matrix coefficient
+/// Matrix coefficient defined as the transpose of a matrix coefficient
 class TransposeMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2135,7 +2135,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the inverse a matrix coefficient.
+/// Matrix coefficient defined as the inverse of a matrix coefficient.
 class InverseMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2158,6 +2158,29 @@ public:
                      const IntegrationPoint &ip);
 };
 
+/// Matrix coefficient defined as the exponential of a matrix coefficient.
+class ExponentialMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   MatrixCoefficient * a;
+
+public:
+   /// Construct the matrix coefficient.  Result is $ \exp(A) $.
+   ExponentialMatrixCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the matrix coefficient at @a ip.
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
 /// Matrix coefficient defined as the outer product of two vector coefficients.
 class OuterProductCoefficient : public MatrixCoefficient
 {
diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 2f46fa4c3e..8848f0f77c 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -532,6 +532,56 @@ MatrixInverse *DenseMatrix::Inverse() const
    return new DenseMatrixInverse(*this);
 }
 
+void DenseMatrix::Exponential()
+{
+   MFEM_ASSERT(Height() == Width() && Height() <= 2,
+               "The matrix must be square and "
+               << "of size less than or equal to 2."
+               << "  Height() = " << Height()
+               << ", Width() = " << Width());
+
+   switch (Height())
+   {
+      case 1:
+         data[0] = std::exp(data[0]);
+
+      case 2:
+         /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
+         /// Note typo in the paper, in the prefactor in the equation under (i)·
+         real_t a = data[0];
+         real_t b = data[1];
+         real_t c = data[2];
+         real_t d = data[3];
+         real_t e = (a - d)*(a - d) + 4*b*c;
+         real_t f = std::exp((a + d)/2.0);
+         real_t g = std::sqrt(std::abs(e)) / 2.0;
+
+         if (e == 0)
+         {
+            data[0] = 1.0 + (a - d)/2.0;
+            data[3] = 1.0 - (a - d)/2.0;
+         }
+         else if (e > 0)
+         {
+            data[0] = std::cosh(g) + (a - d)/2 * std::sinh(g) / g;
+            data[1] = b * std::sinh(g) / g;
+            data[2] = c * std::sinh(g) / g;
+            data[3] = std::cosh(g) - (a - d)/2 * std::sinh(g) / g;
+         }
+         else
+         {
+            data[0] = std::cos(g) + (a - d)/2 * std::sin(g) / g;
+            data[1] = b * std::sin(g) / g;
+            data[2] = c * std::sin(g) / g;
+            data[3] = std::cos(g) - (a - d)/2 * std::sin(g) / g;
+         }
+         for (int i = 0; i < 4; i++)
+         {
+            data[i] *= f;
+         }
+   }
+}
+
 real_t DenseMatrix::Det() const
 {
    MFEM_ASSERT(Height() == Width() && Height() > 0,
diff --git a/linalg/densemat.hpp b/linalg/densemat.hpp
index b0ce139e04..7c6402bda4 100644
--- a/linalg/densemat.hpp
+++ b/linalg/densemat.hpp
@@ -207,6 +207,10 @@ public:
    /// Replaces the current matrix with its square root inverse
    void SquareRootInverse();
 
+   /// Replaces the current matrix with its exponential
+   /// (currently only supports 2x2 matrices)
+   void Exponential();
+
    /// Calculates the determinant of the matrix
    /// (optimized for 2x2, 3x3, and 4x4 matrices)
    real_t Det() const;
diff --git a/tests/unit/linalg/test_matrix_dense.cpp b/tests/unit/linalg/test_matrix_dense.cpp
index 53bc1f3fd4..9f64216b9f 100644
--- a/tests/unit/linalg/test_matrix_dense.cpp
+++ b/tests/unit/linalg/test_matrix_dense.cpp
@@ -719,4 +719,62 @@ TEST_CASE("NNLS", "[DenseMatrix]")
    REQUIRE(sol[4] == MFEM_Approx(2.5));
 }
 
+TEST_CASE("Exponential", "[DenseMatrix]")
+{
+   // case 1
+   DenseMatrix A(2,2);
+   A(0,0) = 5.0;
+   A(0,1) = 3.0;
+   A(1,0) = 0.0;
+   A(1,1) = 5.0;
+   A.Exponential();
+
+   DenseMatrix expA(2,2);
+   expA(0,0) = std::exp(5.0);
+   expA(0,1) = 3.0 * std::exp(5.0);
+   expA(1,0) = 0.0;
+   expA(1,1) = std::exp(5.0);
+
+   A.Print();
+   expA.Print();
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 2
+   A(0,0) = 3.0;
+   A(0,1) = 5.0;
+   A(1,0) = 4.0;
+   A(1,1) = 2.0;
+   A.Exponential();
+
+   expA(0,0) = 4.0 / (9.0 * std::exp(2.0)) + (5.0 * std::exp(7.0)) / 9.0;
+   expA(0,1) = (5.0 * std::exp(7.0)) / 9.0 - 5.0 / (9.0 * std::exp(2.0));
+   expA(1,0) = (4.0 * std::exp(7.0)) / 9.0 - 4.0 / (9.0 * std::exp(2.0));
+   expA(1,1) = 5.0 / (9.0 * std::exp(2.0)) + (4.0 * std::exp(7.0)) / 9.0;
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 3
+   A(0,0) = 10.0;
+   A(0,1) = 2.0;
+   A(1,0) = -2.0;
+   A(1,1) = 8.0;
+   A.Exponential();
+
+   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0) + std::cos(std::sqrt(3.0)));
+   expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0))  / std::sqrt(3.0);
+   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0)) - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+}
+
 #endif // if MFEM_USE_LAPACK

From d453981d3cbb8954776a9e5966f566b95ea7b140 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 4 Apr 2024 16:09:13 -0400
Subject: [PATCH 02/17] style

---
 fem/coefficient.cpp                     | 4 ++--
 tests/unit/linalg/test_matrix_dense.cpp | 6 ++++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index a3878fa139..54ecf0524b 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -1343,8 +1343,8 @@ void ExponentialMatrixCoefficient::SetTime(real_t t)
 }
 
 void ExponentialMatrixCoefficient::Eval(DenseMatrix &M,
-                                    ElementTransformation &T,
-                                    const IntegrationPoint &ip)
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
 {
    a->Eval(M, T, ip);
    M.Exponential();
diff --git a/tests/unit/linalg/test_matrix_dense.cpp b/tests/unit/linalg/test_matrix_dense.cpp
index 9f64216b9f..138ec5e808 100644
--- a/tests/unit/linalg/test_matrix_dense.cpp
+++ b/tests/unit/linalg/test_matrix_dense.cpp
@@ -766,10 +766,12 @@ TEST_CASE("Exponential", "[DenseMatrix]")
    A(1,1) = 8.0;
    A.Exponential();
 
-   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0) + std::cos(std::sqrt(3.0)));
+   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(
+                                   3.0) + std::cos(std::sqrt(3.0)));
    expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
    expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0))  / std::sqrt(3.0);
-   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0)) - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
+   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0)) - std::sin(std::sqrt(
+                                                                       3.0)) / std::sqrt(3.0));
 
    REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
    REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));

From 1f5f30c9c4ae9ce36dd49b03ed10f5839a63bbcf Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 4 Apr 2024 16:27:45 -0400
Subject: [PATCH 03/17] 3x3 abort message

---
 linalg/densemat.cpp | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 8848f0f77c..1d1ed4903b 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -547,7 +547,7 @@ void DenseMatrix::Exponential()
 
       case 2:
          /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
-         /// Note typo in the paper, in the prefactor in the equation under (i)·
+         /// Note typo in the paper, in the prefactor in the equation under (i).
          real_t a = data[0];
          real_t b = data[1];
          real_t c = data[2];
@@ -579,6 +579,10 @@ void DenseMatrix::Exponential()
          {
             data[i] *= f;
          }
+      case 3:
+         MFEM_ABORT("3x3 matrices are not currently supported");
+      default:
+         MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
    }
 }
 

From 4dcb5933a9eea4cecc78cc6210dc97b3eaa89cd6 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 4 Apr 2024 17:11:32 -0400
Subject: [PATCH 04/17] bug in switch

---
 linalg/densemat.cpp | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 1d1ed4903b..77c2cff2e7 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -543,9 +543,11 @@ void DenseMatrix::Exponential()
    switch (Height())
    {
       case 1:
+      {
          data[0] = std::exp(data[0]);
-
+      }
       case 2:
+      {
          /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
          /// Note typo in the paper, in the prefactor in the equation under (i).
          real_t a = data[0];
@@ -579,10 +581,15 @@ void DenseMatrix::Exponential()
          {
             data[i] *= f;
          }
+      }
       case 3:
+      {
          MFEM_ABORT("3x3 matrices are not currently supported");
+      }
       default:
+      {
          MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
+      }
    }
 }
 

From 99db13a3c27e386876139f171765dc68819c93a3 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 4 Apr 2024 17:22:06 -0400
Subject: [PATCH 05/17] missing break;

---
 linalg/densemat.cpp                     |  2 ++
 tests/unit/linalg/test_matrix_dense.cpp | 10 +++++-----
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 77c2cff2e7..4e87cfe836 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -545,6 +545,7 @@ void DenseMatrix::Exponential()
       case 1:
       {
          data[0] = std::exp(data[0]);
+         break;
       }
       case 2:
       {
@@ -581,6 +582,7 @@ void DenseMatrix::Exponential()
          {
             data[i] *= f;
          }
+         break;
       }
       case 3:
       {
diff --git a/tests/unit/linalg/test_matrix_dense.cpp b/tests/unit/linalg/test_matrix_dense.cpp
index 138ec5e808..aebc031abb 100644
--- a/tests/unit/linalg/test_matrix_dense.cpp
+++ b/tests/unit/linalg/test_matrix_dense.cpp
@@ -766,12 +766,12 @@ TEST_CASE("Exponential", "[DenseMatrix]")
    A(1,1) = 8.0;
    A.Exponential();
 
-   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(
-                                   3.0) + std::cos(std::sqrt(3.0)));
+   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0)
+                                + std::cos(std::sqrt(3.0)));
    expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
-   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0))  / std::sqrt(3.0);
-   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0)) - std::sin(std::sqrt(
-                                                                       3.0)) / std::sqrt(3.0));
+   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0))
+                                - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
 
    REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
    REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));

From bc0ab53d1937f6e2a3001fb7587acaed858beded Mon Sep 17 00:00:00 2001
From: Brendan Keith <keith10@llnl.gov>
Date: Fri, 5 Apr 2024 08:45:33 -0400
Subject: [PATCH 06/17] Update linalg/densemat.cpp

Co-authored-by: Dohyun Kim <dhkim.cse@gmail.com>
---
 linalg/densemat.cpp | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 4e87cfe836..17673f4cf7 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -551,13 +551,13 @@ void DenseMatrix::Exponential()
       {
          /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
          /// Note typo in the paper, in the prefactor in the equation under (i).
-         real_t a = data[0];
-         real_t b = data[1];
-         real_t c = data[2];
-         real_t d = data[3];
-         real_t e = (a - d)*(a - d) + 4*b*c;
-         real_t f = std::exp((a + d)/2.0);
-         real_t g = std::sqrt(std::abs(e)) / 2.0;
+         const real_t a = data[0];
+         const real_t b = data[1];
+         const real_t c = data[2];
+         const real_t d = data[3];
+         const real_t e = (a - d)*(a - d) + 4*b*c;
+         const real_t f = std::exp((a + d)/2.0);
+         const real_t g = std::sqrt(std::abs(e)) / 2.0;
 
          if (e == 0)
          {

From 18bee592c419b635679ff5e19010832ec9b3870a Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Fri, 5 Apr 2024 12:38:35 -0400
Subject: [PATCH 07/17] move test out of #ifdef MFEM_USE_LAPACK

---
 tests/unit/linalg/test_matrix_dense.cpp | 122 ++++++++++++------------
 1 file changed, 62 insertions(+), 60 deletions(-)

diff --git a/tests/unit/linalg/test_matrix_dense.cpp b/tests/unit/linalg/test_matrix_dense.cpp
index aebc031abb..9047ce1af2 100644
--- a/tests/unit/linalg/test_matrix_dense.cpp
+++ b/tests/unit/linalg/test_matrix_dense.cpp
@@ -532,6 +532,68 @@ TEST_CASE("MatrixInverse", "[DenseMatrix]")
    }
 
 }
+
+TEST_CASE("Exponential", "[DenseMatrix]")
+{
+   // case 1
+   DenseMatrix A(2,2);
+   A(0,0) = 5.0;
+   A(0,1) = 3.0;
+   A(1,0) = 0.0;
+   A(1,1) = 5.0;
+   A.Exponential();
+
+   DenseMatrix expA(2,2);
+   expA(0,0) = std::exp(5.0);
+   expA(0,1) = 3.0 * std::exp(5.0);
+   expA(1,0) = 0.0;
+   expA(1,1) = std::exp(5.0);
+
+   A.Print();
+   expA.Print();
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 2
+   A(0,0) = 3.0;
+   A(0,1) = 5.0;
+   A(1,0) = 4.0;
+   A(1,1) = 2.0;
+   A.Exponential();
+
+   expA(0,0) = 4.0 / (9.0 * std::exp(2.0)) + (5.0 * std::exp(7.0)) / 9.0;
+   expA(0,1) = (5.0 * std::exp(7.0)) / 9.0 - 5.0 / (9.0 * std::exp(2.0));
+   expA(1,0) = (4.0 * std::exp(7.0)) / 9.0 - 4.0 / (9.0 * std::exp(2.0));
+   expA(1,1) = 5.0 / (9.0 * std::exp(2.0)) + (4.0 * std::exp(7.0)) / 9.0;
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 3
+   A(0,0) = 10.0;
+   A(0,1) = 2.0;
+   A(1,0) = -2.0;
+   A(1,1) = 8.0;
+   A.Exponential();
+
+   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0)
+                                + std::cos(std::sqrt(3.0)));
+   expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0))
+                                - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+}
+
+
 #ifdef MFEM_USE_LAPACK
 
 enum class TestCase { GenEigSPD, GenEigGE, SVD};
@@ -719,64 +781,4 @@ TEST_CASE("NNLS", "[DenseMatrix]")
    REQUIRE(sol[4] == MFEM_Approx(2.5));
 }
 
-TEST_CASE("Exponential", "[DenseMatrix]")
-{
-   // case 1
-   DenseMatrix A(2,2);
-   A(0,0) = 5.0;
-   A(0,1) = 3.0;
-   A(1,0) = 0.0;
-   A(1,1) = 5.0;
-   A.Exponential();
-
-   DenseMatrix expA(2,2);
-   expA(0,0) = std::exp(5.0);
-   expA(0,1) = 3.0 * std::exp(5.0);
-   expA(1,0) = 0.0;
-   expA(1,1) = std::exp(5.0);
-
-   A.Print();
-   expA.Print();
-   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
-   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
-   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
-   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
-
-   // case 2
-   A(0,0) = 3.0;
-   A(0,1) = 5.0;
-   A(1,0) = 4.0;
-   A(1,1) = 2.0;
-   A.Exponential();
-
-   expA(0,0) = 4.0 / (9.0 * std::exp(2.0)) + (5.0 * std::exp(7.0)) / 9.0;
-   expA(0,1) = (5.0 * std::exp(7.0)) / 9.0 - 5.0 / (9.0 * std::exp(2.0));
-   expA(1,0) = (4.0 * std::exp(7.0)) / 9.0 - 4.0 / (9.0 * std::exp(2.0));
-   expA(1,1) = 5.0 / (9.0 * std::exp(2.0)) + (4.0 * std::exp(7.0)) / 9.0;
-
-   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
-   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
-   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
-   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
-
-   // case 3
-   A(0,0) = 10.0;
-   A(0,1) = 2.0;
-   A(1,0) = -2.0;
-   A(1,1) = 8.0;
-   A.Exponential();
-
-   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0)
-                                + std::cos(std::sqrt(3.0)));
-   expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
-   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
-   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0))
-                                - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
-
-   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
-   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
-   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
-   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
-}
-
 #endif // if MFEM_USE_LAPACK

From b0a3350622f488cf6eb1d4c29b25da5d1947fe01 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Fri, 5 Apr 2024 17:40:59 -0400
Subject: [PATCH 08/17] added MatrixArrayVectorCoefficient clas

---
 fem/coefficient.cpp                     | 66 +++++++++++++++++++++++++
 fem/coefficient.hpp                     | 39 +++++++++++++++
 tests/unit/linalg/test_matrix_dense.cpp |  1 -
 3 files changed, 105 insertions(+), 1 deletion(-)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index 54ecf0524b..4cbed5daf5 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -924,6 +924,72 @@ void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
    }
 }
 
+MatrixArrayVectorCoefficient::MatrixArrayVectorCoefficient (int dim)
+   : MatrixCoefficient (dim)
+{
+   Coeff.SetSize(height);
+   ownCoeff.SetSize(height);
+   for (int i = 0; i < (height); i++)
+   {
+      Coeff[i] = NULL;
+      ownCoeff[i] = true;
+   }
+}
+
+void MatrixArrayVectorCoefficient::SetTime(real_t t)
+{
+   for (int i=0; i < height; i++)
+   {
+      if (Coeff[i]) { Coeff[i]->SetTime(t); }
+   }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void MatrixArrayVectorCoefficient::Set(int i, VectorCoefficient * c, bool own)
+{
+   if (ownCoeff[i]) { delete Coeff[i]; }
+   Coeff[i] = c;
+   ownCoeff[i] = own;
+}
+
+MatrixArrayVectorCoefficient::~MatrixArrayVectorCoefficient ()
+{
+   for (int i=0; i < height; i++)
+   {
+      if (ownCoeff[i]) { delete Coeff[i]; }
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(int i, Vector &V,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   if (Coeff[i])
+   {
+      Coeff[i] -> Eval(V, T, ip);
+   }
+   else
+   {
+      V = 0.0;
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   K.SetSize(height, width);
+   for (int i = 0; i < height; i++)
+   {
+      Vector V(width);
+      this->Eval(i, V, T, ip);
+      for (int j = 0; j < width; j++)
+      {
+         K(i,j) = V(j);
+      }
+   }
+}
+
 void MatrixRestrictedCoefficient::SetTime(real_t t)
 {
    if (c) { c->SetTime(t); }
diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 2529ba7589..320b9df3ad 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1334,6 +1334,45 @@ public:
    virtual ~MatrixArrayCoefficient();
 };
 
+/** @brief Matrix coefficient defined by a matrix of vector coefficients.
+    Coefficients that are not set will evaluate to zero in the vector. The
+    coefficient is stored as an Array indexing the rows of the matrix. */
+class MatrixArrayVectorCoefficient : public MatrixCoefficient
+{
+private:
+   Array<VectorCoefficient *> Coeff;
+   Array<bool> ownCoeff;
+
+public:
+   /** @brief Construct a coefficient matrix of dimensions @a dim * @a dim. The
+       actual coefficients still need to be added with Set(). */
+   explicit MatrixArrayVectorCoefficient (int dim);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Get the coefficient located at (i,j) in the matrix.
+   VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
+
+   /** @brief Set the coefficient located at the i-th row of the matrix.
+       By this will take ownership of the Coefficient passed in, but this
+       can be overridden with the @a own parameter. */
+   void Set(int i, VectorCoefficient * c, bool own=true);
+
+   using MatrixCoefficient::Eval;
+
+   /// Evaluate coefficient located at the i-th row of the matrix using integration
+   /// point @a ip.
+   void Eval(int i, Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+   /// Evaluate the matrix coefficient @a ip.
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+   virtual ~MatrixArrayVectorCoefficient();
+};
+
 
 /** @brief Derived matrix coefficient that has the value of the parent matrix
     coefficient where it is active and is zero otherwise. */
diff --git a/tests/unit/linalg/test_matrix_dense.cpp b/tests/unit/linalg/test_matrix_dense.cpp
index 9047ce1af2..2ab423a6d1 100644
--- a/tests/unit/linalg/test_matrix_dense.cpp
+++ b/tests/unit/linalg/test_matrix_dense.cpp
@@ -593,7 +593,6 @@ TEST_CASE("Exponential", "[DenseMatrix]")
    REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
 }
 
-
 #ifdef MFEM_USE_LAPACK
 
 enum class TestCase { GenEigSPD, GenEigGE, SVD};

From 90ecbf2bfb0fa26594fad0cfc77445559bf1a8b3 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Fri, 5 Apr 2024 17:44:29 -0400
Subject: [PATCH 09/17] Adding TraceCoefficient

---
 fem/coefficient.cpp | 21 +++++++++++++++++++++
 fem/coefficient.hpp | 25 +++++++++++++++++++++++++
 2 files changed, 46 insertions(+)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index 4cbed5daf5..864973da5d 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -1107,6 +1107,27 @@ real_t DeterminantCoefficient::Eval(ElementTransformation &T,
    return ma.Det();
 }
 
+TraceCoefficient::TraceCoefficient(MatrixCoefficient &A)
+   : a(&A), ma(A.GetHeight(), A.GetWidth())
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
+               "TraceCoefficient:  "
+               "Argument must be a square matrix.");
+}
+
+void TraceCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+real_t TraceCoefficient::Eval(ElementTransformation &T,
+                              const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   return ma.Trace();
+}
+
 VectorSumCoefficient::VectorSumCoefficient(int dim)
    : VectorCoefficient(dim),
      ACoef(NULL), BCoef(NULL),
diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 320b9df3ad..d6a5743b49 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1800,6 +1800,31 @@ public:
                        const IntegrationPoint &ip);
 };
 
+/// Scalar coefficient defined as the trace of a matrix coefficient
+class TraceCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient * a;
+
+   mutable DenseMatrix ma;
+
+public:
+   /// Construct with the matrix.
+   TraceCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the determinant coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+};
+
 /// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {

From 850f0f7e89e9afa1cdd8f2a31ded9dd403448328 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Fri, 5 Apr 2024 17:45:10 -0400
Subject: [PATCH 10/17] typo

---
 fem/coefficient.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index d6a5743b49..c1a823e3bd 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1820,7 +1820,7 @@ public:
    /// Return the matrix coefficient
    MatrixCoefficient * GetACoef() const { return a; }
 
-   /// Evaluate the determinant coefficient at @a ip.
+   /// Evaluate the trace coefficient at @a ip.
    virtual real_t Eval(ElementTransformation &T,
                        const IntegrationPoint &ip);
 };

From af5a7844a839318a6de865d75d69f96551887242 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Sat, 6 Apr 2024 11:16:53 -0400
Subject: [PATCH 11/17] test for MatrixArrayVectorCoefficient

---
 fem/coefficient.cpp                 |  8 +++-
 fem/coefficient.hpp                 |  7 ++--
 tests/unit/fem/test_coefficient.cpp | 61 +++++++++++++++++++++++++++++
 3 files changed, 72 insertions(+), 4 deletions(-)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index 864973da5d..a38e423d80 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -929,7 +929,7 @@ MatrixArrayVectorCoefficient::MatrixArrayVectorCoefficient (int dim)
 {
    Coeff.SetSize(height);
    ownCoeff.SetSize(height);
-   for (int i = 0; i < (height); i++)
+   for (int i = 0; i < height; i++)
    {
       Coeff[i] = NULL;
       ownCoeff[i] = true;
@@ -947,6 +947,9 @@ void MatrixArrayVectorCoefficient::SetTime(real_t t)
 
 void MatrixArrayVectorCoefficient::Set(int i, VectorCoefficient * c, bool own)
 {
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
    if (ownCoeff[i]) { delete Coeff[i]; }
    Coeff[i] = c;
    ownCoeff[i] = own;
@@ -964,6 +967,9 @@ void MatrixArrayVectorCoefficient::Eval(int i, Vector &V,
                                         ElementTransformation &T,
                                         const IntegrationPoint &ip)
 {
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
    if (Coeff[i])
    {
       Coeff[i] -> Eval(V, T, ip);
diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index c1a823e3bd..4240456cf9 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1334,9 +1334,10 @@ public:
    virtual ~MatrixArrayCoefficient();
 };
 
-/** @brief Matrix coefficient defined by a matrix of vector coefficients.
-    Coefficients that are not set will evaluate to zero in the vector. The
-    coefficient is stored as an Array indexing the rows of the matrix. */
+/** @brief Matrix coefficient defined row-wise by an array of vector
+    coefficients. Rows that are not set will evaluate to zero. The
+    matrix coefficient is stored as an array indexing the rows of
+    the matrix. */
 class MatrixArrayVectorCoefficient : public MatrixCoefficient
 {
 private:
diff --git a/tests/unit/fem/test_coefficient.cpp b/tests/unit/fem/test_coefficient.cpp
index 041a5d0c9a..5128900682 100644
--- a/tests/unit/fem/test_coefficient.cpp
+++ b/tests/unit/fem/test_coefficient.cpp
@@ -303,3 +303,64 @@ TEST_CASE("Piecewise Matrix Coefficient", "[Coefficient]")
       REQUIRE(m.FNorm() == MFEM_Approx(twoNorm));
    }
 }
+
+TEST_CASE("MatrixArrayVectorCoefficient", "[Coefficient]")
+{
+   Vector V1(2), V2(2);
+   V1(0) = 0.0; V1(1) = 1.0;
+   V2(0) = 2.0; V2(1) = 3.0;
+   VectorConstantCoefficient Coef1(V1), Coef2(V2);
+
+   IsoparametricTransformation T;
+   IntegrationPoint ip;
+
+   MatrixArrayVectorCoefficient mavc(2);
+   Vector V(2);
+
+   // Verify zeros for unset rows
+   int row = 0;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   DenseMatrix K(2);
+   mavc.Eval(K, T, ip);
+   REQUIRE(K(0,0) == MFEM_Approx(0.0));
+   REQUIRE(K(0,1) == MFEM_Approx(0.0));
+   REQUIRE(K(1,0) == MFEM_Approx(0.0));
+   REQUIRE(K(1,1) == MFEM_Approx(0.0));
+
+   // Test setting individual rows
+   row = 0;
+   // mavc.Set(row, &Coef1, true); // doesn't work
+   mavc.Set(row, &Coef1, false);
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(1.0));
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   mavc.Set(row, &Coef2, false);
+   row = 0;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(1.0));
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(2.0));
+   REQUIRE(V(1) == MFEM_Approx(3.0));
+
+   mavc.Eval(K, T, ip);
+   REQUIRE(K(0,0) == MFEM_Approx(0.0));
+   REQUIRE(K(0,1) == MFEM_Approx(1.0));
+   REQUIRE(K(1,0) == MFEM_Approx(2.0));
+   REQUIRE(K(1,1) == MFEM_Approx(3.0));
+
+}
\ No newline at end of file

From 821c41fba911ff909fa4ca7b73e19b41d2de8cbb Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Sat, 6 Apr 2024 11:50:53 -0400
Subject: [PATCH 12/17] remove comments

---
 tests/unit/fem/test_coefficient.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/unit/fem/test_coefficient.cpp b/tests/unit/fem/test_coefficient.cpp
index 5128900682..bcfaffe6e9 100644
--- a/tests/unit/fem/test_coefficient.cpp
+++ b/tests/unit/fem/test_coefficient.cpp
@@ -337,7 +337,6 @@ TEST_CASE("MatrixArrayVectorCoefficient", "[Coefficient]")
 
    // Test setting individual rows
    row = 0;
-   // mavc.Set(row, &Coef1, true); // doesn't work
    mavc.Set(row, &Coef1, false);
    mavc.Eval(row, V, T, ip);
    REQUIRE(V(0) == MFEM_Approx(0.0));

From 31a977ac5f36315adaf9071a6149526a0e3da042 Mon Sep 17 00:00:00 2001
From: Brendan Keith <keith10@llnl.gov>
Date: Sun, 7 Apr 2024 22:08:33 -0400
Subject: [PATCH 13/17] Update fem/coefficient.cpp

Co-authored-by: Dohyun Kim <dhkim.cse@gmail.com>
---
 fem/coefficient.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index a38e423d80..bc2777a8f3 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -985,9 +985,9 @@ void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
                                         const IntegrationPoint &ip)
 {
    K.SetSize(height, width);
+   Vector V(width);
    for (int i = 0; i < height; i++)
    {
-      Vector V(width);
       this->Eval(i, V, T, ip);
       for (int j = 0; j < width; j++)
       {

From 0494eb22e6826e9bc80dcd39d815b5b2a0b4199f Mon Sep 17 00:00:00 2001
From: Brendan Keith <keith10@llnl.gov>
Date: Sun, 7 Apr 2024 22:08:55 -0400
Subject: [PATCH 14/17] Update fem/coefficient.hpp

Co-authored-by: Dohyun Kim <dhkim.cse@gmail.com>
---
 fem/coefficient.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 4240456cf9..7a62f235d2 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1368,8 +1368,8 @@ public:
              const IntegrationPoint &ip);
 
    /// Evaluate the matrix coefficient @a ip.
-   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip);
+   void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip) override;
 
    virtual ~MatrixArrayVectorCoefficient();
 };

From 77a3bb103cec74032dfab2134ada0466a4837ee8 Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Sun, 7 Apr 2024 22:12:02 -0400
Subject: [PATCH 15/17] style

---
 fem/coefficient.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 7a62f235d2..0cca620fd8 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1369,7 +1369,7 @@ public:
 
    /// Evaluate the matrix coefficient @a ip.
    void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip) override;
+             const IntegrationPoint &ip) override;
 
    virtual ~MatrixArrayVectorCoefficient();
 };

From 4ebbbc45ae8f4c8ca10cfd934ebd43a0cc80737c Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Sun, 7 Apr 2024 22:20:52 -0400
Subject: [PATCH 16/17] override

---
 fem/coefficient.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 0cca620fd8..6ce9e8c0e6 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1350,7 +1350,7 @@ public:
    explicit MatrixArrayVectorCoefficient (int dim);
 
    /// Set the time for internally stored coefficients
-   void SetTime(real_t t);
+   void SetTime(real_t t) override;
 
    /// Get the coefficient located at (i,j) in the matrix.
    VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }

From 2489c68047548aca88b665ce09cd32db3636bc2d Mon Sep 17 00:00:00 2001
From: Brendan Keith <brendan_keith@brown.edu>
Date: Thu, 23 May 2024 11:40:12 -0400
Subject: [PATCH 17/17] addressing review suggestions

---
 fem/coefficient.cpp | 5 +----
 fem/coefficient.hpp | 2 +-
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
index bc2777a8f3..1c7742c25a 100644
--- a/fem/coefficient.cpp
+++ b/fem/coefficient.cpp
@@ -989,10 +989,7 @@ void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
    for (int i = 0; i < height; i++)
    {
       this->Eval(i, V, T, ip);
-      for (int j = 0; j < width; j++)
-      {
-         K(i,j) = V(j);
-      }
+      K.SetRow(i, V);
    }
 }
 
diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp
index 6ce9e8c0e6..22ecb255f4 100644
--- a/fem/coefficient.hpp
+++ b/fem/coefficient.hpp
@@ -1352,7 +1352,7 @@ public:
    /// Set the time for internally stored coefficients
    void SetTime(real_t t) override;
 
-   /// Get the coefficient located at (i,j) in the matrix.
+   /// Get the vector coefficient located at the i-th row of the matrix
    VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
 
    /** @brief Set the coefficient located at the i-th row of the matrix.