Assemble the bilinear form and the corresponding linear system

Assemble the bilinear form and the corresponding linear system, applying any necessary transformations such as: eliminating boundary conditions, applying conforming constraints for non-conforming AMR, static condensation, etc.

   if (static_cond) { a->EnableStaticCondensation(); }
   a->Assemble();

   SparseMatrix A;
   Vector B, X;
   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

   cout << "Size of linear system: " << A.Height() << endl;

This code first checks if the “static_cond” variable is true. If it is, it calls the EnableStaticCondensation method on the BilinearForm object “a” to enable static condensation.

It then calls the Assemble method on the BilinearForm object “a” to assemble the linear system of equations.

It then declares SparseMatrix object “A”, Vector objects “B” and “X”.

It then calls the FormLinearSystem method on the BilinearForm object “a”, passing it the essential true dofs list, the GridFunction object “x”, the LinearForm object “b”, and the SparseMatrix object “A”, Vector objects “X” and “B” as arguments. This method forms the linear system of equations that needs to be solved.

It then prints the size of the linear system by calling the Height method on the SparseMatrix object “A” and printing the result.

The FormLinearSystem method is used to form the linear system of equations that needs to be solved. It takes the Bilinear form, the linear form, the essential true dofs list, the GridFunction, the matrix and the vectors representing the right-hand side and the unknowns.

This code forms the linear system of equations to be solved, it checks if static condensation is enabled, assembles the system and forms the linear system of equations by calling the FormLinearSystem method on the Bilinear form object. It also prints the size of the linear system.

solve the system Ax=b

#ifndef MFEM_USE_SUITESPARSE
   // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
   //     solve the system Ax=b with PCG.
   GSSmoother M(A);
   PCG(A, M, B, X, 1, 500, 1e-12, 0.0);
#else
   // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
   UMFPackSolver umf_solver;
   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
   umf_solver.SetOperator(A);
   umf_solver.Mult(B, X);
#endif

This code uses a conditional statement to choose the method of solving the linear system based on whether MFEM was compiled with SuiteSparse or not.

If MFEM was not compiled with SuiteSparse, it uses a simple symmetric Gauss-Seidel preconditioner and PCG (Preconditioned Conjugate Gradient) method to solve the linear system. It creates a GSSmoother object, “M”, by passing it the SparseMatrix object “A”, and then calls the PCG method, passing it the SparseMatrix and GSSmoother objects, the Vector objects “B” and “X”, and the number of iterations, tolerance and other parameters as arguments.

If MFEM was compiled with SuiteSparse, it uses the UMFPACK solver to solve the system. It creates an UMFPackSolver object, “umf_solver”, and sets its control parameters such as ordering method. Then it sets the operator (SparseMatrix A) and it calls the Mult method on the UMFPackSolver object, passing it the Vector objects “B” and “X” as arguments.

The GSSmoother class is a simple symmetric Gauss-Seidel smoother, it’s used as a preconditioner in the PCG method. The PCG method is a iterative method to solve a linear system, it’s preconditioned by the GSSmoother. The UMFPackSolver is a solver for sparse linear systems, it’s based on the UMFPACK library.

This code solves the linear system of the PDE, it uses the PCG method preconditioned by the GSSmoother or the UMFPackSolver, depending on whether MFEM was compiled with SuiteSparse or not.

Recover the solution as a finite element grid function

  a->RecoverFEMSolution(X, *b, x);

This code calls the RecoverFEMSolution method on the BilinearForm object “a”, passing it the Vector object “X” that contains the solution of the linear system, the LinearForm object “b”, and the GridFunction object “x” as arguments.

The RecoverFEMSolution method is used to recover the solution of the PDE from the linear system. It takes the solution of the linear system, represented by the Vector object “X”, and the right-hand side of the PDE, represented by the LinearForm object “b”, and it recovers the solution of the PDE, represented by the GridFunction object “x”.

This code recovers the solution of the PDE from the linear system, it uses the RecoverFEMSolution method which takes the solution of the linear system and the right-hand side of the PDE, and it recovers the solution of the PDE on the finite element space.

Compute and print the L^2 norm of the error

cout << "\n|| E_h - E ||_{L^2} = " << x.ComputeL2Error(E) << '\n' << endl;

Computes and outputs the L2 error of the numerical solution.

This code calls the ComputeL2Error method on the GridFunction object “x”, passing it the VectorFunctionCoefficient object “E” as an argument. This method calculates the L2 norm of the difference between the numerical solution represented by the GridFunction object “x” and the exact solution represented by the VectorFunctionCoefficient object “E”.

It then outputs the error by printing it using the “cout” statement.

The ComputeL2Error method is used to calculate the L2 norm of the difference between the numerical solution represented by the GridFunction object and the exact solution represented by the VectorFunctionCoefficient.

This code calculates the L2 error of the numerical solution by comparing it with the exact solution, it uses the ComputeL2Error method which takes the exact solution and the numerical solution and it returns the L2 error.

Save the refined mesh and the solution.

This output can be viewed later using GLVis: “glvis -m refined.mesh -g sol.gf”.

   {
      ofstream mesh_ofs("refined.mesh");
      mesh_ofs.precision(8);
      mesh->Print(mesh_ofs);
      ofstream sol_ofs("sol.gf");
      sol_ofs.precision(8);
      x.Save(sol_ofs);
   }

This code creates two output file stream objects, “mesh_ofs” and “sol_ofs”, by passing it the file names “refined.mesh” and “sol.gf” respectively.

It then sets the precision of these output file stream objects to 8 digits using the precision method.

It then calls the Print method on the Mesh object, passing it the mesh_ofs as an argument. This method writes the refined mesh to the “refined.mesh” file.

It then calls the Save method on the GridFunction object “x”, passing it the sol_ofs as an argument. This method writes the numerical solution to the “sol.gf” file.

The Print method is used to save the refined mesh to a file. The Save method is used to save the numerical solution to a file.

This code saves the refined mesh and the solution to file, it creates two output file stream objects, sets the precision and it uses the Print method to save the refined mesh and the Save method to save the solution.

Send the solution by socket to a GLVis server

   if (visualization)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock(vishost, visport);
      sol_sock.precision(8);
      sol_sock << "solution\n" << *mesh << x << flush;
   }

This code first checks if the visualization variable is true. If it is, it creates a socketstream object, “sol_sock”, by passing it the host name “localhost” and the port number 19916. It then sets the precision of this socketstream object to 8 digits using the precision method.

It then sends the solution to GLVis by writing the string “solution” followed by the mesh and solution data to the socketstream object “sol_sock” using the “<<” operator. The “flush” command is used to ensure that the data is sent immediately.

GLVis is a lightweight tool for visualizing PDE solutions, it can visualize results from MFEM and other finite element libraries. It uses a client-server architecture, the server runs on the machine that performs the simulation, the client runs on the machine that displays the results.

This code sends the solution to the GLVis server, it creates a socketstream object that connects to the GLVis server and it sends the solution data to it.

Free the used memory

   delete a;
   delete sigma;
   delete muinv;
   delete b;
   delete fespace;
   delete fec;
   delete mesh;

This code first deletes the objects created in the reverse order of their creation. It deletes the BilinearForm object “a” first, followed by the Coefficient objects “sigma” and “muinv”, the LinearForm object “b”, the FiniteElementSpace object “fespace”, the FiniteElementCollection object “fec”, and finally the Mesh object “mesh”.

This is important because the objects that were created later in the program may depend on the objects created earlier, so they should be deleted first.

This code releases the dynamically allocated memory by deleting the objects created during the execution of the program in the correct order, it’s important to do this to prevent memory leaks and ensure the correct operation of the program.

Result