Define a finite element space on the mesh.

Here we use vector finite elements, i.e. dim copies of a scalar finite element space. The vector dimension is specified by the last argument of the FiniteElementSpace constructor. For NURBS meshes, we use the (degree elevated) NURBS space associated with the mesh nodes.

 FiniteElementCollection *fec;
   FiniteElementSpace *fespace;
   if (mesh->NURBSext)
   {
      fec = NULL;
      fespace = mesh->GetNodes()->FESpace();
   }
   else
   {
      fec = new H1_FECollection(order, dim);
      fespace = new FiniteElementSpace(mesh, fec, dim);
   }
   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
        << endl << "Assembling: " << flush;

This code sets up the finite element space on which the PDE will be solved.

It first declares a pointer to a FiniteElementCollection object named “fec” and a pointer to a FiniteElementSpace object named “fespace”.

Then it checks whether the mesh has NURBS extension (mesh->NURBSext) or not.

If the mesh has NURBS extension, it sets “fec” to NULL and assigns the FiniteElementSpace of the mesh nodes to “fespace” by calling mesh->GetNodes()->FESpace().

Otherwise, it creates a new H1_FECollection object, passing in the “order” variable and the “dim” variable as parameters. The H1_FECollection is a collection of elements that are suitable for approximating functions in the Sobolev space H1. The “order” variable represents the polynomial degree of the elements in the collection and the “dim” variable represents the dimension of the elements. It then creates a new FiniteElementSpace object, passing in the mesh object, the “fec” object, and the “dim” variable as parameters. The FiniteElementSpace class is used to represent the space of finite elements on a given mesh.

Finally, it prints the number of finite element unknowns in the “fespace” object by calling the fespace->GetTrueVSize() method and also prints a message “Assembling: ” to indicate that the assembly process is about to begin.

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

In this example, the boundary conditions are defined by marking only boundary attribute 1 from the mesh as essential and converting it to a list of true dofs.

   Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
   ess_bdr = 0;
   ess_bdr[0] = 1;
   fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

This code is setting up the essential (or Dirichlet) boundary conditions for the PDE being solved.

It first declares an Array of integers named “ess_tdof_list” and an Array of integers named “ess_bdr” with maximum size of the boundary attributes in the mesh.

It then sets all the elements of the “ess_bdr” array to 0 and sets the first element of the array to 1.

Finally, it calls the fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list) method. This method takes two arguments, the first is the boundary attributes array and the second is the essential true degrees of freedom array. It identifies the essential (Dirichlet) boundary conditions by finding the true degrees of freedom that correspond to the boundary attributes in the ess_bdr array. It then returns the list of the essential true degrees of freedom in the ess_tdof_list array. These true degrees of freedom correspond to the nodes on the boundary where the solution is specified, and they will be set to the appropriate values during the solution process, while the remaining degrees of freedom will be determined by the PDE and the boundary conditions.

Set up the linear form b(.) which corresponds to the right-hand side of the FEM linear system

In this case, b_i equals the boundary integral of f*phi_i where f represents a “pull down” force on the Neumann part of the boundary and phi_i are the basis functions in the finite element fespace. The force is defined by the VectorArrayCoefficient object f, which is a vector of Coefficient objects. The fact that f is non-zero on boundary attribute 2 is indicated by the use of piece-wise constants coefficient for its last component.

   VectorArrayCoefficient f(dim);
   for (int i = 0; i < dim-1; i++)
   {
      f.Set(i, new ConstantCoefficient(0.0));
   }
   {
      Vector pull_force(mesh->bdr_attributes.Max());
      pull_force = 0.0;
      pull_force(1) = -1.0e-2;
      f.Set(dim-1, new PWConstCoefficient(pull_force));
   }

   LinearForm *b = new LinearForm(fespace);
   b->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));
   cout << "r.h.s. ... " << flush;
   b->Assemble();

This code sets up the right-hand side (RHS) of the PDE being solved.

It first creates a VectorArrayCoefficient object “f” with dimension “dim” .

Then it creates a for loop that iterates dim-1 times and sets the coefficient for each component of the vector to a constant coefficient with the value of 0.0.

Then it creates a Vector object named “pull_force” with size of the maximum boundary attributes in the mesh and sets all of its elements to 0.0 and sets the second element to -1.0e-2.

Then it sets the last component of the “f” vector to a PWConstCoefficient object, with the pull_force vector as the parameter. PWConstCoefficient is a class that takes a vector-valued constant coefficient and creates an array of constant coefficients, one for each component of the vector.

It then creates a new LinearForm object named “b” and passing the fespace object as the parameter. The LinearForm class is used to represent the linear form in the weak formulation of a PDE.

It then adds a boundary integrator to the LinearForm object “b” by calling the b->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f)) method and passing in the “f” vector as the parameter. The VectorBoundaryLFIntegrator is an integrator that computes the boundary contribution to the linear form given a VectorArrayCoefficient

Finally, it assembles the linear form by calling the b->Assemble() method. This method computes the contributions to the linear form from all the elements in the finite element mesh and adds them up to form the complete RHS of the PDE. The message “r.h.s. … ” is printed to indicate that the assembly of the RHS is in progress.

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 new GridFunction object named “x” and passing the fespace object as the parameter. The GridFunction class is used to represent a finite element function on a given finite element space.

It then sets all the values of the GridFunction object “x” to 0.0 by calling x = 0.0. This initializes the solution vector to zero, which is often used as the initial guess for iterative solvers before the actual solution is computed.

Set up the bilinear form a(.,.) on the finite element space

Set up the bilinear form a(.,.) on the finite element space corresponding to the linear elasticity integrator with piece-wise // constants coefficient lambda and mu.

   Vector lambda(mesh->attributes.Max());
   lambda = 1.0;
   lambda(0) = lambda(1)*50;
   PWConstCoefficient lambda_func(lambda);
   Vector mu(mesh->attributes.Max());
   mu = 1.0;
   mu(0) = mu(1)*50;
   PWConstCoefficient mu_func(mu);

   BilinearForm *a = new BilinearForm(fespace);
   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));

This code sets up the bilinear form (matrix) of the PDE being solved.

It first creates a Vector object named “lambda” with the size of the maximum attributes in the mesh and sets all its elements to 1.0. Then it sets the first element to 50 times the second element. This is likely to specify different values of the material property “lambda” for the two regions of the mesh.

It then creates a PWConstCoefficient object named “lambda_func” passing the “lambda” vector as the parameter. PWConstCoefficient is a class that takes a vector-valued constant coefficient and creates an array of constant coefficients, one for each component of the vector.

It then creates another Vector object named “mu” with the size of the maximum attributes in the mesh and sets all its elements to 1.0. Then it sets the first element to 50 times the second element. This is likely to specify different values of the material property “mu” for the two regions of the mesh.

It then creates a PWConstCoefficient object named “mu_func” passing the “mu” vector as the parameter.

It then creates a new BilinearForm object named “a” and passing the fespace object as the parameter. The BilinearForm class is used to represent the bilinear form (matrix) in the weak formulation of a PDE.

It then adds a domain integrator to the BilinearForm object “a” by calling the a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func)) method and passing in the “lambda_func” and “mu_func” as the parameters. The ElasticityIntegrator is an integrator that applies the equations of linear elasticity to the bilinear form. The first argument passed to the constructor is the lambda coefficient, the second argument is the mu coefficient. The ElasticityIntegrator uses these coefficients to compute the matrix entries in the bilinear form corresponding to the linear elasticity equations.

This sets up the bilinear form of the PDE, which will be used to compute the matrix entries in the linear system that needs to be solved to find the solution of the PDE.

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.

   cout << "matrix ... " << flush;
   if (static_cond) { a->EnableStaticCondensation(); }
   a->Assemble();

   SparseMatrix A;
   Vector B, X;
   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
   cout << "done." << endl;

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

This code assembles the linear system that needs to be solved to find the solution of the PDE.

It first calls the a->Assemble() method to assemble the bilinear form.

It then checks if the static_cond variable is true or false. If it is true, it enables static condensation on the bilinear form object by calling a->EnableStaticCondensation(). Static condensation is an optimization technique that can be used when solving a PDE using the finite element method. It reduces the size of the linear system to be solved by eliminating some of the unknowns that correspond to internal degrees of freedom of the elements and it can speed up the solution process.

It then calls the a->FormLinearSystem(ess_tdof_list, x, b, A, X, B) method. This method takes the essential true dofs list, the initial guess, the linear form and returns the sparse matrix A, the solution vector X and the right-hand side vector B. This function forms the linear system, i.e. AX=B, corresponding to the given bilinear and linear forms, and applies the essential boundary conditions specified in the ess_tdof_list.

Finally, it prints the size of the linear system by calling A.Height() method. This tells us the number of unknowns in the linear system that needs to be solved to find the solution of the PDE.

Solve the system Ax=b

Define a simple symmetric Gauss-Seidel preconditioner and use it to solve the system Ax=b with PCG. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system

#ifndef MFEM_USE_SUITESPARSE
   // 11. 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-8, 0.0);
#else
   // 11. 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 solves the linear system A*X=B using either PCG or UMFPACK.

It first checks if MFEM_USE_SUITESPARSE is defined or not.

If it is not defined, it uses a simple symmetric Gauss-Seidel preconditioner (GSSmoother) to solve the system A*X=B with PCG. The GSSmoother class is a simple iterative smoother, commonly used as a preconditioner for Krylov subspace methods. It applies a fixed number of symmetric Gauss-Seidel iterations to the input vector. The PCG method, is a iterative method for solving linear systems that uses a preconditioner, in this case GSSmoother, to speed up the convergence. The method takes the matrix A, the preconditioner M, the right-hand side vector B, the initial guess X, the maximum number of iterations, the tolerance, and the relative residual as input parameters.

Otherwise, if MFEM_USE_SUITESPARSE is defined, it uses UMFPACK to solve the system. UMFPack is a sparse direct solver for linear systems of equations. It uses the UMFPackSolver class and sets the operator to A. Then it calls the Mult(B, X) method on the solver to solve the system and store the solution in the X vector.

This code is solving the linear system A*X=B for the solution X, where A is the coefficient matrix, X is the solution vector and B is the right-hand side vector obtained from the previous step.

Recover the solution as a finite element grid function

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

This code is recovering the solution of the PDE from the solution of the linear system obtained in the previous step.

It calls the a->RecoverFEMSolution(X, *b, x) method. This method takes the solution vector X of the linear system, the linear form and the GridFunction as input and recovers the finite element solution defined on the finite element space from the solution of the linear system. It fills the GridFunction x with the nodal values of the finite element solution.

It is important to note that the solution of the linear system obtained in the previous step, is not the same as the solution of the PDE defined on the finite element space. To obtain the solution of the PDE on the finite element space, the nodal values are recovered by this method.

For non-NURBS meshes, make the mesh curved based on the finite element space.

This means that we define the mesh elements through a fespace based transformation of the reference element. This allows us to save the displaced mesh as a curved mesh when using high-order finite element displacement field. We assume that the initial mesh (read from the file) is not higher order curved mesh compared to the chosen FE space.

  if (!mesh->NURBSext)
   {
      mesh->SetNodalFESpace(fespace);
   }

This code sets the finite element space on the mesh after the solution of the PDE is recovered.

It first checks if the mesh is a NURBS extension or not, by calling mesh->NURBSext.

If it is not a NURBS extension, it sets the finite element space on the mesh by calling the mesh->SetNodalFESpace(fespace) method. This method associates the finite element space with the mesh, so that the solution of the PDE can be visualized using the GLVis visualization tool.

It is important to note that if the mesh is a NURBS extension, it already has a finite element space and this step is not necessary.

Save the displaced mesh and the inverted solution

Save the displaced mesh and the inverted solution which gives the // backward displacements to the original grid). This output can be viewed later using GLVis: “glvis -m displaced.mesh -g sol.gf”.

 {
      GridFunction *nodes = mesh->GetNodes();
      *nodes += x;
      x *= -1;
      ofstream mesh_ofs("displaced.mesh");
      mesh_ofs.precision(8);
      mesh->Print(mesh_ofs);
      ofstream sol_ofs("sol.gf");
      sol_ofs.precision(8);
      x.Save(sol_ofs);
   }

This code saves the displacement of the solution of the PDE to a file.

It first creates a GridFunction pointer named “nodes” that points to the nodes of the mesh by calling mesh->GetNodes().

Then it adds the solution of the PDE to the nodes of the mesh by calling *nodes += x. This will update the coordinates of the nodes of the mesh with the displacement of the solution.

Then it multiplies the solution by -1 by calling x *= -1. This is likely to change the direction of the displacement, so it is visualized correctly in the GLVis visualization tool.

It then creates two ofstream objects named “mesh_ofs” and “sol_ofs” and opens two files named “displaced.mesh” and “sol.gf” respectively.

It then sets the precision of the ofstream objects to 8 by calling mesh_ofs.precision(8) and sol_ofs.precision(8).

It then prints the mesh with the displacement to the “displaced.mesh” file by calling mesh->Print(mesh_ofs) and the solution to the “sol.gf” file by calling x.Save(sol_ofs).

This allows the user to visualize the solution and the displacement of the solution using the GLVis visualization tool.

Send the above data by socket to a GLVis server. Use the “n” and “b” keys in GLVis to visualize the displacements

   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 sends the solution to the GLVis visualization tool to be visualized if the visualization variable is true.

It first checks if the visualization variable is true.

If it is true, it creates a char array named “vishost” with the value “localhost” and an int variable named “visport” with the value 19916. These values are the host and port of the GLVis server.

It then creates a socketstream object named “sol_sock” passing the “vishost” and “visport” as the parameters.

It then sets the precision of the socketstream object to 8 by calling sol_sock.precision(8).

It then sends the solution to the GLVis server by calling sol_sock << “solution\n” << *mesh << x << flush;. The “solution\n” string tells the GLVis server that the following data is a solution, then it sends the mesh and the solution using the << operator. The flush at the end of the line makes sure that the data is sent immediately to the GLVis server.

This allows the user to visualize the solution using the GLVis visualization tool.

Free the used memory

   delete a;
   delete b;
   if (fec)
   {
      delete fespace;
      delete fec;
   }
   delete mesh;

This code cleans up the memory by deleting the objects that were dynamically allocated in the previous steps.

It first deletes the bilinear and linear form objects by calling delete a; and delete b;

Then it checks if the fec variable is not null. If it is not null, it deletes the finite element space and the finite element collection objects by calling delete fespace; and delete fec;

Finally, it deletes the mesh object by calling delete mesh;

This is important to prevent memory leaks and to properly release the memory that was allocated for these objects during the execution of the program.