Merge branch 'master' into Eikonal

examples/ex41.cpp

@@ -0,0 +1,384 @@
+//                                MFEM Example 41
+//
+// Compile with: make ex41
+//
+// Sample runs: ex41 -o 2
+//              ex41 -o 1 -r 4
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              Eikonal equation,
+//
+//                      |∇u| = 1 in Ω,  u = g on ∂Ω.
+//
+//              This example constructs a fast converging sequence,
+//
+//                      uₖ → u  as k → \infty,
+//
+//              by using in Newton's method to solve the sequence of nonlinear
+//              saddle-point problems
+//
+//               Find ψₖ ∈ H(div,Ω) and uₖ ∈ L²(Ω) such that
+//               ( Zₖ(ψₖ) , τ ) + ( uₖ , ∇⋅τ ) = 0                    ∀ τ ∈ H(div,Ω)
+//               ( ∇⋅ψₖ , v )                 = ( -1 + ∇⋅ψₖ₋₁ , v )   ∀ v ∈ L²(Ω)
+//
+//              where Zₖ(ψ) = ψ / ( 1/αₖ + |ψ|² )^{1/2} and αₖ > 0.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height, true),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t tol = 1e-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 (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(&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();
+   int sdim = mesh.SpaceDimension();
+
+  //  MFEM_ASSERT(mesh.bdr_attributes.Size(),
+  //   "This example does not currently support meshes"
+  //   " without boundary attributes."
+  //  )
+
+   bool zero_average = not mesh.bdr_attributes.Size();
+   if (zero_average)
+   {
+      cout << "\nThe domain has no boundary. "
+           << "Solving for zero-average solution.\n" << endl;
+   }
+
+   // 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.
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   L2_FECollection L2fec(order, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec);
+
+   cout << "Number of Hdiv finite element unknowns: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L2 finite element unknowns: "
+        << L2fes.GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e., conforming) essential boundary dofs.
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   if (zero_average)
+   {
+      offsets.Append(1);
+   }
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 6. Define an initial guess for the solution.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 7. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_psi_gf;
+
+   delta_psi_gf.MakeRef(&RTfes,x,offsets[0]);
+   u_gf.MakeRef(&L2fes,x,offsets[1]);
+   delta_psi_gf = 0.0;
+
+   GridFunction psi_old_gf(&RTfes);
+   GridFunction psi_gf(&RTfes);
+   GridFunction u_old_gf(&L2fes);
+   u_old_gf = 0.0;
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   psi_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   SparseMatrix A21(1, L2fes.GetVSize());
+   if (zero_average)
+   {
+      Array<int> rows(L2fes.GetVSize());
+      Vector values(L2fes.GetVSize());
+      for (int i = 0; i < L2fes.GetVSize(); i++)
+      {
+        rows[i] = i;
+      }
+      values = 1.0;
+      A21.AddRow(0,rows,values);
+   }
+   A21.Finalize();
+   SparseMatrix *A12 = Transpose(A21);
+
+   // 10. Iterate
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   for (k = 0; k < max_it; k++)
+   {
+      GridFunction u_tmp(&L2fes);
+      u_tmp = u_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         ConstantCoefficient alpha_cf(alpha);
+
+         LinearForm b0,b1;
+         b0.Update(&RTfes,rhs.GetBlock(0),0);
+         b1.Update(&L2fes,rhs.GetBlock(1),0);
+
+         ZCoefficient Z(sdim, psi_gf, alpha);
+         DZCoefficient DZ(sdim, psi_gf, alpha);
+
+         ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+         b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_Z));
+         b0.Assemble();
+
+         DivergenceGridFunctionCoefficient div_psi_cf(&psi_gf);
+         DivergenceGridFunctionCoefficient div_psi_old_cf(&psi_old_gf);
+         SumCoefficient psi_old_minus_psi(div_psi_old_cf, div_psi_cf, 1.0, -1.0);
+         b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+         b1.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+         b1.Assemble();
+
+         BilinearForm a00(&RTfes);
+         a00.AddDomainIntegrator(new VectorFEMassIntegrator(DZ));
+         ConstantCoefficient eps(1e-1);
+         a00.AddDomainIntegrator(new VectorFEMassIntegrator(eps));
+         a00.Assemble();
+         a00.Finalize();
+         SparseMatrix &A00 = a00.SpMat();
+
+         MixedBilinearForm a10(&RTfes,&L2fes);
+         a10.AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+         a10.Assemble();
+        //  Array<int> dof_elim_array(1);
+        //  int dof_elim = 10;
+        //  dof_elim_array[0] = dof_elim;
+         a10.Finalize();
+         SparseMatrix &A10 = a10.SpMat();
+        //  A10.EliminateRow(dof_elim);
+
+         SparseMatrix *A01 = Transpose(A10);
+
+         BilinearForm a11(&L2fes);
+         ConstantCoefficient neg_eps(-1e-2);
+         a11.AddDomainIntegrator(new MassIntegrator(neg_eps));
+         a11.Assemble(false);
+         a11.Finalize();
+         SparseMatrix &A11 = a11.SpMat();
+
+        //  BilinearForm a11(&L2fes);
+        //  a11.AddDomainIntegrator(new MassIntegrator(zero));
+        //  a11.Assemble(false);
+        //  SparseMatrix A11;
+        //  a11.FormSystemMatrix(dof_elim_array, A11);
+        //  A11.Set(dof_elim, dof_elim, 1.0);
+        //  rhs[offsets[1] + dof_elim] = 0.0;
+
+        //  BlockOperator A(offsets);
+        //  A.SetBlock(0,0,&A00);
+        //  A.SetBlock(1,0,&A10);
+        //  A.SetBlock(0,1,A01);
+
+        //  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);
+
+         BlockMatrix A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,&A11);
+         if (zero_average)
+         {
+            A.SetBlock(1,2,A12);
+            A.SetBlock(2,1,&A21);
+         }
+
+         SparseMatrix * A_mono = A.CreateMonolithic();
+         UMFPackSolver umf(*A_mono);
+         umf.Mult(rhs,x);
+
+         delta_psi_gf.MakeRef(&RTfes, x.GetBlock(0), 0);
+         u_gf.MakeRef(&L2fes, x.GetBlock(1), 0);
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         real_t gamma = 0.9;
+         delta_psi_gf *= gamma;
+         psi_gf += delta_psi_gf;
+
+         if (visualization)
+         {
+            // sol_sock << "solution\n" << mesh << psi_gf << "window_title 'Discrete solution'"
+            sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         delete A01;
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      // alpha *= 2.0;
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                                                const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+    K(i,i) = phi;
+    for (int j = 0; j < height; j++)
+      {
+        K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}

examples/ex41b.cpp

@@ -0,0 +1,341 @@
+//                                MFEM Example 41
+//
+// Compile with: make ex41
+//
+// Sample runs: ex41 -o 2
+//              ex41 -o 2 -r 4
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              Eikonal equation,
+//
+//                      |∇u| = 1 in Ω,  u = g on ∂Ω.
+//
+//              This example constructs a fast converging sequence,
+//
+//                      uₖ → u  as k → \infty,
+//
+//              by using in Newton's method to solve the sequence of nonlinear
+//              saddle-point problems
+//
+//               Find qₖ ∈ H¹₀(Ω) and uₖ ∈ H¹₀(Ω) such that
+//               ( ϕₖ(|∇qₖ|) ∇qₖ , ∇w ) + ( ∇uₖ , ∇w ) = 0           ∀ w ∈ H¹₀(Ω)
+//               ( ∇qₖ , ∇v )        = ( -1 , v ) + ( ∇qₖ₋₁ , ∇v )   ∀ v ∈ H¹₀(Ω)
+//
+//              where ϕₖ(s) = 1 / ( 1/αₖ + s² )^{1/2} and αₖ > 0.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *q;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, GridFunction &q_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), q(&q_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *q;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, GridFunction &q_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height, true),  q(&q_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t tol = 1e-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 (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(&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();
+   int sdim = mesh.SpaceDimension();
+
+   // 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);
+
+   cout << "Number of dofs: "
+        << H1fes.GetTrueVSize() * 2 << endl;
+
+   // 5. Determine the list of true (i.e., conforming) essential boundary dofs.
+   Array<int> ess_tdof_list;
+   if (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      H1fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+   else
+   {
+    ess_tdof_list.Append(0);
+   }
+   
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = H1fes.GetVSize();
+   offsets[2] = H1fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 6. Define an initial guess for the solution.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 7. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_q_gf;
+
+   delta_q_gf.MakeRef(&H1fes,x,offsets[0]);
+   u_gf.MakeRef(&H1fes,x,offsets[1]);
+   delta_q_gf = 0.0;
+
+   GridFunction q_old_gf(&H1fes);
+   GridFunction q_gf(&H1fes);
+   GridFunction u_old_gf(&H1fes);
+   q_old_gf = 0.0;
+   u_old_gf = 0.0;
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   q_gf = 0.0;
+   q_old_gf = q_gf;
+   u_old_gf = u_gf;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 10. Iterate
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   for (k = 0; k < max_it; k++)
+   {
+      GridFunction u_tmp(&H1fes);
+      u_tmp = u_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      ConstantCoefficient alpha_cf(alpha);
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         LinearForm b0,b1;
+         b0.Update(&H1fes,rhs.GetBlock(0),0);
+         b1.Update(&H1fes,rhs.GetBlock(1),0);
+
+         ZCoefficient Z(sdim, q_gf, alpha);
+         DZCoefficient DZ(sdim, q_gf, alpha);
+
+         ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+         b0.AddDomainIntegrator(new DomainLFGradIntegrator(neg_Z));
+         b0.Assemble();
+
+         GradientGridFunctionCoefficient grad_q_cf(&q_gf);
+         GradientGridFunctionCoefficient grad_q_old_cf(&q_old_gf);
+         VectorSumCoefficient grad_q_old_minus_q(grad_q_old_cf, grad_q_cf, 1.0, -1.0);
+         b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+         b1.AddDomainIntegrator(new DomainLFGradIntegrator(grad_q_old_minus_q));
+         b1.Assemble();
+
+         BilinearForm a00(&H1fes);
+        //  a00.AddDomainIntegrator(new DiffusionIntegrator());
+         a00.AddDomainIntegrator(new DiffusionIntegrator(DZ));
+         a00.Assemble();
+         a00.EliminateVDofs(ess_tdof_list, mfem::Operator::DIAG_ZERO);
+        //  a00.EliminateVDofs(ess_tdof_list,x.GetBlock(0),rhs.GetBlock(0),
+        //                           mfem::Operator::DIAG_ONE);
+         a00.Finalize();
+         SparseMatrix &A00 = a00.SpMat();
+
+         BilinearForm a10(&H1fes);
+         a10.AddDomainIntegrator(new DiffusionIntegrator());
+         a10.Assemble();
+         a10.EliminateVDofs(ess_tdof_list,x.GetBlock(0),rhs.GetBlock(1),
+                                  mfem::Operator::DIAG_ONE);
+         a10.Finalize();
+         SparseMatrix &A10 = a10.SpMat();
+
+         SparseMatrix *A01 = Transpose(A10);
+
+        //  BlockOperator A(offsets);
+        //  A.SetBlock(0,0,&A00);
+        //  A.SetBlock(1,0,&A10);
+        //  A.SetBlock(0,1,A01);
+
+        //  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);
+
+         BlockMatrix A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,0,&A10);
+
+         SparseMatrix * A_mono = A.CreateMonolithic();
+         UMFPackSolver umf(*A_mono);
+         umf.Mult(rhs,x);
+
+         delta_q_gf.MakeRef(&H1fes, x.GetBlock(0), 0);
+         u_gf.MakeRef(&H1fes, x.GetBlock(1), 0);
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         real_t gamma = 1.0;
+         delta_q_gf *= gamma;
+         q_gf += delta_q_gf;
+
+         if (visualization)
+         {
+            sol_sock << "solution\n" << mesh << u_tmp << "window_title 'Discrete solution'"
+            // sol_sock << "solution\n" << mesh << q_gf << "window_title 'Discrete solution'"
+            // sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         delete A01;
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      q_old_gf = q_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= 2.0;
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << H1fes.GetTrueVSize() * 2
+             << endl;
+
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(q != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector gradq(vdim);
+   q->GetGradient(T,gradq);
+   real_t norm = gradq.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   V = gradq;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                                                const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(q != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector gradq(height);
+   q->GetGradient(T,gradq);
+   real_t norm = gradq.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+    K(i,i) = phi;
+    for (int j = 0; j < height; j++)
+      {
+        K(i,j) -= gradq(i) * gradq(j) * pow(phi, 3);
+      }
+   }
+}

examples/ex41c.cpp

@@ -0,0 +1,353 @@
+//                                MFEM Example 41
+//
+// Compile with: make ex41
+//
+// Sample runs: ex41 -o 2
+//              ex41 -o 1 -r 4
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              Eikonal equation,
+//
+//                      |∇u| = 1 in Ω,  u = g on ∂Ω.
+//
+//              This example constructs a fast converging sequence,
+//
+//                      uₖ → u  as k → \infty,
+//
+//              by using in Newton's method to solve the sequence of nonlinear
+//              saddle-point problems
+//
+//               Find ψₖ ∈ L²(Ω)ⁿ and uₖ ∈ H¹₀(Ω) such that
+//               ( Zₖ(ψₖ) , τ ) + ( ∇uₖ , τ ) = 0                    ∀ τ ∈ L²(Ω)ⁿ
+//                ( ψₖ , ∇v )                = ( -1 , v) + ( ψₖ₋₁ , ∇v )   ∀ v ∈ H¹₀(Ω)
+//
+//              where Zₖ(ψ) = ψ / ( 1/αₖ + |ψ|² )^{1/2} and αₖ > 0.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height, true),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t tol = 1e-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 (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(&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();
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(mesh.bdr_attributes.Size(),
+      "This example does not currently support meshes"
+      " without boundary attributes."
+   )
+
+   // 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.
+   L2_FECollection L2fec(order, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec, sdim);
+
+   H1_FECollection H1fec(order, dim);
+   FiniteElementSpace H1fes(&mesh, &H1fec);
+
+   cout << "Number of L2 finite element unknowns: "
+        << L2fes.GetTrueVSize() << endl;
+   cout << "Number of H1 finite element unknowns: "
+        << H1fes.GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e., conforming) essential boundary dofs.
+   Array<int> ess_vdof_list;
+   ess_vdof_list.SetSize(H1fes.GetTrueVSize());
+   if (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      H1fes.GetEssentialVDofs(ess_bdr, ess_vdof_list);
+   }
+   
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = L2fes.GetVSize();
+   offsets[2] = H1fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 6. Define an initial guess for the solution.
+   ConstantCoefficient one(-1.0);
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 7. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_psi_gf;
+
+   delta_psi_gf.MakeRef(&L2fes,x,offsets[0]);
+   u_gf.MakeRef(&H1fes,x,offsets[1]);
+   delta_psi_gf = 0.0;
+
+   GridFunction psi_old_gf(&L2fes);
+   GridFunction psi_gf(&L2fes);
+   GridFunction u_old_gf(&H1fes);
+   u_old_gf = 0.0;
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   psi_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 10. Iterate
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   for (k = 0; k < max_it; k++)
+   {
+      GridFunction u_tmp(&H1fes);
+      u_tmp = u_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         ConstantCoefficient alpha_cf(alpha);
+
+         LinearForm b0,b1;
+         b0.Update(&L2fes,rhs.GetBlock(0),0);
+         b1.Update(&H1fes,rhs.GetBlock(1),0);
+
+         ZCoefficient Z(sdim, psi_gf, alpha);
+         DZCoefficient DZ(sdim, psi_gf, alpha);
+
+         ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+         b0.AddDomainIntegrator(new VectorDomainLFIntegrator(neg_Z));
+         b0.Assemble();
+
+         VectorGridFunctionCoefficient psi_cf(&psi_gf);
+         VectorGridFunctionCoefficient psi_old_cf(&psi_old_gf);
+         VectorSumCoefficient psi_old_minus_psi(psi_old_cf, psi_cf, 1.0, -1.0);
+
+         b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+         b1.AddDomainIntegrator(new DomainLFGradIntegrator(psi_old_minus_psi));
+         b1.Assemble();
+
+         BilinearForm a00(&L2fes);
+         a00.AddDomainIntegrator(new VectorMassIntegrator(DZ));
+         // ConstantCoefficient eps(1e-2);
+         // a00.AddDomainIntegrator(new VectorMassIntegrator(eps));
+         a00.Assemble();
+         a00.Finalize();
+         SparseMatrix &A00 = a00.SpMat();
+
+         MixedBilinearForm a01(&H1fes,&L2fes);
+         a01.AddDomainIntegrator(new GradientIntegrator());
+         a01.Assemble();
+         a01.EliminateEssentialBCFromTrialDofs(ess_vdof_list,x.GetBlock(1),rhs.GetBlock(0));
+         a01.Finalize();
+         SparseMatrix &A01 = a01.SpMat();
+
+         SparseMatrix *A10 = Transpose(A01);
+
+         BilinearForm a11(&H1fes);
+         a11.AddDomainIntegrator(new MassIntegrator(zero));
+         a11.Assemble(false);
+         a11.EliminateEssentialBCFromDofs(ess_vdof_list,x.GetBlock(1),rhs.GetBlock(1));
+         a11.Finalize();
+         SparseMatrix &A11 = a11.SpMat();
+
+        //  BlockOperator A(offsets);
+        //  A.SetBlock(0,0,&A00);
+        //  A.SetBlock(1,0,&A10);
+        //  A.SetBlock(0,1,A01);
+
+        //  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);
+
+         BlockMatrix A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,A10);
+         A.SetBlock(0,1,&A01);
+         A.SetBlock(1,1,&A11);
+
+         SparseMatrix * A_mono = A.CreateMonolithic();
+         UMFPackSolver umf(*A_mono);
+         umf.Mult(rhs,x);
+
+         delta_psi_gf.MakeRef(&L2fes, x.GetBlock(0), 0);
+         u_gf.MakeRef(&H1fes, x.GetBlock(1), 0);
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         real_t gamma = 0.1;
+         delta_psi_gf *= gamma;
+         psi_gf += delta_psi_gf;
+
+         if (visualization)
+         {
+            // sol_sock << "solution\n" << mesh << psi_gf << "window_title 'Discrete solution'"
+            sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         delete A10;
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      // alpha *= 2.0;
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << L2fes.GetTrueVSize() + H1fes.GetTrueVSize()
+             << endl;
+
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                                                const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/alpha + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+    K(i,i) = phi;
+    for (int j = 0; j < height; j++)
+      {
+        K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}

fem/bilininteg.cpp

@@ -768,16 +768,17 @@ void GradientIntegrator::AssembleElementMatrix2(
    ElementTransformation &Trans,  DenseMatrix &elmat)
 {
    dim = test_fe.GetDim();
+   int spaceDim = Trans.GetSpaceDim();
    int trial_dof = trial_fe.GetDof();
    int test_dof = test_fe.GetDof();
    real_t c;
    Vector d_col;
 
    dshape.SetSize(trial_dof, dim);
-   gshape.SetSize(trial_dof, dim);
-   Jadj.SetSize(dim);
+   gshape.SetSize(trial_dof, spaceDim);
+   Jadj.SetSize(dim, spaceDim);
    shape.SetSize(test_dof);
-   elmat.SetSize(dim * test_dof, trial_dof);
+   elmat.SetSize(spaceDim * test_dof, trial_dof);
 
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
                                                             Trans);
@@ -804,7 +805,7 @@ void GradientIntegrator::AssembleElementMatrix2(
       }
       shape *= c;
 
-      for (int d = 0; d < dim; ++d)
+      for (int d = 0; d < spaceDim; ++d)
       {
          gshape.GetColumnReference(d, d_col);
          MultVWt(shape, d_col, elmat_comp);