Linear elasticity

Table of Contents

• Introduction
• Setup and initialization
• Problem data and material parameters
• Grid creation
• Setting up the Dirichlet boundary
• Preparing thread-parallelism
• Selection of the finite element basis
• Matrix, vector, and constraints data structures
• Definition of the variational problem
• Global assembly of the problem
• Algebraic solution
• Writing the solution to VTK
• Error handling

Full code of this example This example on gitlab

Introduction

The following example shows how to solve a linear elasticity problem with Dirichlet boundary conditions on a subset of the boundary.

Setup and initialization

After including the required headers the main() function first creates a simple logger that allows to write log messages to std::cout. The logger will prefix log messages with the time passed since the program was started and since the last log message using the given format string. Then the command line arguments are parsed. They allow to set the number of grid refinement steps, the number of in parallel algorithms, and parameters of the iterative linear solver. Afterwards MPI is initialized, passing the command line arguments.

int main (int argc, char *argv[]) try
{
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto refinements  = parameters.get<std::size_t>("refinements", 1);
  auto threadCount  = parameters.get<std::size_t>("threads", std::thread::hardware_concurrency());
#if not(HAVE_SUITESPARSE_CHOLMOD)
  auto tol          = parameters.get<double>("tol", 1e-4);
  auto maxIter      = parameters.get<std::size_t>("maxIter", 1000);
  auto verbosity    = parameters.get<std::size_t>("verbosity", 2);
#endif
 
  log("Command line processed");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");

Problem data and material parameters

The example solves a linear elasticity problem in 3 dimensions. For convenience we first define the spatial dimension, and introduce aliases for corresponding local vector and square matrix types. for a steel beam with length 1m and a quadratic profile of thickness 1cm. First we define the spatial dimension and extend of the domain.

  const int dim = 3;
 
  using LocalVector = Dune::FieldVector<double,dim>;
  using LocalMatrix = Dune::FieldMatrix<double,dim,dim>;

We use a linear St.Venant-Kirchhoff model which is parameterized by the Lamé parameters \(\lambda\) und \(\mu\). Here we compute them from Youngs modulus \(E_{Young}\) and the Poisson ration \(\nu\). Furthermore we define the density \(\rho\) and gravitational acceleration vector \(g\).

  double E_Young = 200e9; // in Pa = N/m^2 = 10^9 GPa
  double nu = 0.3;
 
  double lambda = E_Young*nu/((1.-2.*nu) * (1.+nu));
  double mu = E_Young/(2*(1+nu));
 
  double rho = 7.85e3;    // in kg/m^3 = 10^3 g/cm^3
  double g_n = 9.80665;   // in m/s^2 = N/kg
  auto g = LocalVector({0.0, -g_n*rho, 0.0});

Now we specify the extend of the computational domain, an indicator function for the Dirichlet boundary, and the Dirichlet boundary values.

  double thickness = .01;
  double length = 1.0;
 
  auto dirichletIndicator = [](auto x) { return x[0] <= 1e-10; };
  auto dirichletValues = [](const auto& x){ return LocalVector(0.0); };

Grid creation

Since we will use a uniform grid we select the YaspGrid implementation which supports structured uniform grids only. In contrast to an implementation supporting unstructured grids this is more efficient since all grid quantities can be computed on the fly. To create the grid we use a StructuredGridFactory which only requires to specify two opposite corners of the domain and the number of elements in each direction. The createCubeGrid() function returns a std::unique_ptr referring to the newly creates grid. For simplicity we also store a reference to this grid. The coarse grid is then refined globally and a leaf grid view is obtained.

  using Grid = Dune::YaspGrid<dim>;
 
  auto lowerCorner = LocalVector({0, 0, 0});
  auto upperCorner = LocalVector({length, thickness, thickness});
  auto elementCount = std::array{(unsigned int)(length/thickness), 1u, 1u};
 
  auto gridPtr = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerCorner, upperCorner, elementCount);
  auto& grid = *gridPtr;
 
  log("Grid created");
 
  grid.globalRefine(refinements);
 
  log(Dune::formatString("Grid refined %d times", refinements));
 
  auto gridView = grid.leafGridView();

Setting up the Dirichlet boundary

To specify a Dirichlet boundary condition on a subset of the boundary, we create a Dune::Fufem::BoundaryPatch for this boundary subset and insert all boundary intersections that are identified by the indicator function. Here we use the Dune::Fufem::Boundary class which provides a convenient way to traverse all boundary intersections.

  auto dirichletPatch = Dune::Fufem::BoundaryPatch(gridView);
  for (const auto& intersection : Dune::Fufem::Boundary(gridView))
    if (dirichletIndicator(intersection.geometry().center()))
      dirichletPatch.insertFace(intersection);
 
  log("Boundary patch set up");

Preparing thread-parallelism

Dune-fufem supports thread-parallel assembly based on a coloring of the grid view. To use this feature we need to compute such a coloring. The coloredElementRange() function from the dune-assembler module computes such a coloring and returns a partition of the grid view's elements according to the coloring. The ColoredRangeExecutor class then encapsulates the information for executing parallel algorithms based on this partition.

  auto coloredRange = Dune::Assembler::Experimental::coloredElementRange(gridView);
 
  log("Grid coloring computed");
 
  auto executor = Dune::Assembler::Experimental::ColoredRangeExecutor(coloredRange, threadCount);
 
  log(Dune::formatString("Using parallelism with %d threads", threadCount));

Selection of the finite element basis

As a next step we create a global finite element basis on the GridView. The helper function lagrange<2>() declares that we want to use Lagrange finite elements of order two. Since we want to compute a 3d vector field, we take a three-fold cartesian product of this space and use makeBasis() to create the corresponding basis on the given GridView.

The resulting basis will contain three copies of each scalar Lagrange basis function, one for each component. The additional argument of the power<dim>(...) function allows to control how these basis functions are indexed. Here we use the blockedInterleaved strategy which will append the index digits 0, 1, 2 to the three copies of the scalar basis function with index i leading to multi-indices (i,0), (i,1), (i,2).

  const int order = 2;
 
  using namespace Dune::Functions::BasisFactory;
  auto basis = makeBasis(gridView, power<dim>(lagrange<order>(), blockedInterleaved()));
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));

Matrix, vector, and constraints data structures

Since the basis' indices are multi-indices with two entries, we introduce aliases for blocked vector and matrix types with small 3d blocks and create a solution and a right-hand-side vector and a stiffness matrix.

  using Vector = Dune::BlockVector<LocalVector>;
  using BitVector = std::vector<std::array<char,dim>>;
  using Matrix = Dune::BCRSMatrix<LocalMatrix>;
 
  auto sol = Vector();
  auto rhs = Vector();
  auto stiffnessMatrix = Matrix();

To handle the Dirichlet boundary conditions we then create an object for storing affine constraints with respect to the basis and then fill the object by computing the actual constraints values on the boundary patch from the Dirichlet value function.

  auto constraints = Dune::Fufem::makeAffineConstraints<BitVector, Vector>(basis);
 
  Dune::Fufem::computeBoundaryConstraints(constraints, basis, dirichletValues, dirichletPatch);
 
  log("Boundary constraints computed");

Definition of the variational problem

We define the local assembler using the expression language from Dune::Fufem::Forms namespace. After including the namespace, we first obtain expressions for trial and test functions associated to the basis. To define the integrant of the bilinear form for a St. Venant-Kirchhoff material, we first introduce the short-cut E for the symmetric gradient operator. Using a local identity matrix we can then define the stress tensor \(\sigma\) in dependence of the trial function. Here we make use of the fact that the dot(...) function of Dune::Fufem::Forms implements the Frobenius inner product when applied to matrices.

  using namespace Dune::Fufem::Forms;
 
  auto u = trialFunction(basis);
  auto v = testFunction(basis);
 
  auto Id = Dune::ScaledIdentityMatrix<double,dim>{1.0};
 
  auto E = [&](const auto& v) {
    return symmetrize(grad(v));
  };
 
  auto sigma = 2*mu*E(u) + lambda*Id*trace(E(u));
 
  auto a = integrate(dot(sigma, E(v)));
  auto l = integrate(dot(g,v));

Note

Alternatively we could use (and this is even a bit faster)

auto C = RangeOperator([&](const auto& e) {
return 2*mu*e + lambda*LocalMatrix(Id)*trace(e);
});
auto a = integrate(dot(C(E(u)), E(v)));

Global assembly of the problem

For global assembly we use the dune-assembler module. This module allows to use different linear algebra frameworks through a unified matrix- and vector-backend interface. To use the matrix- and vector-classes from dune-istl, we wrap them into corresponding ISTL-backends.

  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);

Next we create a global assembler object using the same basis for the test and trial space. In order to use thread-parallel assembly, we additionally pass the parallel executor which encapsulates the grid coloring and requested thread count to the assembler.

  auto assembler = Dune::Assembler::Assembler(basis, basis, executor);
 
  log("Assembler set up");

As a first step we need to assemble the matrix pattern. To this end we obtain a patternBuilder object from the matrix backend and initialize its size according to the test and trial basis. Then we use the global assemblers assembleMatrixPattern() method to assemble the pattern. While local matrices are often fully populated, this is not necessarily the case. Hence the local assembler also defines a local pattern and thus has to be passed when assembling the pattern, too.

  auto patternBuilder = matrixBackend.patternBuilder();
 
  patternBuilder.resize(basis, basis);
 
  assembler.assembleMatrixPattern(a, patternBuilder);
 
  log("Matrix pattern assembled");

Now we are ready to assemble the problem. To this end we let the patternBuilder initialize the matrix size and pattern and then call the actual assembly method for the matrix.

  patternBuilder.setupMatrix();
 
  log("Matrix set up");
 
  assembler.assembleMatrixEntries(a, matrixBackend);
 
  log("Matrix assembled");

Similarly we assemble the right hand side by first resizing the vector according to the test basis and initializing it with zero and then calling the corresponding assembly method.

  rhsBackend.resize(basis);
  rhsBackend.setZero();
 
  assembler.assembleVectorEntries(l, rhsBackend);
 
  log("RHS assembled");

Once the system is assembled, we can constrain it to the affine subspace satisfying the constraints that we defined earlier.

  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
 
  log("Boundary condition incorporated into system");

Algebraic solution

Before solving the linear system, we have to initialize the solution vector. While we could directly use the corresponding resize method and assignment operator, we here demonstrate how the ISTLVectorBackend from dune-assembler can be used, since this approach would also work for basis using multi-indices and corresponding nested vector containers.

  Dune::Assembler::ISTLVectorBackend(sol).resize(basis);
  Dune::Assembler::ISTLVectorBackend(sol).setZero();

We will solve the linear system using the sparse direct solver Choldmod, if available. Otherwise we use a CG solver with incomplete LDLT decomposition as preconditioner.

#if HAVE_SUITESPARSE_CHOLMOD
  auto solver = Dune::Cholmod<Vector>();
  solver.setMatrix(stiffnessMatrix);
#else
  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
  auto preconditioner = Dune::SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix, 1.0);
 
  log("Preconditioner created");
 
  auto solver = Dune::CGSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  tol,
                  maxIter,
                  verbosity
                );
#endif
 
  log("Solver created");
 
  auto solverStatistics = Dune::InverseOperatorResult();
 
  solver.apply(sol, rhs, solverStatistics);
 
  log("Linear system solved");

Writing the solution to VTK

Finally we can write the solution to a VTK file understood by the ParaView visualization program using the infrastructure from the dune-vtk module. To this end we need to create finite element functions by combining the basis with the coefficient vector. Here we can again use Dune::Fufem::Forms expressions and bind them to a coefficient vector using the expr|sol syntax resulting in grid functions that can later be passed to the writer object. Additionally to the deformation u and the tensors E(u) sigma that have already been defined, we will also write the volumetric, deviatoric, and von Mises stress.

  auto sigma_vol = 1./3. * trace(sigma);
  auto sigma_dev = sigma - sigma_vol*Id;
  auto sqrt = RangeOperator([](const auto& z) { return std::sqrt(z); });
  auto sigma_von_Mises_sol = sqrt(3./2.*Dune::Fufem::Forms::dot(sigma_dev|sol, sigma_dev|sol));

Note

In the notion of Dune::Fufem::Forms expressions like u, E(u), sigma represent linear operators. Binding them to a coefficient vector corresponds to applying the operator to the function associated with the coefficient vector and thus results in a grid function (also called bound operator in the following).

Note

While one can either combine several linear operators and then bind the result or combine several already bound operators, non-linear operations can only be applied after binding. Hence we have to bind sigma_dev, before the application of the nonlinear dot(...) operation. Furthermore we have to make std::sqrt Dune::Fufem::Forms-aware by wrapping it into a Dune::Fufem::Forms::RangeOperator. Applying this to a bound operator produces another bound operator where the wrapped function is applied pointwise.

Now we create a writer capable of writing unstructured grids and discontinuous polynomial data of the desired order, pass it all the functions of interest, and finally write the result to a file.

  auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::DiscontinuousLagrangeDataCollector(gridView, order));
  vtkWriter.addPointData(u|sol, Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::vector, dim));
  vtkWriter.addPointData(E(u)|sol, Dune::VTK::FieldInfo("Eu", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma|sol, Dune::VTK::FieldInfo("sigma", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma_dev|sol, Dune::VTK::FieldInfo("sigma_dev", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma_vol|sol, Dune::VTK::FieldInfo("sigma_vol", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(sigma_von_Mises_sol, Dune::VTK::FieldInfo("sigma_von_Mises", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write("linear-elasticity");
 
  log("Solution written to vtk file");

Error handling

In case an exception is thrown within main() it is followed by a catch block that will print the error message and exit.

}
catch (Dune::Exception& e) {
  std::cout << e.what() << std::endl;
  return 1;
}