Description

This example code solves a simple linear elasticity problem describing a multi-material cantilever beam.

Specifically, we approximate the weak form of

-\nabla \cdot  \sigma(\textit{\textbf{u}}) =0

where

\sigma(\textit{\textbf{u}})=\lambda \nabla \cdot \textit{\textbf{u}} + \mu \left(   \nabla \textit{\textbf{u}}+\nabla  \textit{\textbf{u}} ^T\right)

is the stress tensor corresponding to displacement field $\textit{\textbf{u}}$, and lambda and mu // are the material Lame constants. The boundary conditions are $\textit{\textbf{u}}=0$ on the fixed part of the boundary with attribute 1, and

\sigma(\textit{\textbf{u}}) \cdot \textit{\textbf{n}}=\textit{\textbf{f}} 

on the remainder with $\textit{\textbf{f}}$ being a constant pull down // vector on boundary elements with attribute 2, and zero otherwise. The geometry of the domain is assumed to be as follows:

The example demonstrates the use of high-order and NURBS vector finite element spaces with the linear elasticity bilinear form, meshes with curved elements, and the definition of piece-wise constant and vector coefficient objects. Static condensation is also illustrated.

The purpose of ex2 in the MFEM (Modular Finite Element Methods) library is likely to demonstrate the use of the library to solve a simple linear elliptic partial differential equation (PDE) using the finite element method. The specific PDE and boundary conditions will depend on the particular example, but it is likely to be a simple model problem that is easy to understand and can be used to demonstrate the basic functionality of the library. Additionally, it could also serve as a tutorial for the users to understand how to use the library to solve other similar PDE problems.

Parse command-line options

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

int main(int argc, char *argv[])
{
   // 1. Parse command-line options.
   const char *mesh_file = "../data/beam-tri.mesh";
   int order = 1;
   bool static_cond = false;
   bool visualization = 1;

   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(&static_cond, "-sc", "--static-condensation", "-no-sc",
                  "--no-static-condensation", "Enable static condensation.");
   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);

The code you provided is a C++ program that uses the MFEM library to solve a linear elliptic PDE. The program uses the “OptionsParser” class to parse command-line arguments passed to the program. The arguments allow the user to specify the mesh file to use, the finite element order (polynomial degree), whether to enable static condensation and whether to enable GLVis visualization.

The variable “mesh_file” is a string that stores the name of the mesh file to be used, which is set to “../data/beam-tri.mesh” by default.

The variable “order” is an integer that stores the finite element order (polynomial degree) and is set to 1 by default.

The variable “static_cond” is a boolean that determines whether to enable static condensation and is set to false by default.

The variable “visualization” is a boolean that determines whether to enable GLVis visualization and is set to true by default.

The program uses the “OptionsParser” class to parse the command-line arguments and set the values of these variables accordingly. If the user does not provide valid arguments or if there is an error in the arguments, the program will print the usage instructions and exit.

Read the mesh from the given mesh file

Read the mesh from the given mesh file. We can handle triangular, quadrilateral, tetrahedral or hexahedral elements with the same code.

   Mesh *mesh = new Mesh(mesh_file, 1, 1);
   int dim = mesh->Dimension();

   if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
   {
      cerr << "\nInput mesh should have at least two materials and "
           << "two boundary attributes! (See schematic in ex2.cpp)\n"
           << endl;
      return 3;
   }

The code you provided creates a new Mesh object by passing in the name of the mesh file, “mesh_file”, and the values 1, 1 as parameters. The Mesh class is part of the MFEM library and is used to represent the finite element mesh on which the PDE will be solved. The first parameter passed to the constructor is the name of the mesh file, which is specified by the user via the command-line arguments. The two values 1, 1 are additional parameters that are used in this specific example for creating a new Mesh object.

It then defines an integer variable “dim” and assigns the dimension of the mesh to it by calling the mesh->Dimension() method.

It then checks if the maximum attribute in the mesh is less than 2 or the maximum boundary attribute is less than 2, if that is the case it prints an error message saying that the input mesh should have at least 2 materials and 2 boundary attributes and returns 3 as the error code. This check is likely to ensure that the input mesh has the necessary number of materials and boundary attributes required by the specific example.

Select the order of the finite element discretization space

Select the order of the finite element discretization space. For NURBS meshes, we increase the order by degree elevation.

   if (mesh->NURBSext)
   {
      mesh->DegreeElevate(order, order);
   }

This code checks if the mesh has NURBS extension by checking if mesh->NURBSext is true or not. The NURBS extension is an optional part of MFEM that allows for the representation and manipulation of Non-Uniform Rational B-Splines (NURBS) in the library.

If the mesh has NURBS extension, it calls the DegreeElevate() method on the mesh object, passing in the values of the “order” variable as the first and second parameters. The DegreeElevate method is used to elevate the degree of the NURBS elements of the mesh. The first parameter represents the degree elevation in the u-direction, and the second parameter represents the degree elevation in the v-direction. By passing in the “order” variable as the parameter, it means the degree of NURBS element is elevated to the order specified by the user.

This is likely to increase the accuracy of the solution by increasing the polynomial degree of the finite elements used to discretize the PDE. However, this step is only performed if the mesh has NURBS extension and if not, it will be skipped.

Refine the mesh to increase the resolution

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 5,000 elements.

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

This code is performing uniform refinement on the mesh. The refinement is done by repeatedly calling the UniformRefinement() method on the mesh object. The number of times the method is called, or the number of refinement levels, is determined by the variable “ref_levels”.

The value of ref_levels is computed by the formula: ref_levels = (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);

Here:

• log(5000./mesh->GetNE()) is the base-2 logarithm of the ratio of a target number of elements (5000) to the current number of elements in the mesh.
• log(2.) is the base-2 logarithm of 2
• dim is the dimension of the mesh.
• floor() is a C++ function that rounds down to the nearest integer.

This formula aims to make the number of elements in the mesh close to 5000, while taking into account the dimension of the mesh. The larger the dimension, the more elements are needed to represent the same physical area, therefore the more refinement is needed.

The for loop then iterates over “ref_levels” and calls the UniformRefinement() method on the mesh object each time. This method sub-divides each element of the mesh by creating new nodes at the midpoints of the edges and subdividing the element into smaller ones. This process increases the number of elements in the mesh and can help to improve the accuracy of the solution by providing a finer discretization of the domain.