Descriptions

MFEM’s ex1 is a basic example of solving a linear elliptic PDE (partial differential equation) problem using the MFEM library. The problem is a simple scalar diffusion equation on a 2D mesh. The equation is defined as:

-\nabla \cdot (-k \nabla u) =f

where $u$ is the solution, $k$ is a scalar coefficient and $f$ is a source term.

The example demonstrates the use of MFEM to set up the finite element discretization of the problem, assemble the linear system, solve the system, and visualize the solution. The example also demonstrates the use of MFEM’s data structures and basic operations such as mesh refinement, finite element space definition, and linear form assembly.

The example uses a 2D triangular mesh of a star shape and solve the problem using a linear finite element space and a constant coefficient. The example uses the PCG iterative solver method with a Gauss-Seidel smoother as a preconditioner to solve the system and shows the result of the solution.

The example is a good starting point for understanding the basic usage of MFEM and how to use it to solve a simple PDE problem. It can also be used as a template to solve other similar problems.

Parse command-line options

   const char *mesh_file = "../data/star.mesh";
   int order = 1;
   bool static_cond = false;
   bool pa = false;
   bool fa = false;
   const char *device_config = "cpu";
   bool visualization = true;
   bool algebraic_ceed = false;

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh",
                  "Mesh file to use.");
   args.AddOption(&order, "-o", "--order",
                  "Finite element order (polynomial degree) or -1 for"
                  " isoparametric space.");
   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                  "--no-static-condensation", "Enable static condensation.");
   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                  "--no-partial-assembly", "Enable Partial Assembly.");
   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
                  "--no-full-assembly", "Enable Full Assembly.");
   args.AddOption(&device_config, "-d", "--device",
                  "Device configuration string, see Device::Configure().");

The above code is a part of the ex1 example program of MFEM (Modular Finite Element Methods). It sets up the command-line options for the program, allowing the user to specify various parameters such as the mesh file to use, the finite element order, and various assembly options.

• The variable “mesh_file” is the path to the mesh file that will be used to solve the linear diffusion problem.
• The variable “order” is the finite element order (polynomial degree) or -1 for isoparametric space.
• The variable “static_cond” is a boolean that enables or disables static condensation.
• The variable “pa” is a boolean that enables or disables partial assembly.
• The variable “fa” is a boolean that enables or disables full assembly.
• The variable “device_config” is the device configuration string, which specifies whether to use the CPU or GPU for computations.
• The variable “visualization” is a boolean that enables or disables visualization of the solution.
• The variable “algebraic_ceed” is a boolean that enables or disables algebraic_ceed.

The OptionsParser class is used to parse the command-line arguments, and the AddOption() method is used to add the different options to the parser. The user can specify the options when running the program.

#ifdef MFEM_USE_CEED
   args.AddOption(&algebraic_ceed, "-a", "--algebraic", "-no-a", "--no-algebraic",
                  "Use algebraic Ceed solver");
#endif

This code is checking if MFEM_USE_CEED is defined and if it is, it adds an additional command-line option, “-a” or “–algebraic” to the options parser. This option allows the user to enable or disable the use of the algebraic CEED solver. CEED (Efficient Linear Algebra for Complex Applications Development) is a library for solving partial differential equations that can accelerate the performance of finite element and boundary element solvers.

The ifdef statement is checking whether the macro MFEM_USE_CEED is defined, if it is defined then the option to use algebraic ceed is added to the command line arguments and if it is not defined then it will not add the option. It means that the option to use algebraic ceed solver will only be available if the library is installed and the macro is defined.

 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
 args.Parse();

This code adds another command-line option, “-vis” or “–visualization” to the options parser, which allows the user to enable or disable GLVis visualization. GLVis (Graphics Library for Visualization) is a lightweight tool for interactively visualizing large-scale finite element discretizations.

The option is added to the OptionsParser using the AddOption() method, which takes a boolean variable, “visualization” as an argument. This variable will be set to true if the user specifies the “-vis” or “–visualization” option on the command line, and false if the user specifies the “-no-vis” or “–no-visualization” option.

After adding the options, the Parse() method is called, which parses the command-line arguments and sets the variables accordingly.

   if (!args.Good())
   {
      args.PrintUsage(cout);
      return 1;
   }
   args.PrintOptions(cout);

This code checks if the command-line options passed to the program are valid using the Good() method of the OptionsParser class. If the options are not valid, the program will display the usage information for the options using the PrintUsage() method and then return 1, indicating that an error occurred.

The PrintOptions(cout) method is used to print the current state of the options, that is the values of the variables after parsing the command line arguments. This can be useful for debugging and understanding what options the user has specified.

The first function will print the usage of the command line options in case the arguments passed are not correct and the program will exit with error code 1.

The second function will print the current state of the options, which can be helpful for user to understand what options they passed and its corresponding values.

Enable hardware devices based on command line options.

Enable hardware devices such as GPUs, and programming models such as CUDA, OCCA, RAJA and OpenMP based on command line options.

This code creates a Device object, which is used to configure the device that will be used for computations. The device can be either the CPU or GPU, depending on the device configuration string passed as an argument to the constructor. The Device class is provided by MFEM library and it is used for the device management.

The Device class constructor takes a string argument, “device_config”, which specifies the device to use. The string can be either “cpu” or “gpu” or any other string depending on the configuration of the device.

The Print() method of the Device class is then called, which prints information about the device configuration. This can be helpful to understand the current settings of the device and also for debugging purpose.

Read the mesh from the given mesh file.

We can handle triangular, quadrilateral, tetrahedral, hexahedral, surface and volume meshes with the same code.

   Mesh mesh(mesh_file, 1, 1);
   int dim = mesh.Dimension();

This code creates a Mesh object, which is used to represent the finite element discretization of the problem. The Mesh class is provided by the MFEM library, and it is used to read and store the finite element mesh data.

The Mesh class constructor takes three arguments:

• The first argument, “mesh_file”, is the path to the mesh file that will be read to create the Mesh object.
• The second argument, 1, is the refinement level, which determines how many times the mesh will be refined. A value of 1 means that the mesh will not be refined.
• The third argument, 1, is the order of the elements in the mesh.

The mesh object is created by reading the input mesh file and it is a representation of the finite element discretization of the problem.

The dimension of the mesh is then obtained by calling the Dimension() method of the Mesh class, which returns the number of spatial dimensions of the mesh, either 2 or 3. The value is then stored in the variable “dim”. This information is required for the further processing of the problem.

Refine the mesh to increase the resolution.

In this example we do ‘ref_levels’ of uniform refinement. We choose ‘ref_levels’ to be the largest number that gives a final mesh with no more than 50,000 elements.

   {
      int ref_levels =
         (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
      for (int l = 0; l < ref_levels; l++)
      {
         mesh.UniformRefinement();
      }
   }

This block of code refines the mesh a certain number of times to increase the resolution of the finite element discretization. The number of refinements is determined by the target number of elements in the mesh, 50000, and the initial number of elements in the mesh obtained from the mesh.GetNE() method.

The refinement level is calculated as the floor of the log base 2 of the ratio of the target number of elements to the initial number of elements, divided by the dimension of the mesh. The floor function is used to round down to the nearest integer value.

The for loop then iteratively calls the UniformRefinement() method of the Mesh class, which refines the mesh uniformly by splitting each element into smaller elements. The loop iterates ref_levels times, where ref_levels is the number of refinements calculated earlier.

This block of code will refine the mesh using the UniformRefinement method until the number of elements reaches the target number of elements, 50000. The refinement is performed only in the spatial dimensions of the mesh. The resulting mesh will have a higher resolution and more accurate solution.

Define a finite element space on the mesh.

Here we use continuous Lagrange finite elements of the specified order. If order < 1, we instead use an isoparametric/isogeometric space.

   FiniteElementCollection *fec;
   bool delete_fec;
   if (order > 0)
   {
      fec = new H1_FECollection(order, dim);
      delete_fec = true;
   }
   else if (mesh.GetNodes())
   {
      fec = mesh.GetNodes()->OwnFEC();
      delete_fec = false;
      cout << "Using isoparametric FEs: " << fec->Name() << endl;
   }
   else
   {
      fec = new H1_FECollection(order = 1, dim);
      delete_fec = true;
   }
   FiniteElementSpace fespace(&mesh, fec);
   cout << "Number of finite element unknowns: "
        << fespace.GetTrueVSize() << endl;

This code creates a FiniteElementCollection object and a FiniteElementSpace object, which are used to define the finite element discretization of the problem.

The FiniteElementCollection object (fec) is created based on the value of the “order” variable, which is passed as an argument to the program. If the value of “order” is greater than 0, a new H1_FECollection object is created with the specified order and dimension, and the variable delete_fec is set to true. H1_FECollection is a collection of conforming, continuous H1-compatible finite elements.

If the value of “order” is less than or equal to 0, the program checks if the mesh object has nodes. If the mesh object has nodes, the program uses the finite element collection of the mesh’s nodes and sets the delete_fec variable to false. If the mesh object does not have nodes, the program creates a new H1_FECollection object with order 1 and dimension, and the variable delete_fec is set to true.

The FiniteElementSpace object (fespace) is then created by passing the mesh object and the finite element collection object to its constructor. The FiniteElementSpace object is used to represent the finite element space of the problem.

The size of the unknowns is obtained by calling the GetTrueVSize() method of the FiniteElementSpace class, which returns the number of degrees of freedom in the finite element space. This information is printed to the console.

In summary, this block of code is initializing the finite element collection and space that will be used to discretize the linear diffusion problem. It uses the order and dimension provided as command line arguments and creates the corresponding finite element collection. It also creates the finite element space and prints the number of unknowns in the space.

Determine the list of true (i.e. conforming) essential boundary dofs.

In this example, the boundary conditions are defined by marking all the boundary attributes from the mesh as essential (Dirichlet) and converting them to a list of true dofs.

   Array<int> ess_tdof_list;
   if (mesh.bdr_attributes.Size())
   {
      Array<int> ess_bdr(mesh.bdr_attributes.Max());
      ess_bdr = 1;
      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
   }

This code creates an Array of integers, “ess_tdof_list”, which stores the essential (or Dirichlet) boundary true degrees of freedom of the finite element space. Essential boundary conditions are conditions that are imposed on the solution at the boundary of the domain. In this case, it is assumed that the boundary attributes are defined in the mesh file, and Dirichlet boundary conditions are associated with attribute 1.

The first line creates an Array “ess_bdr” of integers of size mesh.bdr_attributes.Max(). It is initialized with all 1s. The second line calls the GetEssentialTrueDofs() method of the FiniteElementSpace class, passing the “ess_bdr” array and the “ess_tdof_list” array as arguments. The GetEssentialTrueDofs() method extracts the true degrees of freedom on the specified essential (or Dirichlet) boundary, and stores the result in the “ess_tdof_list” array.

The ess_tdof_list array will be used later to impose the Dirichlet boundary conditions on the linear system. The GetEssentialTrueDofs() method takes as input an array that defines the attribute of each boundary and returns an array with the corresponding dofs.

Set up the linear form b

Set up the linear form b(.) which corresponds to the right-hand side of the FEM linear system, which in this case is (1,phi_i) where phi_i are the basis functions in the finite element fespace.

   LinearForm b(&fespace);
   ConstantCoefficient one(1.0);
   b.AddDomainIntegrator(new DomainLFIntegrator(one));
   b.Assemble();

This code creates a LinearForm object, “b”, which is used to represent the right-hand side (RHS) of the linear system. The LinearForm class is provided by the MFEM library, and it is used to represent linear forms in finite element spaces.

The LinearForm object is created by passing a pointer to the finite element space, “fespace”, to its constructor.

A ConstantCoefficient object, “one”, is created with a value of 1.0. The ConstantCoefficient object is used to represent a constant coefficient function in the finite element space.

A DomainLFIntegrator object is then created by passing the “one” ConstantCoefficient object to its constructor. The DomainLFIntegrator class is used to perform the integration of a linear form over a domain.

The AddDomainIntegrator() method of the LinearForm object is then called, passing the DomainLFIntegrator object as an argument. This adds the integrator to the linear form.

Finally, the Assemble() method of the LinearForm object is called. This assembles the linear form, i.e., it computes the values of the linear form at the degrees of freedom of the finite element space.

In summary, this block of code creates a LinearForm that represents the right-hand side of the linear system, it uses a ConstantCoefficient object with a value of 1.0 and a DomainLFIntegrator to integrate the linear form over the domain, and then it assembles the linear form. It means that it computes the values of the linear form at the degrees of freedom of the finite element space.

Define the solution vector x as a finite element grid function

Define the solution vector x as a finite element grid function corresponding to fespace. Initialize x with initial guess of zero, which satisfies the boundary conditions.

   GridFunction x(&fespace);
   x = 0.0;

This code creates a GridFunction object, “x”, which is used to represent the solution of the linear system. The GridFunction class is provided by the MFEM library and it is used to represent functions defined on a finite element space.

The GridFunction object is created by passing a pointer to the finite element space, “fespace”, to its constructor.

The GridFunction object “x” is then assigned the value of 0.0. This initializes the GridFunction to 0.0 at all degrees of freedom in the finite element space.

The GridFunction object “x” represents the unknown solution of the linear system, which will be computed later. This block of code initializes the GridFunction to 0.0 and it will be used to store the solution later.