fixed BC issue

Adding coefficients from 2x2_matrix_exponential

examples/CMakeLists.txt

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

examples/ex40.cpp

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

examples/makefile

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

fem/coefficient.hpp

@@ -1352,17 +1352,17 @@ public:
    /// Set the time for internally stored coefficients
    void SetTime(real_t t) override;
 
-   /// Get the coefficient located at (i,j) in the matrix.
+   /// Get the coefficient located at the iᵗʰ row of the matrix.
    VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
 
-   /** @brief Set the coefficient located at the i-th row of the matrix.
+   /** @brief Set the coefficient located at the iᵗʰ row of the matrix.
        By this will take ownership of the Coefficient passed in, but this
        can be overridden with the @a own parameter. */
    void Set(int i, VectorCoefficient * c, bool own=true);
 
    using MatrixCoefficient::Eval;
 
-   /// Evaluate coefficient located at the i-th row of the matrix using integration
+   /// Evaluate coefficient located at the iᵗʰ row of the matrix using integration
    /// point @a ip.
    void Eval(int i, Vector &V, ElementTransformation &T,
              const IntegrationPoint &ip);