fixed BC issue

Co-authored-by: Dohyun Kim <dhkim.cse@gmail.com>

examples/CMakeLists.txt

@@ -45,6 +45,7 @@ list(APPEND ALL_EXE_SRCS
   ex37.cpp
   ex38.cpp
   ex39.cpp
+  ex40.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -87,6 +88,7 @@ if (MFEM_USE_MPI)
     ex36p.cpp
     ex37p.cpp
     ex39p.cpp
+    ex40p.cpp
   )
 endif()
 

examples/ex40.cpp

@@ -0,0 +1,411 @@
+//                                MFEM Example 40
+//
+// Compile with: make ex40
+//
+// Sample runs: ex40 -o 2
+//              ex40 -o 2 -r 4
+//
+// Description: This example code demonstrates to how to use MFEM to solve
+//              the Monge–Ampère equation
+//
+//                      det(∇²u) = f in Ω,   u = 0 on ∂Ω.
+//
+//              This example highlights the ExponentialMatrixCoefficient
+//              class, which is used in Newton's method to solve the
+//              variational formulation
+//
+//               Find M ∈ H₀(div,Ω)ⁿ and u ∈ H₀¹(Ω) such that
+//                 (exp(M), N) + (∇u, ∇⋅N) = 0           ∀ N ∈ H₀(div,Ω)ⁿ
+//                 (tr(M), v)              = (ln f, v)   ∀ v ∈ H₀¹(Ω)
+//
+//              where n is the spatial dimension of the domain Ω.
+//
+//
+//              The linearized subproblem is
+//
+//               Find δM ∈ H₀(div,Ω)ⁿ and u ∈ H₀¹(Ω) such that
+//                 (exp(M) δM, N) + (∇u, ∇⋅N) = -(exp(M), N)        ∀ N ∈ H₀(div,Ω)ⁿ
+//                 (tr(δM), v)                = (ln f - tr(M), v)   ∀ v ∈ H₀¹(Ω)
+//
+//
+//              (exp(M) δM, N)    :::   VectorFEMassIntegrator
+//              (∇u, ∇⋅N)         :::   MixedGradDivIntegrator
+//              (tr(δM), v)       :::   MixedDotProductIntegrator
+//              (exp(M), N)       :::   VectorFEDomainLFIntegrator
+//              (ln f - tr(M), v) :::   DomainLFIntegrator
+//
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+real_t exact_solution(const Vector &pt);
+void exact_solution_gradient(const Vector &pt, Vector &grad);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/disc-nurbs.mesh";
+   // const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int max_it = 10;
+   int ref_levels = 1;
+   real_t tol = 1e-5;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   if (dim != 2)
+   {
+      MFEM_ABORT("Example 40 currently only supports 2D problems")
+   }
+
+   // 3. Postprocess the mesh.
+   // 3A. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 3B. Interpolate the geometry after refinement to control geometry error.
+   // NOTE: Minimum second-order interpolation is used to improve the accuracy.
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection H1fec(order, dim);
+   FiniteElementSpace H1fes(&mesh, &H1fec);
+
+   RT_FECollection RTfec(order-1, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   cout << "Number of H¹ degrees of freedom: "
+        << H1fes.GetTrueVSize() << endl;
+   cout << "Number of H(div) degrees of freedom: "
+        << RTfes.GetTrueVSize() * dim << endl;
+
+   Array<int> offsets(4);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = RTfes.GetVSize();
+   offsets[3] = H1fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 5. Determine the list of true (i.e., conforming) essential boundary dofs.
+   Array<int> ess_bdr;
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   // 6. Define constants to be used later.
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+   Vector V1(2), V2(2);
+   V1(0) = 1.0; V1(1) = 0.0;
+   V2(0) = 0.0; V2(1) = 1.0;
+   VectorConstantCoefficient onezero(V1);
+   VectorConstantCoefficient zeroone(V2);
+   ScalarVectorProductCoefficient neg_onezero(-1.0, onezero);
+   ScalarVectorProductCoefficient neg_zeroone(-1.0, zeroone);
+
+   // 7. Define the solution vectors as finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction delta_M1_gf, delta_M2_gf, delta_u_gf;
+
+   delta_M1_gf.MakeRef(&RTfes,x,offsets[0]);
+   delta_M2_gf.MakeRef(&RTfes,x,offsets[1]);
+   delta_u_gf.MakeRef(&H1fes,x,offsets[2]);
+
+   GridFunction M1_gf(&RTfes);
+   GridFunction M2_gf(&RTfes);
+   GridFunction u_gf(&H1fes);
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   FunctionCoefficient exact_coef(exact_solution);
+   VectorFunctionCoefficient exact_grad_coef(dim,exact_solution_gradient);
+   ConstantCoefficient ln_rhs_coef(0.0);
+   u_gf.ProjectCoefficient(exact_coef);
+   // u_gf.ProjectCoefficient(zero);
+   M1_gf = 0.0;
+   M2_gf = 0.0;
+
+   delta_M1_gf = 0.0;
+   delta_M2_gf = 0.0;
+   delta_u_gf = 0.0;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 10. Iterate
+   int k;
+   for (k = 0; k < max_it; k++)
+   {
+      mfem::out << "\nITERATION " << k+1 << endl;
+
+      LinearForm b0,b1,b2;
+      b0.Update(&RTfes,rhs.GetBlock(0),0);
+      b1.Update(&RTfes,rhs.GetBlock(1),0);
+      b2.Update(&H1fes,rhs.GetBlock(2),0);
+
+      VectorGridFunctionCoefficient M1(&M1_gf);
+      VectorGridFunctionCoefficient M2(&M2_gf);
+
+      MatrixArrayVectorCoefficient M(dim);
+      M.Set(0, &M1, false);
+      M.Set(1, &M2, false);
+      ExponentialMatrixCoefficient exp_M(M);
+
+      MatrixVectorProductCoefficient exp_M1(exp_M, onezero);
+      MatrixVectorProductCoefficient exp_M2(exp_M, zeroone);
+      InnerProductCoefficient exp_M11(exp_M1, onezero);
+      InnerProductCoefficient exp_M12(exp_M1, zeroone);
+      InnerProductCoefficient exp_M21(exp_M2, onezero);
+      InnerProductCoefficient exp_M22(exp_M2, zeroone);
+
+      GradientGridFunctionCoefficient grad_u(&u_gf);
+      InnerProductCoefficient neg_dudx(neg_onezero, grad_u);
+      ScalarVectorProductCoefficient neg_exp_M1(-1.0, exp_M1);
+      b0.AddDomainIntegrator(new VectorFEDomainLFDivIntegrator(neg_dudx));
+      b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_exp_M1));
+      b0.Assemble();
+
+      InnerProductCoefficient neg_dudy(neg_zeroone, grad_u);
+      b1.AddDomainIntegrator(new VectorFEDomainLFDivIntegrator(neg_dudy));
+      ScalarVectorProductCoefficient neg_exp_M2(-1.0, exp_M2);
+      b1.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_exp_M2));
+      b1.Assemble();
+
+      InnerProductCoefficient M11(M1, onezero);
+      InnerProductCoefficient M22(M2, zeroone);
+      SumCoefficient trace_M(M11, M22);
+      SumCoefficient rhs2(ln_rhs_coef, trace_M, 1.0, -1.0);
+      b2.AddDomainIntegrator(new DomainLFIntegrator(rhs2));
+      b2.Assemble();
+
+      cout << "b0.Norml2() = " << b0.Norml2() << endl;
+      cout << "b1.Norml2() = " << b1.Norml2() << endl;
+      cout << "b2.Norml2() = " << b2.Norml2() << endl;
+
+      BilinearForm a00(&RTfes);
+      a00.AddDomainIntegrator(new VectorFEMassIntegrator());
+      // a00.AddDomainIntegrator(new VectorFEMassIntegrator(exp_M11));
+      a00.Assemble();
+      a00.EliminateEssentialBC(ess_bdr,x.GetBlock(0),rhs.GetBlock(0),mfem::Operator::DIAG_ONE);
+      a00.Finalize();
+      SparseMatrix &A00 = a00.SpMat();
+
+      BilinearForm a01(&RTfes);
+      a01.AddDomainIntegrator(new VectorFEMassIntegrator(zero));
+      // a01.AddDomainIntegrator(new VectorFEMassIntegrator(exp_M12));
+      a01.Assemble();
+      a01.EliminateEssentialBC(ess_bdr,mfem::Operator::DIAG_ZERO);
+      a01.Finalize();
+      SparseMatrix &A01 = a01.SpMat();
+
+      MixedBilinearForm a02(&H1fes,&RTfes);
+      a02.AddDomainIntegrator(new MixedGradDivIntegrator(neg_onezero));
+      a02.Assemble(false);
+      a02.EliminateTrialDofs(ess_bdr,x.GetBlock(2),rhs.GetBlock(0));
+      a02.EliminateTestDofs(ess_bdr);
+      a02.Finalize();
+      SparseMatrix &A02 = a02.SpMat();
+
+      BilinearForm a10(&RTfes);
+      a10.AddDomainIntegrator(new VectorFEMassIntegrator(zero));
+      // a10.AddDomainIntegrator(new VectorFEMassIntegrator(exp_M21));
+      a10.Assemble();
+      a10.EliminateEssentialBC(ess_bdr,mfem::Operator::DIAG_ZERO);
+      a10.Finalize();
+      SparseMatrix &A10 = a10.SpMat();
+
+      BilinearForm a11(&RTfes);
+      a11.AddDomainIntegrator(new VectorFEMassIntegrator());
+      // a11.AddDomainIntegrator(new VectorFEMassIntegrator(exp_M22));
+      a11.Assemble();
+      a11.EliminateEssentialBC(ess_bdr,x.GetBlock(1),rhs.GetBlock(1),mfem::Operator::DIAG_ONE);
+      a11.Finalize();
+      SparseMatrix &A11 = a11.SpMat();
+
+      MixedBilinearForm a12(&H1fes,&RTfes);
+      a12.AddDomainIntegrator(new MixedGradDivIntegrator(neg_zeroone));
+      a12.Assemble(false);
+      a12.EliminateTrialDofs(ess_bdr,x.GetBlock(2),rhs.GetBlock(1));
+      a12.EliminateTestDofs(ess_bdr);
+      a12.Finalize();
+      SparseMatrix &A12 = a12.SpMat();
+
+      MixedBilinearForm a20(&RTfes,&H1fes);
+      a20.AddDomainIntegrator(new MixedDotProductIntegrator(onezero));
+      a20.Assemble();
+      a20.EliminateTrialDofs(ess_bdr,x.GetBlock(0),rhs.GetBlock(2));
+      a20.EliminateTestDofs(ess_bdr);
+      a20.Finalize();
+      SparseMatrix &A20 = a20.SpMat();
+
+      MixedBilinearForm a21(&RTfes,&H1fes);
+      a21.AddDomainIntegrator(new MixedDotProductIntegrator(zeroone));
+      a21.Assemble();
+      a21.EliminateTrialDofs(ess_bdr,x.GetBlock(1),rhs.GetBlock(2));
+      a21.EliminateTestDofs(ess_bdr);
+      a21.Finalize();
+      SparseMatrix &A21 = a21.SpMat();
+
+      BilinearForm a22(&H1fes);
+      // a22.AddDomainIntegrator(new MassIntegrator(neg_one));
+      a22.AddDomainIntegrator(new MassIntegrator(zero));
+      a22.Assemble(false);
+      a22.EliminateEssentialBC(ess_bdr,x.GetBlock(2),rhs.GetBlock(2),mfem::Operator::DIAG_ONE);
+      a22.Finalize();
+      SparseMatrix &A22 = a22.SpMat();
+
+      cout << "b0.Norml2() = " << b0.Norml2() << endl;
+      cout << "b1.Norml2() = " << b1.Norml2() << endl;
+      cout << "b2.Norml2() = " << b2.Norml2() << endl;
+
+      // BlockOperator A(offsets);
+      // A.SetBlock(0,0,&A00);
+      // A.SetBlock(0,1,&A01);
+      // A.SetBlock(0,2,&A02);
+      // A.SetBlock(1,0,&A10);
+      // A.SetBlock(1,1,&A11);
+      // A.SetBlock(1,2,&A12);
+      // A.SetBlock(2,0,&A20);
+      // A.SetBlock(2,1,&A21);
+      // A.SetBlock(2,2,&A22);
+
+      // BlockDiagonalPreconditioner prec(offsets);
+      // prec.SetDiagonalBlock(0,new GSSmoother(A00));
+      // prec.SetDiagonalBlock(1,new GSSmoother(A11));
+      // prec.SetDiagonalBlock(1,new GSSmoother(A22));
+      // prec.owns_blocks = 1;
+
+      // GMRES(A,prec,rhs,x,1,10000,500,1e-12,0.0);
+
+      BlockMatrix A(offsets);
+      A.SetBlock(0,0,&A00);
+      A.SetBlock(0,1,&A01);
+      A.SetBlock(0,2,&A02);
+      A.SetBlock(1,0,&A10);
+      A.SetBlock(1,1,&A11);
+      A.SetBlock(1,2,&A12);
+      A.SetBlock(2,0,&A20);
+      A.SetBlock(2,1,&A21);
+      A.SetBlock(2,2,&A22);
+
+      SparseMatrix * A_mono = A.CreateMonolithic();
+      UMFPackSolver umf(*A_mono);
+      umf.Mult(rhs,x);
+
+      delta_M1_gf.MakeRef(&RTfes, x.GetBlock(0), 0);
+      delta_M2_gf.MakeRef(&RTfes, x.GetBlock(1), 0);
+      delta_u_gf.MakeRef(&H1fes, x.GetBlock(2), 0);
+
+      real_t Newton_update_size = delta_u_gf.ComputeL2Error(zero);
+
+      real_t gamma = 0.3;
+      delta_M1_gf *= gamma;
+      delta_M2_gf *= gamma;
+      delta_u_gf *= gamma;
+      M1_gf += delta_M1_gf;
+      M2_gf += delta_M2_gf;
+      u_gf += delta_u_gf;
+
+      if (visualization)
+      {
+         // sol_sock << "solution\n" << mesh << delta_M1_gf << "window_title 'Discrete solution'"
+         sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                  << flush;
+         mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << Newton_update_size <<
+                   endl;
+      }
+
+      // if (Newton_update_size < tol || k == max_it-1)
+      // {
+      //    break;
+      // }
+
+      real_t H1_error = u_gf.ComputeH1Error(&exact_coef,&exact_grad_coef);
+      real_t L2_error = u_gf.ComputeL2Error(exact_coef);
+      mfem::out << "L2-error  (|| u - uₕᵏ||)       = " << L2_error << endl;
+      // mfem::out << "H1-error  (|| u - uₕᵏ||)       = " << H1_error << endl;
+
+      cin.get();
+
+   }
+
+   mfem::out << "\n Total iterations: " << k+1
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() * 2 + H1fes.GetTrueVSize()
+             << endl;
+
+   // 11. Exact solution.
+   // if (visualization)
+   // {
+   //    socketstream err_sock(vishost, visport);
+   //    err_sock.precision(8);
+
+   //    GridFunction error_gf(&H1fes);
+   //    error_gf.ProjectCoefficient(exact_coef);
+   //    error_gf -= u_gf;
+
+   //    err_sock << "solution\n" << mesh << error_gf << "window_title 'Error'"  <<
+   //             flush;
+   // }
+
+   return 0;
+}
+
+real_t exact_solution(const Vector &pt)
+{
+   real_t x = pt(0), y = pt(1);
+   return (x*x + y*y) / 2.0 - 4.0;
+}
+
+void exact_solution_gradient(const Vector &pt, Vector &grad)
+{
+   real_t x = pt(0), y = pt(1);
+
+   grad(0) = x;
+   grad(1) = y;
+}

examples/makefile

@@ -23,11 +23,11 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
    ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p
+   ex37p ex39p ex40p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
    ex24p ex25p ex26p ex34p ex35p

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
 {

fem/coefficient.cpp

@@ -924,6 +924,78 @@ 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)
+{
+   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;
+}
+
+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)
+{
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
+   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);
+   Vector V(width);
+   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);
+      }
+   }
+}
+
 void MatrixRestrictedCoefficient::SetTime(real_t t)
 {
    if (c) { c->SetTime(t); }
@@ -1041,6 +1113,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),
@@ -1326,6 +1419,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),

fem/coefficient.hpp

@@ -1334,6 +1334,46 @@ public:
    virtual ~MatrixArrayCoefficient();
 };
 
+/** @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:
+   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) override;
+
+   /// Get the coefficient located at the iᵗʰ row of the matrix.
+   VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
+
+   /** @brief Set the coefficient located at the iᵗʰ 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ᵗʰ 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.
+   void Eval(DenseMatrix &K, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+
+   virtual ~MatrixArrayVectorCoefficient();
+};
+
 
 /** @brief Derived matrix coefficient that has the value of the parent matrix
     coefficient where it is active and is zero otherwise. */
@@ -1761,6 +1801,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 trace 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
 {
@@ -2112,7 +2177,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 +2200,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 +2223,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
 {

linalg/densemat.cpp

@@ -532,6 +532,69 @@ 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]);
+         break;
+      }
+      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).
+         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)
+         {
+            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;
+         }
+         break;
+      }
+      case 3:
+      {
+         MFEM_ABORT("3x3 matrices are not currently supported");
+      }
+      default:
+      {
+         MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
+      }
+   }
+}
+
 real_t DenseMatrix::Det() const
 {
    MFEM_ASSERT(Height() == Width() && Height() > 0,

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;

tests/unit/fem/test_coefficient.cpp

@@ -303,3 +303,63 @@ 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, 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));
+
+}

tests/unit/linalg/test_matrix_dense.cpp

@@ -532,6 +532,67 @@ 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};