sync alpha rule manager & cleanup

examples/denoise.cpp

@@ -0,0 +1,593 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_)
+      : VectorCoefficient(vdim), psi(&psi_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_)
+      : MatrixCoefficient(height),  psi(&psi_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+
+bool CheckVectorComponents(const mfem::GridFunction &gf, double limit)
+{
+    const double* data = gf.GetData();
+    const int size = gf.Size();
+
+    for (int i = 0; i < size; ++i)
+    {
+        if (std::abs(data[i]) > limit)
+        {
+            const int vdim = gf.FESpace()->GetVDim();
+            int dof_index = i / vdim;
+            int component_index = i % vdim;
+
+            std::cout << "--> Condition VIOLATED at DOF #" << dof_index
+                      << ", component " << component_index
+                      << ". Value: " << data[i]
+                      << ", Limit: " << limit << std::endl;
+
+            return false;
+        }
+    }
+
+    return true;
+}
+
+class RotationCoefficient : public MatrixCoefficient {
+public:
+    RotationCoefficient() : MatrixCoefficient(2) {}
+
+    virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                      const IntegrationPoint &ip) {
+        M(0,0) =  0;  M(0,1) = -1;  //  [0, -1]
+        M(1,0) =  1;  M(1,1) =  0;  //  [1, 0]
+    }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int order_l2 = 1; 
+   int max_it = 10;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-6;
+
+   int ex = 1; 
+
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   // args.AddOption(&mesh_file, "-m", "--mesh",
+                  // "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order for RT space");
+   args.AddOption(&order_l2, "-o2", "--order", "FEM order for L2 vec space"); 
+   args.AddOption(&ex, "-ex", "--example", "example number"); 
+   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(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   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);
+
+   Mesh mesh = Mesh::MakeCartesian2D(1, 1, Element::Type::TRIANGLE, false);
+   const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
+
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   H1_FECollection h1fec(order_l2, sdim);
+   FiniteElementSpace h1fes(&mesh, &h1fec);
+
+   L2_FECollection L2fec(order_l2, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec, 2);
+
+   Array<int> ess_tdof_list_rt;
+   RTfes.GetBoundaryTrueDofs(ess_tdof_list_rt);
+   
+   Array<int> ess_tdof_list_h1; 
+   h1fes.GetBoundaryTrueDofs(ess_tdof_list_h1); 
+
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+
+   // ess_tdof_list = 1;
+   // if (mesh.bdr_attributes.Size())
+   // {
+      // Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      // ess_bdr = 1;
+
+      // RTfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   // }
+
+   cout << "Number of H(div) dofs: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L² dofs: "
+        << L2fes.GetTrueVSize() << endl;
+   cout << "Number of H1 dofs: " 
+         << h1fes.GetTrueVSize() << endl; 
+
+   Array<int> offsets({0, RTfes.GetVSize(), h1fes.GetVSize(), L2fes.GetVSize(), 1}); 
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   GridFunction p_gf(&RTfes, x.GetBlock(0)), vphi_gf(&h1fes, x.GetBlock(1)), delta_psi_gf(&L2fes, x.GetBlock(2));
+
+   GridFunction psi_old_gf(&L2fes);
+   GridFunction psi_gf(&L2fes);
+   GridFunction p_old_gf(&RTfes);
+
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   p_gf = 0.0;
+   psi_old_gf = psi_gf;
+   p_old_gf = p_gf;
+
+   VectorGridFunctionCoefficient psi_old_cf(&psi_old_gf), psi_cf(&psi_gf), p_vc(&p_gf); 
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock, varphi_sol_sock, true_sock, varphi_true_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+      
+      // sol_sock  << "keys jlA\n";
+       // turn off perspective & light 
+      // sol_sock << "keys cmA\n"; // colorbar + mesh + anti-alias
+   
+      true_sock.open(vishost,visport);
+      true_sock.precision(8);
+
+      varphi_sol_sock.open(vishost, visport); 
+      varphi_sol_sock.precision(8); 
+   
+      varphi_true_sock.open(vishost, visport); 
+      varphi_true_sock.precision(8); 
+   }
+
+   ConstantCoefficient one_cf(1.0); 
+   ConstantCoefficient neg_one(-1.0);
+   VectorConstantCoefficient zero_vec_cf(Vector({0., 0.}));
+   ConstantCoefficient zero_cf(0.0);
+	VectorConstantCoefficient one_vec_cf(Vector({1., 1.})); 
+	VectorConstantCoefficient neg_one_vec_cf(Vector({-1., -1.})); 
+
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+
+   ConstantCoefficient alpha_cf((real_t) alpha); 
+   ProductCoefficient neg_alpha_cf(neg_one, alpha_cf); 
+
+   ZCoefficient Z(sdim, psi_gf);
+   DZCoefficient DZ(sdim, psi_gf);
+   ScalarMatrixProductCoefficient neg_DZ(-1.0, DZ);
+
+   VectorSumCoefficient psi_newton_res(psi_old_cf, psi_cf, 1., -1.);
+
+   LinearForm b0(&RTfes, rhs.GetBlock(0).GetData()), b2(&L2fes, rhs.GetBlock(2).GetData());
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(psi_newton_res));
+
+   VectorFunctionCoefficient f_coeff(2, [ex](const Vector &x, Vector &u) {
+      // NOTE: constant example 
+      // u(0) = 0.5; 
+      // u(1) = 0.0;
+
+      // NOTE: linear example 
+      // u(0) = x(0); 
+      // u(1) = -1.0 * x(1);
+
+      if (ex == 1) {
+         // NOTE: trig example
+         // u(0) = cosh(M_PI*x(0)) * sin(M_PI*x(1));  
+         // u(1) = sinh(M_PI*x(0)) * cos(M_PI*x(1)); 
+
+         // u /= cosh(M_PI);
+         // u(0) = pow(x(1), 2) * (1. - 2./3 * x(1)); 
+         // u(1) = pow(x(0), 2) * (-1. + 2./3 * x(0)); 
+         // u*= -4.;
+
+         u(0) = sin(2 * M_PI * x(0)) * cos(2*M_PI*x(1)); 
+         u(1) = cos(2*M_PI*x(0))*sin(2*M_PI*x(1)); 
+         u *= -.9;
+
+      }
+      else if (ex == 2) { 
+         // NOTE: trig example 2
+         u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+         u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+
+         u *= (1. + 2. * pow(M_PI, 2)); 
+      }
+   });
+
+   // mfem::Coefficient *f_rhs = nullptr; 
+// 
+   // f_rhs = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+            // {
+               // return 1./4. * sin(M_PI*x(0)) * sin (M_PI*x(1)); 
+            // });
+
+   // FunctionCoefficient f_rhs([](const mfem::Vector &x)
+            // {
+               // return 1./4. * cos(M_PI*x(0)) * cos(M_PI*x(1)); 
+            // }); 
+
+   ScalarVectorProductCoefficient alpha_f_cf(alpha_cf, f_coeff); 
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(alpha_f_cf));
+
+   
+   // ProductCoefficient alpha_f_cf(alpha_cf, f_rhs); 
+   // b0.AddDomainIntegrator(new VectorFEDomainLFDivIntegrator(alpha_f_cf)); 
+
+   b2.AddDomainIntegrator(new VectorDomainLFIntegrator(Z));
+
+
+   RotationCoefficient R;
+   ScalarMatrixProductCoefficient neg_R(neg_one, R); 
+   ScalarMatrixProductCoefficient alpha_R(alpha_cf, R);
+
+   // BlockMatrix A(offsets); 
+
+   BilinearForm a11(&h1fes); 
+   a11.AddDomainIntegrator(new DiffusionIntegrator(neg_one)); 
+
+   // ConstantCoefficient eps_cf(-1e-6);
+
+
+   a11.Assemble(false);
+   a11.Finalize(false);
+
+   SparseMatrix &A11 = a11.SpMat();
+
+   BilinearForm a22(&L2fes);
+   a22.AddDomainIntegrator(new VectorMassIntegrator(neg_DZ)); 
+   // a22.AddDomainIntegrator(new VectorMassIntegrator(eps_cf)); 
+
+
+   // Avg 0 condition
+   int dof_h1(h1fes.GetTrueVSize());
+   LinearForm avg0_data(&h1fes);
+   avg0_data.AddDomainIntegrator(new DomainLFIntegrator(one_cf));
+   avg0_data.Assemble();
+   Array<int> avg0_i({0, dof_h1}), avg0_j(dof_h1);
+   std::iota(avg0_j.begin(), avg0_j.end(), 0);
+   SparseMatrix avg0(avg0_i.GetData(), avg0_j.GetData(), avg0_data.GetData(), 1,
+                     dof_h1, false, false, true);
+   auto avg0T = *Transpose(avg0);
+
+   
+   int k;
+   int total_iterations = 0;
+   real_t increment_p = 0.1;
+   GridFunction p_tmp(&RTfes);
+
+   for (k = 0; k < max_it; k++)
+   {
+      p_tmp = p_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b2.Assemble();
+
+         BlockMatrix A(offsets); 
+
+         BilinearForm a00(&RTfes);
+         a00.AddDomainIntegrator(new DivDivIntegrator(alpha_cf)); 
+         a00.SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+
+         a00.Assemble();
+         a00.EliminateEssentialBC(ess_tdof_list_rt, x.GetBlock(0), rhs.GetBlock(0), mfem::Operator::DIAG_ONE); 
+         a00.Finalize(); 
+
+         SparseMatrix &A00 = a00.SpMat(); 
+
+         a22.Assemble(false);
+         a22.Finalize(false);
+
+         SparseMatrix &A22 = a22.SpMat();
+
+         MixedBilinearForm a01(&h1fes, &RTfes); 
+         a01.AddDomainIntegrator(new MixedVectorGradientIntegrator(alpha_R)); 
+
+         a01.Assemble(false); 
+         a01.EliminateTestDofs(ess_tdof_list_rt); 
+         a01.Finalize(false); 
+
+         SparseMatrix &A01 = a01.SpMat();
+         SparseMatrix *A10 = Transpose(A01);
+         *A10 *= 1.0/alpha; 
+
+         // A.SetBlock(0, 1, &A01); 
+         // A.SetBlock(1, 0, A10); 
+
+         MixedBilinearForm a20(&RTfes, &L2fes);
+         a20.AddDomainIntegrator(new VectorFEMassIntegrator());
+
+         a20.Assemble();
+         a20.EliminateTrialDofs(ess_tdof_list_rt, x.GetBlock(0), rhs.GetBlock(0)); 
+
+         a20.Finalize();
+         SparseMatrix &A20 = a20.SpMat();
+         SparseMatrix *A02 = Transpose(A20);
+
+         // NOTE: this does not work because VectorFEMassIntegrator expects the test & trial spaces to be in a specific order :/ 
+         // MixedBilinearForm a02(&L2fes, &RTfes); 
+         // a02.AddDomainIntegrator(new VectorFEMassIntegrator()); 
+         // a02.Assemble(false); 
+         // a02.EliminateTestDofs(ess_bdr); 
+         // a02.Finalize(); 
+         // SparseMatrix &A02 = a02.SpMat(); 
+         // SparseMatrix *A20 = Transpose(A02); 
+
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,A10);
+         A.SetBlock(0,1,&A01);
+         A.SetBlock(1,1,&A11);
+         A.SetBlock(2,0,&A20); 
+         A.SetBlock(0,2,A02); 
+         A.SetBlock(2,2,&A22);
+         A.SetBlock(3, 1, &avg0); 
+         A.SetBlock(1, 3, &avg0T); 
+
+         // TODO: correct the schur preconditioner
+         // Vector A22_diag(a22.Height()); 
+         // A22.GetDiag(A22_diag); 
+         // A22_diag.Reciprocal(); 
+         // SparseMatrix *S = Mult_AtDA(A20, A22_diag); 
+
+         // prec.SetDiagonalBlock(2, new DSmoother(A22));          
+// #ifndef MFEM_USE_SUITESPARSE 
+         BlockDiagonalPreconditioner prec(offsets);
+
+         prec.SetDiagonalBlock(0,new GSSmoother(A00));
+         prec.SetDiagonalBlock(1,new GSSmoother(A11));
+         prec.SetDiagonalBlock(2,new DSmoother(A22));
+
+         prec.owns_blocks = 3;
+
+         GMRES(A,prec,rhs,x,0,10000,500,1e-12,0.0);
+
+         // prec.SetDiagonalBlock(0, new GSSmoother(*S)); 
+// #else
+         // prec.SetDiagonalBlock(0, new UMFPackSolver(*S));
+         // SparseMatrix *A_mono = A.CreateMonolithic(); 
+         // UMFPackSolver umf(*A_mono); 
+         // umf.Mult(rhs, x); 
+// #endif
+
+         // delete S; 
+
+         p_tmp -= p_gf;
+         real_t Newton_update_size = p_tmp.ComputeL2Error(zero_vec_cf);
+         p_tmp = p_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         // a11.Update();
+         // a22.Update(); 
+         b0.Update();
+         b2.Update();
+
+         if (visualization)
+         {
+            GridFunction p_vec(&L2fes); 
+            p_vec.ProjectCoefficient(p_vc); 
+
+            sol_sock << "solution\n" << mesh << p_vec << "window_title 'Discrete solution '" << flush;            
+
+            varphi_sol_sock << "solution\n" << mesh << vphi_gf << "window_title 'Discrete varphi '" << flush; 
+         }
+
+         mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+
+         if (Newton_update_size < increment_p)
+         {
+            break;
+         }
+
+      }
+
+      p_tmp = p_gf;
+      p_tmp -= p_old_gf;
+      increment_p = p_tmp.ComputeL2Error(zero_vec_cf);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_p << endl;
+
+      p_old_gf = p_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_p < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+      alpha_cf.constant = alpha; 
+   }
+   // delete A01;
+   // delete A10; 
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   VectorFunctionCoefficient exact_coeff(2, [ex](const Vector &x, Vector &u) {      
+      // NOTE: constant example 
+      // u(0) = 0.5; 
+      // u(1) = 0.0; 
+      
+      // NOTE: linear example 
+      // u(0) = x(0); 
+      // u(1) = -1.0 * x(1);
+
+      if (ex == 1) { 
+      // NOTE: trig example 
+         // u(0) = cosh(M_PI*x(0)) * sin(M_PI*x(1));  
+         // u(1) = sinh(M_PI*x(0)) * cos(M_PI*x(1)); 
+         // 
+         // u /= cosh(M_PI);
+         // u(0) = - sin(M_PI*x(0))*cos(M_PI*x(1)); 
+         // u(1) = cos(M_PI*x(0)) * sin(M_PI*x(1)); 
+         // u *= M_PI; 
+
+         // u(0) = x(0) * (1. - x(0))*(1 - 2.*x(1)); 
+         // u(1) = x(1) * (1. - x(1))*(2.*x(0) - 1.); 
+         // u*= 4.;
+
+         u(0) = sin(2 * M_PI * x(0)) * cos(2*M_PI*x(1)); 
+         u(1) = -cos(2*M_PI*x(0))*sin(2*M_PI*x(1)); 
+         u *= 0.9;
+      }
+      else if (ex == 2) { 
+         // NOTE trig example 2 
+         u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+         u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+      }
+   });
+
+   FunctionCoefficient varphi_exact_coeff([](const Vector &x) {
+      return 1 / (2. * M_PI) * sin(2*M_PI*x(0)) * sin(2*M_PI*x(1)); 
+   }); 
+   // GridFunctionCoefficient vphi_exact_gf(&varphi_exact_coeff); 
+   if (visualization) { 
+      GridFunction vphi_exact_gf(&h1fes); 
+      vphi_exact_gf.ProjectCoefficient(varphi_exact_coeff); 
+
+      GridFunction exact_vec(&L2fes); 
+      exact_vec.ProjectCoefficient(exact_coeff); 
+
+      true_sock << "solution\n" << mesh <<  exact_vec << "window_title 'True solution '" << flush; 
+      varphi_true_sock << "solution\n" << mesh << vphi_exact_gf << "window_title 'varphi True Solution '" << flush; 
+   }
+
+    if (CheckVectorComponents(p_gf, 1.0))
+    {
+        std::cout << "Result: SUCCESS. All components are within the limit." << std::endl;
+    }
+    else
+    {
+        std::cout << "Result: FAILURE. At least one component is outside the limit." << std::endl;
+    }
+
+   real_t l2_error = p_gf.ComputeL2Error(exact_coeff); 
+
+   cout << "L2 error: " << l2_error << endl; 
+
+   mfem::Coefficient *div_u_exact = nullptr;
+
+   // NOTE: for constant, linear, trig examples, div p = 0  
+    if (ex == 1)
+   {
+      // For ex=1, the divergence is zero.
+      div_u_exact = new mfem::ConstantCoefficient(0.0);
+   }
+   else if (ex == 2) { 
+   // NOTE: for trig example2, div p != 0 
+      div_u_exact = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+      {
+         return -2. * M_PI * sin(M_PI * x(0)) * sin(M_PI * x(1));
+      });
+   }
+
+   real_t hdiv_error = p_gf.ComputeDivError(div_u_exact);  
+
+   cout << "div error: " << hdiv_error << endl; 
+
+   return 0;
+}
+
+// NOTE: 2D ONLY 
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   V.SetSize(2); 
+
+   for (int i = 0; i < psi_vals.Size(); ++i) { V(i) = tanh(psi_vals(i) / 2.); }
+}
+
+// NOTE: 2D ONLY 
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   K.SetSize(2);
+   K = 0.0;
+   for (int i = 0; i < 2; ++i) { K(i, i) = (1. - pow(tanh(psi_vals(i) / 2.), 2)) / 2.; }
+}

examples/dual_L2.cpp

@@ -0,0 +1,443 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_)
+      : VectorCoefficient(vdim), psi(&psi_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_)
+      : MatrixCoefficient(height),  psi(&psi_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+
+bool CheckVectorComponents(const mfem::GridFunction &gf, double limit)
+{
+    const double* data = gf.GetData();
+    const int size = gf.Size();
+
+    for (int i = 0; i < size; ++i)
+    {
+        if (std::abs(data[i]) > limit)
+        {
+            const int vdim = gf.FESpace()->GetVDim();
+            int dof_index = i / vdim;
+            int component_index = i % vdim;
+
+            std::cout << "--> Condition VIOLATED at DOF #" << dof_index
+                      << ", component " << component_index
+                      << ". Value: " << data[i]
+                      << ", Limit: " << limit << std::endl;
+
+            return false;
+        }
+    }
+
+    return true;
+}
+
+int main(int argc, char *argv[])
+{
+   // const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int order_l2 = 1; 
+   int max_it = 10;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-6;
+
+   int ex = 1; 
+
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   // args.AddOption(&mesh_file, "-m", "--mesh",
+                  // "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order for RT space");
+   args.AddOption(&order_l2, "-o2", "--order", "FEM order for L2 vec space"); 
+   args.AddOption(&ex, "-ex", "--example", "example number"); 
+   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(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   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);
+
+   // Mesh mesh = Mesh::MakeCartesian2D(1, 1, Element::Type::TRIANGLE, false);
+   Mesh mesh = Mesh::MakeCartesian2D(1, 1, Element::Type::QUADRILATERAL, false);
+   
+   const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
+
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   L2_FECollection L2fec(order_l2, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec, 2);
+
+   cout << "Number of H(div) dofs: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L² dofs: "
+        << L2fes.GetTrueVSize() << endl;
+
+   Array<int> offsets({0, RTfes.GetVSize(), L2fes.GetVSize()});
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   GridFunction p_gf(&RTfes, x.GetBlock(0)), delta_psi_gf(&L2fes, x.GetBlock(1));
+
+   GridFunction psi_old_gf(&L2fes);
+   GridFunction psi_gf(&L2fes);
+   GridFunction p_old_gf(&RTfes);
+
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   p_gf = 0.0;
+   psi_old_gf = psi_gf;
+   p_old_gf = p_gf;
+
+   VectorGridFunctionCoefficient psi_old_cf(&psi_old_gf), psi_cf(&psi_gf), p_vc(&p_gf); 
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock, true_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+      
+      // sol_sock  << "keys jlA\n";
+       // turn off perspective & light 
+      // sol_sock << "keys " << "cmA" << endl; // colorbar + mesh + anti-alias
+   
+      true_sock.open(vishost,visport);
+      true_sock.precision(8);
+   }
+
+   ConstantCoefficient neg_one(-1.0);
+   VectorConstantCoefficient zero_vec_cf(Vector({0., 0.}));
+   ConstantCoefficient zero_cf(0.0);
+	VectorConstantCoefficient one_vec_cf(Vector({1., 1.})); 
+
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+
+
+   ConstantCoefficient alpha_cf((real_t) alpha); 
+   ProductCoefficient neg_alpha_cf(neg_one, alpha_cf); 
+
+   ZCoefficient Z(sdim, psi_gf);
+   DZCoefficient DZ(sdim, psi_gf);
+   ScalarMatrixProductCoefficient neg_DZ(-1.0, DZ);
+
+   VectorSumCoefficient psi_newton_res(psi_old_cf, psi_cf, 1., -1.);
+
+   LinearForm b0(&RTfes, rhs.GetBlock(0).GetData()), b1(&L2fes, rhs.GetBlock(1).GetData());
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(psi_newton_res));
+
+   VectorFunctionCoefficient f_coeff(2, [ex](const Vector &x, Vector &u) {
+      // NOTE: constant example 
+      // u(0) = 0.5; 
+      // u(1) = 0.0;
+
+      // NOTE: linear example 
+      // u(0) = x(0); 
+      // u(1) = -1.0 * x(1);
+
+      if (ex == 1) {
+         // NOTE: trig example
+         u(0) = cosh(M_PI*x(0)) * sin(M_PI*x(1));  
+         u(1) = sinh(M_PI*x(0)) * cos(M_PI*x(1)); 
+
+         u /= cosh(M_PI);
+      }
+      else if (ex == 2) { 
+         // NOTE: trig example 2
+         u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+         u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+
+         u *= (1. + 2. * pow(M_PI, 2)); 
+      }
+      else if (ex == 3) { 
+         // u(0) = 1.; 
+         // u(1) = -2. * pow(x(1), 3) + 3. * pow(x(1), 2) + 12. * x(1) - 6.; 
+
+         u(0) = x(0); 
+         u(1) = -x(1); 
+      }
+   });
+
+   ScalarVectorProductCoefficient alpha_f_cf(alpha_cf, f_coeff); 
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(alpha_f_cf));
+
+   b1.AddDomainIntegrator(new VectorDomainLFIntegrator(Z));
+
+   BilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new DivDivIntegrator(alpha_cf)); 
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(alpha_cf)); 
+   
+   a00.Assemble();
+   a00.Finalize(); 
+   SparseMatrix &A00 = a00.SpMat(); 
+
+   MixedBilinearForm a10(&RTfes, &L2fes);
+   a10.AddDomainIntegrator(new VectorFEMassIntegrator());
+   a10.Assemble(false);
+   a10.Finalize(false);
+   SparseMatrix &A10 = a10.SpMat();
+   SparseMatrix *A01 = Transpose(A10);
+
+   BilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new VectorMassIntegrator(neg_DZ)); 
+
+   int k;
+   int total_iterations = 0;
+   real_t increment_p = 0.1;
+   GridFunction p_tmp(&RTfes);
+
+   for (k = 0; k < max_it; k++)
+   {
+      p_tmp = p_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b1.Assemble();
+
+         a11.Assemble(false);
+         a11.Finalize(false);
+         SparseMatrix &A11 = a11.SpMat();
+
+         BlockMatrix A(offsets); 
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,&A11);
+
+#ifndef MFEM_USE_SUITESPARSE 
+         BlockDiagonalPreconditioner prec(offsets);
+         prec.SetDiagonalBlock(0,new GSSmoother(A00));
+         prec.SetDiagonalBlock(1,new GSSmoother(A11));
+         prec.owns_blocks = 1;
+
+         GMRES(A,prec,rhs,x,0,10000,500,1e-12,0.0);
+#else
+         SparseMatrix *A_mono = A.CreateMonolithic(); 
+         UMFPackSolver umf(*A_mono); 
+         umf.Mult(rhs, x); 
+#endif
+         
+         p_tmp -= p_gf;
+         real_t Newton_update_size = p_tmp.ComputeL2Error(zero_vec_cf);
+         p_tmp = p_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a11.Update();
+         b0.Update();
+         b1.Update();
+
+         if (visualization)
+         {
+            GridFunction p_vec(&L2fes); 
+            p_vec.ProjectCoefficient(p_vc); 
+
+            sol_sock << "solution\n" << mesh << p_vec << "window_title 'Discrete solution '" << flush;            
+         }
+
+         mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+
+         if (Newton_update_size < increment_p)
+         {
+            break;
+         }
+
+      }
+
+      p_tmp = p_gf;
+      p_tmp -= p_old_gf;
+      increment_p = p_tmp.ComputeL2Error(zero_vec_cf);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_p << endl;
+
+      p_old_gf = p_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_p < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+      alpha_cf.constant = alpha; 
+   }
+   delete A01;
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   VectorFunctionCoefficient exact_coeff(2, [ex](const Vector &x, Vector &u) {      
+      // NOTE: constant example 
+      // u(0) = 0.5; 
+      // u(1) = 0.0; 
+      
+      // NOTE: linear example 
+      // u(0) = x(0); 
+      // u(1) = -1.0 * x(1);
+
+      if (ex == 1) { 
+      // NOTE: trig example 
+         u(0) = cosh(M_PI*x(0)) * sin(M_PI*x(1));  
+         u(1) = sinh(M_PI*x(0)) * cos(M_PI*x(1)); 
+         
+         u /= cosh(M_PI);
+      }
+      else if (ex == 2) { 
+         // NOTE trig example 2 
+         u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+         u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+      }
+      else if (ex == 3) { 
+         // u(0) = 1.; 
+         // u(1)= 3. * pow(x(1), 2) - 2. * pow(x(1), 3); 
+
+         u(0) = x(0); 
+         u(1) = -x(1); 
+      }
+   });
+
+   GridFunction exact_vec(&L2fes); 
+   exact_vec.ProjectCoefficient(exact_coeff); 
+   if (visualization) { 
+      true_sock << "solution\n" << mesh <<  exact_vec << "window_title 'True solution '" << flush; 
+   }
+
+    if (CheckVectorComponents(p_gf, 1.0))
+    {
+        std::cout << "Result: SUCCESS. All components are within the limit." << std::endl;
+    }
+    else
+    {
+        std::cout << "Result: FAILURE. At least one component is outside the limit." << std::endl;
+    }
+
+   real_t l2_error = p_gf.ComputeL2Error(exact_coeff); 
+
+   cout << "L2 error: " << l2_error << endl; 
+
+   mfem::Coefficient *div_u_exact = nullptr;
+
+   // NOTE: for constant, linear, trig examples, div p = 0  
+    if (ex == 1)
+   {
+      // For ex=1, the divergence is zero.
+      div_u_exact = new mfem::ConstantCoefficient(0.0);
+   }
+   else if (ex == 2) { 
+   // NOTE: for trig example2, div p != 0 
+      div_u_exact = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+      {
+         return -2. * M_PI * sin(M_PI * x(0)) * sin(M_PI * x(1));
+      });
+   }
+   else if (ex == 3) { 
+      // div_u_exact = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+         // {
+            // return 6. * (x(1) - pow(x(1), 2)); 
+         // });
+      div_u_exact = new mfem::ConstantCoefficient(0.0);
+   }
+
+   real_t hdiv_error = p_gf.ComputeDivError(div_u_exact);  
+
+   cout << "div error: " << hdiv_error << endl; 
+
+   return 0;
+}
+
+// NOTE: 2D ONLY 
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   V.SetSize(2); 
+
+   for (int i = 0; i < psi_vals.Size(); ++i) { V(i) = tanh(psi_vals(i) / 2.); }
+}
+
+// NOTE: 2D ONLY 
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   K.SetSize(2);
+   K = 0.0;
+   for (int i = 0; i < psi_vals.Size(); ++i) { K(i, i) = (1. - pow(tanh(psi_vals(i) / 2.), 2)) / 2.; }
+}

examples/dual_L2_3d.cpp

@@ -0,0 +1,553 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <functional>
+
+using namespace std;
+using namespace mfem;
+
+real_t factorial(int n) {
+    if (n <= 1) return 1.0;
+    real_t result = 1.0;
+    for (int i = 2; i <= n; ++i) {
+        result *= i;
+    }
+    return result;
+}
+
+class AlphaRuleManager {
+private:
+    int rule_type;
+    real_t C, m_real, q, r;
+    int m_int;
+    real_t alpha_prev;  // For rule 4
+    std::function<real_t(int)> update_func;
+
+public:
+    AlphaRuleManager(int rule, real_t C_val, real_t m_r, real_t q_val, real_t r_val, int m_i)
+        : rule_type(rule), C(C_val), m_real(m_r), q(q_val), r(r_val), m_int(m_i), alpha_prev(0.0) {
+        setupRule();
+    }
+
+    void setupRule() {
+        switch (rule_type) {
+            case 0:
+                break;
+            case 1: // alpha_k = C * k * (k+1) * ... * (k+m)
+                update_func = [this](int k) {
+                    real_t product = 1.0;
+                    for (int i = 0; i <= m_int; ++i) {
+                        product *= (k + 1 + i);
+                    }
+                    return C * product;
+                };
+                break;
+            case 2: // alpha_k = C * pow(m, k-1)
+                update_func = [this](int k) {
+                    return C * pow(m_real, k);
+                };
+                break;
+            case 3: // alpha_{k+1} = C * k * k!
+                update_func = [this](int k) {
+                    return C * (k + 1) * factorial(k + 1);
+                };
+                break;
+            case 4: // alpha_{k+1} = r^{1/(q-1)} * m^{q^k} - alpha_k
+                update_func = [this](int k) {
+                    real_t alpha_new = pow(r, 1.0 / (q - 1.0)) * pow(m_real, pow(q, k + 1)) - alpha_prev;
+                    alpha_prev = alpha_new;
+                    return alpha_new;
+                };
+                break;
+            default:
+                update_func = [](int k) { return 1.0; }; // fallback
+                break;
+        }
+    }
+
+    real_t updateAlpha(int iteration, real_t current_alpha, real_t growth_rate) {
+        if (rule_type == 0) {
+            return current_alpha * max(growth_rate, 1_r);
+        } else {
+            return update_func(iteration);
+        }
+    }
+
+    real_t getInitialAlpha() {
+        if (rule_type == 4) {
+            return pow(r, 1.0 / (q - 1.0));
+        }
+        return 1.0; 
+    }
+};
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_)
+      : VectorCoefficient(vdim), psi(&psi_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_)
+      : MatrixCoefficient(height),  psi(&psi_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+
+bool CheckVectorComponents(const mfem::GridFunction &gf, double limit)
+{
+    const double* data = gf.GetData();
+    const int size = gf.Size();
+
+    for (int i = 0; i < size; ++i)
+    {
+        if (std::abs(data[i]) > limit)
+        {
+            const int vdim = gf.FESpace()->GetVDim();
+            int dof_index = i / vdim;
+            int component_index = i % vdim;
+
+            std::cout << "--> Condition VIOLATED at DOF #" << dof_index
+                      << ", component " << component_index
+                      << ". Value: " << data[i]
+                      << ", Limit: " << limit << std::endl;
+
+            return false;
+        }
+    }
+
+    return true;
+}
+
+int main(int argc, char *argv[])
+{
+   // const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int order_l2 = 1; 
+   int max_it = 10;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-6;
+
+   int ex = 1; 
+   int alpha_rule = 0;  // 0 = original rule, 1-4 = new rules
+   
+   // Parameters for alpha update rules
+   real_t C = 1.0;      // Constant for all rules
+   int m_int = 2;       // Integer parameter for rules 1 and 3
+   real_t m_real = 2.0; // Real parameter for rule 2
+   real_t q = 2.0;      // Parameter for rule 4
+   real_t r = 2.0;      // Parameter for rule 4
+
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   // args.AddOption(&mesh_file, "-m", "--mesh",
+                  // "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order for RT space");
+   args.AddOption(&order_l2, "-o2", "--order", "FEM order for L2 vec space"); 
+   args.AddOption(&ex, "-ex", "--example", "example number"); 
+   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(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   args.AddOption(&alpha_rule, "-ar", "--alpha-rule",
+                  "Alpha update rule (0=original, 1-4=new rules)");
+   args.AddOption(&C, "-C", "--constant",
+                  "Constant C for alpha update rules");
+   args.AddOption(&m_int, "-mi", "--m-int",
+                  "Integer parameter m for alpha rules 1 and 3");
+   args.AddOption(&m_real, "-mr", "--m-real",
+                  "Real parameter m for alpha rule 2");
+   args.AddOption(&q, "-q", "--q-param",
+                  "Parameter q for alpha rule 4");
+   args.AddOption(&r, "-r-param", "--r-param",
+                  "Parameter r for alpha rule 4");
+   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);
+
+   // NOTE: based on step size update rules in PG paper appendix
+   if (alpha_rule < 0 || alpha_rule > 4) {
+      mfem::out << "Error: alpha_rule must be between 0 and 4" << endl;
+      return 1;
+   }
+   // if (alpha_rule == 2 && m_real <= 1.0) {
+      // mfem::out << "Error: for rule 2, m_real must be > 1" << endl;
+      // return 1;
+   // }
+   // if (alpha_rule == 4 && (q <= 1.0 || r <= 1.0)) {
+      // mfem::out << "Error: for rule 4, q and r must be > 1" << endl;
+      // return 1;
+   // }
+
+//    Mesh mesh = Mesh::MakeCartesian2D(1, 1, Element::Type::TRIANGLE, false);
+   
+   Mesh mesh = Mesh::MakeCartesian3D(1, 1, 1, Element::Type::TETRAHEDRON); 
+   const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
+
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   L2_FECollection L2fec(order_l2, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec, 3);
+
+   cout << "Number of H(div) dofs: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L² dofs: "
+        << L2fes.GetTrueVSize() << endl;
+
+   Array<int> offsets({0, RTfes.GetVSize(), L2fes.GetVSize()});
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   GridFunction p_gf(&RTfes, x.GetBlock(0)), delta_psi_gf(&L2fes, x.GetBlock(1));
+
+   GridFunction psi_old_gf(&L2fes);
+   GridFunction psi_gf(&L2fes);
+   GridFunction p_old_gf(&RTfes);
+
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   p_gf = 0.0;
+   psi_old_gf = psi_gf;
+   p_old_gf = p_gf;
+
+   VectorGridFunctionCoefficient psi_old_cf(&psi_old_gf), psi_cf(&psi_gf), p_vc(&p_gf); 
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock, true_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+      
+      // sol_sock  << "keys jlA\n";
+       // turn off perspective & light 
+      // sol_sock << "keys " << "cmA" << endl; // colorbar + mesh + anti-alias
+   
+      true_sock.open(vishost,visport);
+      true_sock.precision(8);
+   }
+
+   ConstantCoefficient neg_one(-1.0);
+   VectorConstantCoefficient zero_vec_cf(Vector({0., 0.}));
+   ConstantCoefficient zero_cf(0.0);
+	VectorConstantCoefficient one_vec_cf(Vector({1., 1.})); 
+
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+
+
+   ConstantCoefficient alpha_cf((real_t) alpha); 
+   ProductCoefficient neg_alpha_cf(neg_one, alpha_cf); 
+
+   ZCoefficient Z(sdim, psi_gf);
+   DZCoefficient DZ(sdim, psi_gf);
+   ScalarMatrixProductCoefficient neg_DZ(-1.0, DZ);
+
+   VectorSumCoefficient psi_newton_res(psi_old_cf, psi_cf, 1., -1.);
+
+   LinearForm b0(&RTfes, rhs.GetBlock(0).GetData()), b1(&L2fes, rhs.GetBlock(1).GetData());
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(psi_newton_res));
+
+   VectorFunctionCoefficient f_coeff(3, [ex](const Vector &x, Vector &u) {
+
+      if (ex == 1) {
+         // NOTE: trig example
+        real_t a = x(0), b = x(1), c = x(2); 
+        u(0) = std::cos(M_PI * a) * std::sin(M_PI * b) * std::sin(M_PI * c);
+        u(1) = std::sin(M_PI * a) * std::cos(M_PI * b) * std::sin(M_PI * c);
+        u(2) = std::sin(M_PI * a) * std::sin(M_PI * b) * std::cos(M_PI * c); 
+
+        u *= (1.0 + 3.0 * M_PI * M_PI);
+    }
+      else if (ex == 2) { 
+         // NOTE: trig example 2
+        //  u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+        //  u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+// 
+        //  u *= (1. + 2. * pow(M_PI, 2)); 
+        u(1) = cos(M_PI * x(1)) * sin (M_PI * x(2)); 
+        u(2) = cos(M_PI * x(2)) * sin (M_PI * x(1));
+
+        u /= 1 + 2. * pow(M_PI, 2); 
+        u(0) = 1.; 
+      }
+   });
+
+   ScalarVectorProductCoefficient alpha_f_cf(alpha_cf, f_coeff); 
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(alpha_f_cf));
+
+   b1.AddDomainIntegrator(new VectorDomainLFIntegrator(Z));
+
+   BilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new DivDivIntegrator(alpha_cf)); 
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(alpha_cf)); 
+   
+   a00.Assemble();
+   a00.Finalize(); 
+   SparseMatrix &A00 = a00.SpMat(); 
+
+   MixedBilinearForm a10(&RTfes, &L2fes);
+   a10.AddDomainIntegrator(new VectorFEMassIntegrator());
+   a10.Assemble(false);
+   a10.Finalize(false);
+   SparseMatrix &A10 = a10.SpMat();
+   SparseMatrix *A01 = Transpose(A10);
+
+   BilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new VectorMassIntegrator(neg_DZ)); 
+
+   int k;
+   int total_iterations = 0;
+   real_t increment_p = 0.1;
+   GridFunction p_tmp(&RTfes);
+
+   AlphaRuleManager alpha_manager(alpha_rule, C, m_real, q, r, m_int);
+   alpha = alpha_manager.getInitialAlpha(); 
+   alpha_cf.constant = alpha; 
+   
+   mfem::out << "Using alpha update rule " << alpha_rule << " with initial alpha = " << alpha << endl;
+
+   for (k = 0; k < max_it; k++)
+   {
+      p_tmp = p_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b1.Assemble();
+
+         a11.Assemble(false);
+         a11.Finalize(false);
+         SparseMatrix &A11 = a11.SpMat();
+
+         BlockMatrix A(offsets); 
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,&A11);
+
+#ifndef MFEM_USE_SUITESPARSE 
+         BlockDiagonalPreconditioner prec(offsets);
+         prec.SetDiagonalBlock(0,new GSSmoother(A00));
+         prec.SetDiagonalBlock(1,new GSSmoother(A11));
+         prec.owns_blocks = 1;
+
+         GMRES(A,prec,rhs,x,0,10000,500,1e-12,0.0);
+#else
+         SparseMatrix *A_mono = A.CreateMonolithic(); 
+         UMFPackSolver umf(*A_mono); 
+         umf.Mult(rhs, x); 
+#endif
+         
+         p_tmp -= p_gf;
+         real_t Newton_update_size = p_tmp.ComputeL2Error(zero_vec_cf);
+         p_tmp = p_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a11.Update();
+         b0.Update();
+         b1.Update();
+
+         if (visualization)
+         {
+            GridFunction p_vec(&L2fes); 
+            p_vec.ProjectCoefficient(p_vc); 
+
+            sol_sock << "solution\n" << mesh << p_vec << "window_title 'Discrete solution '" << flush;            
+         }
+
+         mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+
+         if (Newton_update_size < increment_p)
+         {
+            break;
+         }
+
+      }
+
+      p_tmp = p_gf;
+      p_tmp -= p_old_gf;
+      increment_p = p_tmp.ComputeL2Error(zero_vec_cf);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_p << endl;
+      mfem::out << "Current alpha = " << alpha << " (rule " << alpha_rule << ")" << endl;
+
+      p_old_gf = p_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_p < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      // Update alpha according to the specified rule
+      alpha = alpha_manager.updateAlpha(k, alpha, growth_rate);
+      alpha_cf.constant = alpha; 
+   }
+   delete A01;
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   VectorFunctionCoefficient exact_coeff(3, [ex](const Vector &x, Vector &u) {      
+      // NOTE: constant example 
+      // u(0) = 0.5; 
+      // u(1) = 0.0; 
+      
+      // NOTE: linear example 
+      // u(0) = x(0); 
+      // u(1) = -1.0 * x(1);
+
+      if (ex == 1) { 
+      // NOTE: trig example
+        real_t a = x(0), b = x(1), c = x(2); 
+
+        u(0) = cos(M_PI*a) * sin(M_PI * b) * sin(M_PI*c); 
+        u(1) = sin(M_PI*a)*cos(M_PI*b)*sin(M_PI*c); 
+        u(2) = sin(M_PI*a)*sin(M_PI*b)*cos(M_PI*c);  
+      }
+      else if (ex == 2) { 
+         // NOTE trig example 2 
+        // u(0) = cos(M_PI * x(0)) * sin (M_PI * x(1)); 
+        // u(1) = cos(M_PI * x(1)) * sin (M_PI * x(0));
+        u(0) = 1.; 
+        u(1) = cos(M_PI * x(1)) * sin (M_PI * x(2)); 
+        u(2) = cos(M_PI * x(2)) * sin (M_PI * x(1));
+      }
+   });
+
+   GridFunction exact_vec(&L2fes); 
+   exact_vec.ProjectCoefficient(exact_coeff); 
+   if (visualization) { 
+      true_sock << "solution\n" << mesh <<  exact_vec << "window_title 'True solution '" << flush; 
+   }
+
+    if (CheckVectorComponents(p_gf, 1.0))
+    {
+        std::cout << "Result: SUCCESS. All components are within the limit." << std::endl;
+    }
+    else
+    {
+        std::cout << "Result: FAILURE. At least one component is outside the limit." << std::endl;
+    }
+
+   real_t l2_error = p_gf.ComputeL2Error(exact_coeff); 
+
+   cout << "L2 error: " << l2_error << endl; 
+
+   mfem::Coefficient *div_u_exact = nullptr;
+
+   // NOTE: for constant, linear, trig examples, div p = 0  
+    if (ex == 1)
+   {
+      // For ex=1, the divergence is zero.
+    //   div_u_exact = new mfem::ConstantCoefficient(0.0);
+    div_u_exact = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+      {
+        real_t a = x(0), b = x(1), c = x(2); 
+
+        real_t sin_product = std::sin(M_PI * a) * std::sin(M_PI * b) * std::sin(M_PI * c);
+        real_t scale_factor = -3.0 * M_PI;
+
+        return scale_factor * sin_product;
+      });
+   }
+   else if (ex == 2) { 
+   // NOTE: for trig example2, div p != 0 
+      div_u_exact = new mfem::FunctionCoefficient([](const mfem::Vector &x)
+      {
+         return -2. * M_PI * sin(M_PI * x(2)) * sin(M_PI * x(1));
+      });
+   }
+
+   real_t hdiv_error = p_gf.ComputeDivError(div_u_exact);  
+
+   cout << "div error: " << hdiv_error << endl; 
+
+   return 0;
+}
+
+// NOTE: 2D ONLY 
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(3);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   V.SetSize(3); 
+
+   for (int i = 0; i < psi_vals.Size(); ++i) { V(i) = tanh(psi_vals(i) / 2.); }
+}
+
+// NOTE: 2D ONLY 
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(3);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   K.SetSize(3);
+   K = 0.0;
+   for (int i = 0; i < psi_vals.Size(); ++i) { K(i, i) = (1. - pow(tanh(psi_vals(i) / 2.), 2)) / 2.; }
+}

examples/imre.cpp

@@ -0,0 +1,348 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// NOTE: 2D ONLY 
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(GridFunction &psi_)
+      : VectorCoefficient(2), psi(&psi_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+// NOTE: 2D ONLY  
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(GridFunction &psi_)
+      : MatrixCoefficient(2),  psi(&psi_){ }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+class RotationCoefficient : public MatrixCoefficient {
+public:
+    RotationCoefficient() : MatrixCoefficient(2) {}
+
+    virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                      const IntegrationPoint &ip) {
+        M(0,0) =  0;  M(0,1) = -1;  //  [0, -1]
+        M(1,0) =  1;  M(1,1) =  0;  //  [1, 0]
+    }
+};
+
+int main(int argc, char *argv[])
+{
+   int order = 1;
+   int max_it = 10;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t beta = 1.0; 
+   real_t tol = 1e-5;
+   real_t growth_rate = 2.; 
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   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(&alpha, "-step", "--step",
+                  "Step size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+               "Geometric step size growth rate, alpha = r**k");
+   args.AddOption(&beta, "-reg", "--regularization",
+                  "Image regularization term beta");
+   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);
+
+   // NOTE : pixel-like mesh? 
+   Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::Type::QUADRILATERAL, false, 2., 2.);
+   // NOTE: shift to [-1, 1]x[-1, 1]
+   mesh.Transform([](const Vector &x, Vector &newx){newx=x; newx -= 1.; }); 
+   const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // NOTE: incorrectly set all orders to the same initially
+   RT_FECollection rtfec(order, dim); 
+   H1_FECollection h1fec(order, sdim);
+   L2_FECollection l2fec(order, sdim);
+   FiniteElementSpace rtfes(&mesh, &rtfec); 
+   FiniteElementSpace h1fes(&mesh, &h1fec);
+   FiniteElementSpace l2fes_vec(&mesh, &l2fec, sdim);
+
+   // NOTE: markers for 0 normal trace BC for RT elements; 
+   // this is an essential BC for H(div) 
+   Array<int> ess_tdof_list;
+   if (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      rtfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   int dof_rt(rtfes.GetTrueVSize()), dof_h1(h1fes.GetTrueVSize()), dof_l2(l2fes_vec.GetTrueVSize());
+
+   Array<int> offsets({0, dof_rt, dof_h1, dof_l2, 1});
+   offsets.PartialSum();
+   BlockVector x(offsets), x_old(offsets), x_old_newt(offsets), rhs(offsets);
+   BlockMatrix A(offsets);
+
+   x = 0.;
+
+   GridFunction p(&rtfes, x.GetBlock(0)); 
+   GridFunction vphi(&h1fes, x.GetBlock(1)), 
+                psi(&l2fes_vec, x.GetBlock(2));
+
+   GridFunction p_old(&rtfes, x_old.GetBlock(0)); 
+   GridFunction vphi_old(&h1fes, x_old.GetBlock(1)),
+                psi_old(&l2fes_vec, x_old.GetBlock(2));
+
+   GridFunction p_old_newt(&rtfes, x_old_newt.GetBlock(0)); 
+   GridFunction vphi_old_newt(&h1fes, x_old_newt.GetBlock(1)),
+                psi_old_newt(&l2fes_vec, x_old_newt.GetBlock(2));
+
+   // Discrete coefficients
+   VectorGridFunctionCoefficient psi_old_cf(&psi_old), psi_cf(&psi), p_cf(&p), vphi_cf(&vphi); 
+   
+   // entropy coefficients 
+   VectorSumCoefficient psi_newton_res(psi_old_cf, psi_cf, 1., -1.);
+
+   ZCoefficient Z(psi);
+   DZCoefficient DZ(psi);
+
+   // Other coefficients
+   ConstantCoefficient one_cf(1.), neg_one_cf(-1.);
+   ConstantCoefficient beta_cf(beta); 
+   ConstantCoefficient alpha_cf(alpha); 
+
+   // TODO: meaningful RHS? 
+   VectorConstantCoefficient one_vec_cf(Vector({1., 1.})); 
+   VectorConstantCoefficient alpha_vec_cf(Vector({alpha, alpha})); 
+   ScalarVectorProductCoefficient beta_vec_cf(beta_cf, one_vec_cf); 
+
+   ScalarMatrixProductCoefficient negDZ(neg_one_cf, DZ); 
+
+   // bilinear forms 
+   BilinearForm p_newton(&rtfes), vphi_newtonC(&h1fes), DZform(&l2fes_vec);
+   MixedBilinearForm vphi_newton(&h1fes, &rtfes), psi_newton(&l2fes_vec, &rtfes);
+
+   p_newton.AddDomainIntegrator(new DivDivIntegrator(alpha_cf)); 
+
+   RotationCoefficient R;
+   vphi_newton.AddDomainIntegrator(new MixedVectorGradientIntegrator(R)); 
+
+   vphi_newtonC.AddDomainIntegrator(new DiffusionIntegrator(neg_one_cf)); 
+
+   psi_newton.AddDomainIntegrator(new VectorMassIntegrator()); 
+
+   DZform.AddDomainIntegrator(new VectorMassIntegrator(negDZ)); 
+
+   // apply 0 essential boundary condition to H(div) space 
+   p_newton.SetDiagonalPolicy(mfem::Operator::DIAG_ONE); 
+   p_newton.Assemble(); 
+   p_newton.EliminateEssentialBC(ess_tdof_list, x.GetBlock(0), rhs.GetBlock(0), mfem::Operator::DIAG_ONE); 
+   p_newton.Finalize(true); 
+   A.SetBlock(0, 0, &p_newton.SpMat()); 
+
+   cout << "assembled A blocK" << endl; 
+
+   vphi_newton.Assemble(); 
+   // vphi_newton.EliminateTrialDofs(ess_tdof_list, x.GetBlock(1), rhs.GetBlock(1)); 
+
+   // NOTE: eliminate H(div) essential BC for B^T block test functions <--> trial functions for B block 
+   vphi_newton.EliminateTestDofs(ess_tdof_list); 
+   vphi_newton.Finalize(true); 
+
+   cout << "assembled vphi" << endl; 
+
+   auto vphi_newtonT = *Transpose(vphi_newton.SpMat()); 
+   A.SetBlock(1, 0, &vphi_newtonT); 
+   A.SetBlock(0, 1, &vphi_newton.SpMat()); 
+   
+   vphi_newtonC.Assemble(); 
+   cout << "C block" << endl; 
+   vphi_newtonC.Finalize(true); 
+   cout << "C block final" << endl; 
+
+   A.SetBlock(1, 1, &vphi_newtonC.SpMat()); 
+
+   DZform.Assemble(); 
+   cout << "DZform" << endl; 
+   DZform.Finalize(true); 
+   A.SetBlock(2, 2, &DZform.SpMat()); 
+
+   // average-value 0 constraint
+   LinearForm avg0_data(&h1fes); 
+   avg0_data.AddDomainIntegrator(new DomainLFIntegrator(one_cf)); 
+   avg0_data.Assemble(); 
+   Array<int> avg0_i({0, dof_h1}), avg0_j(dof_h1); 
+   std::iota(avg0_j.begin(), avg0_j.end(), 0); 
+   SparseMatrix avg0(
+      avg0_i.GetData(), avg0_j.GetData(), avg0_data.GetData(), 1, dof_h1, false, false, true
+   );
+   auto avg0T = *Transpose(avg0); 
+
+   A.SetBlock(3, 1, &avg0); 
+   A.SetBlock(1, 3, &avg0T); 
+
+   // linear forms 
+   LinearForm image(&rtfes), prox_res_lf(&l2fes_vec, rhs.GetBlock(0).GetData()), psi_newton_lf(&l2fes_vec, rhs.GetBlock(2).GetData()); 
+   
+   image.AddDomainIntegrator(new VectorDomainLFIntegrator(beta_vec_cf)); // NOTE: currently just beta * 1. 
+
+   prox_res_lf.AddDomainIntegrator(new VectorDomainLFIntegrator(psi_newton_res)); 
+   psi_newton_lf.AddDomainIntegrator(new VectorDomainLFIntegrator(Z));
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   int prox_it = 0;
+   double prox_residual = tol + 1;
+
+   const double newt_tol = 1e-06;
+   double newt_residual = newt_tol + 1;
+
+   const int newt_max_it = 20;
+   int newt_it = 0;
+
+   // prox loop 
+   while (prox_residual > tol && prox_it < max_it) { 
+      prox_it++; 
+      x_old = x; 
+      cout << "prox res lf assemble" << endl; 
+      prox_res_lf.Assemble(); 
+
+      prox_res_lf.Add(-alpha, image);       
+
+      newt_residual = newt_tol + 1; 
+      newt_it = 0; 
+
+      // newton solve for \delta \psi^k 
+      while (newt_residual > newt_tol && newt_it < newt_max_it) { 
+            std::cout << "\tIteration " << newt_it++ << ": ";
+            x_old_newt = x;
+
+            psi_newton.Assemble(false); 
+
+            cout << "psinewt assemble" << endl; 
+            // NOTE: eliminate H(div) essential BC for D^T block test functions <--> trial functions for D block 
+            psi_newton.EliminateTestDofs(ess_tdof_list);
+            cout << "psi newton elim test" << endl; 
+            psi_newton.Finalize(false); 
+
+            auto psi_newtonT = *Transpose(vphi_newton.SpMat()); 
+
+            A.SetBlock(0, 2, &psi_newton.SpMat()); 
+            A.SetBlock(2, 0, &psi_newtonT); 
+
+            // TODO: add multiplication of new stepsize for A & BT blocks 
+
+            DZform.Assemble(false); 
+            DZform.Finalize(false); 
+
+            prox_res_lf.Assemble(); 
+
+            A.SetBlock(2, 2, &DZform.SpMat()); 
+
+            SparseMatrix *A_mono = A.CreateMonolithic(); 
+            UMFPackSolver umf(*A_mono); 
+            umf.Mult(rhs, x); 
+
+            const double newt_residual_psi = psi_old_newt.ComputeL2Error(psi_cf);
+            const double newt_residual_p = p_old_newt.ComputeL2Error(p_cf); 
+            const double newt_residual_vphi = vphi_old_newt.ComputeL2Error(vphi_cf); 
+
+            newt_residual = sqrt(  pow(newt_residual_p, 2)
+                                 + pow(newt_residual_psi, 2) 
+                                 + pow(newt_residual_vphi, 2));
+
+            cout << "Newton iteration residual" << newt_residual << endl; 
+
+            psi_newton.Update(); 
+      }
+      const double prox_residual_psi = psi_old.ComputeL2Error(psi_cf);
+      const double prox_residual_p = p_old.ComputeL2Error(p_cf); 
+      const double prox_residual_vphi = vphi_old.ComputeL2Error(vphi_cf); 
+
+      prox_residual = sqrt(  pow(prox_residual_p, 2)
+                           + pow(prox_residual_psi, 2) 
+                           + pow(prox_residual_vphi, 2));
+
+      cout << "Prox Iteration: " << prox_it << ": " << prox_residual <<
+                " (" << prox_residual_psi << ", " << prox_residual_p << ", " << prox_residual_vphi << ")" << endl;
+      
+      if (visualization)
+      {
+         sol_sock << "solution\n" << mesh << p << flush;
+
+      }
+
+      alpha *= growth_rate; 
+      growth_rate = min(alpha, 1.e09); 
+   }
+
+}
+
+// NOTE: 2D ONLY 
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   for (int i = 0; i < psi_vals.Size(); ++i) { V(i) = tanh(psi_vals(i)); }
+}
+
+// NOTE: 2D ONLY 
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+
+   Vector psi_vals(2);
+   psi->GetVectorValue(T, ip, psi_vals);
+
+   K = 0.;  
+   for (int i = 0; i < 2; ++i) { K(i, i) = (1. - pow(tanh(psi_vals(i)), 2)) / 2.; }
+}

examples/makefile

@@ -23,7 +23,7 @@ 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 ex40
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40 imre dual_L2 dual_L2_3d denoise
 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 \

examples/plot_rates.py

@@ -0,0 +1,253 @@
+import subprocess
+import re
+import matplotlib.pyplot as plt
+import numpy as np
+
+def run_simulation_and_get_errors():
+    executable_path = './dual_L2' 
+    
+    refinement_levels = range(1, 5) 
+
+    l2_errors = []
+    div_errors = []
+    
+    for refs in refinement_levels:
+        try:
+            command = [
+                executable_path, 
+                '--refs', f'{refs}',
+                '-mi', '1000',
+                '-tol', '1e-14', 
+                '-o', '2',
+                '-o2', '1',
+                '-no-vis'
+            ]
+            print(f"Running command: {' '.join(command)}")
+
+            result = subprocess.run(
+                command, 
+                capture_output=True, 
+                text=True, 
+                check=True
+            )
+
+            output = result.stdout
+            
+            l2_match = re.search(r"L2 error: (\d+\.\d+e[+-]\d+|\d+\.\d+)", output)
+            div_match = re.search(r"div error: (\d+\.\d+e[+-]\d+|\d+\.\d+)", output)
+
+            if l2_match and div_match:
+                l2_error = float(l2_match.group(1))
+                div_error = float(div_match.group(1))
+                
+                l2_errors.append(l2_error)
+                div_errors.append(div_error)
+                
+                print(f"  Refinement: {refs} -> L2 Error: {l2_error:.4e}, Div Error: {div_error:.4e}")
+            else:
+                print(f"  ERROR: Could not parse errors for refinement level {refs}.")
+                print("  --- Full Output ---")
+                print(output)
+                print("  -------------------")
+
+        except subprocess.CalledProcessError as e:
+            print(f"FATAL ERROR: The executable returned a non-zero exit code for refinement level {refs}.")
+            print(f"  Return Code: {e.returncode}")
+            print(f"  Stdout: {e.stdout}")
+            print(f"  Stderr: {e.stderr}")
+            return
+        except Exception as e:
+            print(f"An unexpected error occurred: {e}")
+            return
+
+    print("\n--- Simulation Runs Complete ---")
+
+    if not l2_errors or not div_errors:
+        print("No data was collected. Cannot generate plot.")
+        return
+
+    print("Generating convergence plot...")
+    
+    h_values = [1 / (2**r) for r in refinement_levels]
+    # h_values = refinement_levels
+
+    plt.style.use('seaborn-v0_8-dark')
+    fig, ax = plt.subplots(figsize=(10, 8))
+
+    l2_line, = ax.loglog(h_values, l2_errors, 'o-', label='L2 Error', markersize=8, linewidth=2)
+    div_line, = ax.loglog(h_values, div_errors, 'o-', label='Divergence Error', markersize=8, linewidth=2)
+
+    # if l2_errors:
+        # h_squared = [l2_errors[0] * (h / h_values[0])**2 for h in h_values]
+        # ax.loglog(h_values, h_squared, 'k--', label=r'$O(h^2)$')
+    # if div_errors:
+        # h_linear = [div_errors[0] * (h / h_values[0])**2 for h in h_values]
+        # ax.loglog(h_values, h_linear, 'k:', label=r'$O(h)$')
+
+    if len(h_values) > 3:
+        def draw_slope_triangle(h_vals, error_vals, line):
+            line_color = line.get_color()
+
+            observed_slope = (np.log(error_vals[3]) - np.log(error_vals[2])) / (np.log(h_vals[3]) - np.log(h_vals[2]))
+
+            h_pos = [h_vals[2], h_vals[3]] 
+            y_pos = (error_vals[2] + error_vals[3]) * 1.15
+            
+            tri_x = [h_pos[0], h_pos[1], h_pos[1], h_pos[0]]
+            tri_y = [y_pos, y_pos, y_pos * (h_pos[1]/h_pos[0])**observed_slope, y_pos]
+
+            ax.plot(tri_x, tri_y, color=line_color, linestyle='--')
+            
+            ax.fill(tri_x, tri_y, color=line_color, alpha=0.2)
+            
+            ax.text((h_pos[0] + h_pos[1]) / 2, tri_y[0] * 0.9, '1', color=line_color, ha='center', va='top', fontsize=12, fontweight='bold')
+            ax.text(h_pos[1] * 1.1, (tri_y[1] + tri_y[2]) / 2, f'{observed_slope:.1f}', color=line_color, ha='left', va='center', fontsize=12, fontweight='bold')
+
+        draw_slope_triangle(h_values, l2_errors, l2_line)
+        
+        draw_slope_triangle(h_values, div_errors, div_line)
+
+
+    ax.set_xlabel('Mesh size (h)', fontsize=14)
+    ax.set_ylabel('Error', fontsize=14)
+    ax.set_title('Convergence Rates', fontsize=16, fontweight='bold')
+    
+    ax.tick_params(axis='both', which='major', labelsize=12)
+    ax.grid(True, which="both", ls="--", c='0.7')
+    
+    ax.invert_xaxis()
+
+    legend = ax.legend(fontsize=12, frameon=True, facecolor='white', framealpha=0.8)
+    legend.get_frame().set_edgecolor('black')
+
+    plt.tight_layout()
+    plt.show()
+
+
+def draw_slope_triangle(ax, h_vals, error_vals, line):
+    try:
+        observed_slope = (np.log(error_vals[3]) - np.log(error_vals[2])) / (np.log(h_vals[3]) - np.log(h_vals[2]))
+    except (IndexError, ValueError) as e:
+        print(f"  Warning: Could not calculate slope. Not enough data points? Error: {e}")
+        return
+
+    line_color = line.get_color()
+    h_pos = [h_vals[2], h_vals[3]] 
+    y_pos = (error_vals[2] + error_vals[3]) * 1.15
+    
+    tri_x = [h_pos[0], h_pos[1], h_pos[1], h_pos[0]]
+    tri_y = [y_pos, y_pos, y_pos * (h_pos[1]/h_pos[0])**observed_slope, y_pos]
+
+    ax.plot(tri_x, tri_y, color=line_color, linestyle='--')
+    ax.fill(tri_x, tri_y, color=line_color, alpha=0.2)
+    
+    ax.text((h_pos[0] + h_pos[1]) / 2, tri_y[0] * 0.9, '1', color=line_color, ha='center', va='top', fontsize=12, fontweight='bold')
+    ax.text(h_pos[1] * 1.05, (tri_y[1] + tri_y[2]) / 2, f'{observed_slope:.1f}', color=line_color, ha='left', va='center', fontsize=12, fontweight='bold')
+
+
+def run_and_plot_case(case_config):
+    executable_path = './dual_L2' 
+    refinement_levels = range(1, 5)
+
+    l2_errors = []
+    div_errors = []
+    
+    print(f"\n--- Starting Case: {case_config['title']} ---")
+
+    for refs in refinement_levels:
+        try:
+            command = [
+                executable_path, 
+                '--refs',  f'{refs}',
+                '-ex', str(case_config['problem_num']),
+                '-o', str(case_config['rt_order']),
+                '-o2', '1', 
+                '-tol', '1e-14', 
+                '-mi', '1000',
+                '-no-vis'
+            ]
+            print(f"  Running command: {' '.join(command)}")
+
+            result = subprocess.run(command, capture_output=True, text=True, check=True)
+            output = result.stdout
+            
+            l2_match = re.search(r"L2 error: (\d+\.\d+e[+-]\d+|\d+\.\d+)", output)
+            div_match = re.search(r"div error: (\d+\.\d+e[+-]\d+|\d+\.\d+)", output)
+
+            if l2_match and div_match:
+                l2_errors.append(float(l2_match.group(1)))
+                div_errors.append(float(div_match.group(1)))
+            else:
+                print(f"  ERROR: Could not parse errors for refinement level {refs}.")
+        
+        except (FileNotFoundError, subprocess.CalledProcessError, Exception) as e:
+            print(f"  FATAL ERROR during simulation run: {e}")
+            return 
+
+    if not l2_errors or not div_errors:
+        print("  No data collected for this case. Skipping plot.")
+        return
+
+    h_values = [1 / (2**r) for r in refinement_levels]
+    plt.style.use('seaborn-v0_8-dark')
+    fig, ax = plt.subplots(figsize=(8, 6))
+
+    l2_line, = ax.loglog(h_values, l2_errors, 'o-', label='L2 Error', markersize=8, linewidth=2)
+    div_line, = ax.loglog(h_values, div_errors, 'o-', label='Divergence Error', markersize=8, linewidth=2)
+
+    if len(h_values) > 3:
+        draw_slope_triangle(ax, h_values, l2_errors, l2_line)
+        draw_slope_triangle(ax, h_values, div_errors, div_line)
+
+    ax.set_xlabel('Mesh size (h)', fontsize=14)
+    ax.set_ylabel('Error', fontsize=14)
+    ax.set_title(case_config['title'], fontsize=16, fontweight='bold')
+    ax.tick_params(axis='both', which='major', labelsize=12)
+    ax.grid(True, which="both", ls="--", c='0.7')
+    ax.invert_xaxis()
+    legend = ax.legend(fontsize=12, frameon=True, facecolor='white', framealpha=0.8)
+    legend.get_frame().set_edgecolor('black')
+    plt.tight_layout()
+    # plt.show()
+
+    output_filename = case_config['filename']
+    try:
+        plt.savefig(output_filename, dpi=300, bbox_inches='tight')
+        print(f"  Successfully saved plot to '{output_filename}'")
+    except Exception as e:
+        print(f"  Error saving plot: {e}")
+        
+    plt.close(fig) 
+
+if __name__ == '__main__':
+    simulation_cases = [
+        {
+            'problem_num': 1, 
+            'rt_order': 1, 
+            'title': 'Convergence: Trig Example #1 (Div-Free), RT Order 1',
+            'filename': 'trig1_div_free_rt1.png'
+        },
+        {
+            'problem_num': 1, 
+            'rt_order': 2, 
+            'title': 'Convergence: Trig Example #1 (Div-Free), RT Order 2',
+            'filename': 'trig1_div_free_rt2.png'
+        },
+        {
+            'problem_num': 2, 
+            'rt_order': 1, 
+            'title': 'Convergence: Trig Example #2 (Non Div-Free), RT Order 1',
+            'filename': 'trig2_non_div_free_rt1.png'
+        },
+        {
+            'problem_num': 2, 
+            'rt_order': 2, 
+            'title': 'Convergence: Trig Example #2 (Non Div-Free), RT Order 2',
+            'filename': 'trig2_non_div_free_rt2.png'
+        }
+    ]
+
+    for sim in simulation_cases: 
+        run_and_plot_case(sim)
+
+    # run_simulation_and_get_errors()