Hyperelasticity (C++)

Background

See the section Hyperelasticity background for some mathematical background on this demo.

Implementation

The implementation is split in two files: a form file containing the definition of the variational forms expressed in UFL and a C++ file containing the actual solver.

Running this demo requires the files: main.cpp, HyperElasticity.ufl and CMakeLists.txt.

UFL form file

The UFL file is implemented in Hyperelasticity.ufl, and the explanation of the UFL file can be found at here.

C++ program

The main solver is implemented in the main.cpp file.

At the top, we include the DOLFIN header file and the generated header file “HyperElasticity.h” containing the variational forms and function spaces. For convenience we also include the DOLFIN namespace.

#include <dolfin.h>
#include "HyperElasticity.h"

using namespace dolfin;

We begin by defining two classes, deriving from SubDomain for later use when specifying domains for the boundary conditions.

// Sub domain for clamp at left end
class Left : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return (std::abs(x[0]) < DOLFIN_EPS) && on_boundary;
  }
};

// Sub domain for rotation at right end
class Right : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return (std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary;
  }
};

We also define two classes, deriving from Expression, for later use when specifying values for the boundary conditions.

// Dirichlet boundary condition for clamp at left end
class Clamp : public Expression
{
public:

  Clamp() : Expression(3) {}

  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 0.0;
    values[1] = 0.0;
    values[2] = 0.0;
  }

};

// Dirichlet boundary condition for rotation at right end
class Rotation : public Expression
{
public:

  Rotation() : Expression(3) {}

  void eval(Array<double>& values, const Array<double>& x) const
  {
    const double scale = 0.5;

    // Center of rotation
    const double y0 = 0.5;
    const double z0 = 0.5;

    // Large angle of rotation (60 degrees)
    double theta = 1.04719755;

    // New coordinates
    double y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta);
    double z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta);

    // Rotate at right end
    values[0] = 0.0;
    values[1] = scale*(y - x[1]);
    values[2] = scale*(z - x[2]);
  }
};

Next:

int main()
{

Inside the main function, we begin by defining a tetrahedral mesh of the domain and the function space on this mesh. Here, we choose to create a unit cube mesh with 25 ( = 24 + 1) verices in one direction and 17 ( = 16 + 1) vertices in the other two directions. With this mesh, we initialize the (finite element) function space defined by the generated code.

// Create mesh and define function space
auto mesh = std::make_shared<UnitCubeMesh>(24, 16, 16);
auto V = std::make_shared<HyperElasticity::FunctionSpace>(mesh);

Now, the Dirichlet boundary conditions can be created using the class DirichletBC, the previously initialized FunctionSpace V and instances of the previously listed classes Left (for the left boundary) and Right (for the right boundary), and Clamp (for the value on the left boundary) and Rotation (for the value on the right boundary).

// Define Dirichlet boundaries
auto left = std::make_shared<Left>();
auto right = std::make_shared<Right>();

// Define Dirichlet boundary functions
auto c = std::make_shared<Clamp>();
auto r = std::make_shared<Rotation>();

// Create Dirichlet boundary conditions
DirichletBC bcl(V, c, left);
DirichletBC bcr(V, r, right);
std::vector<const DirichletBC*> bcs = {{&bcl, &bcr}};

The two boundary conditions are collected in the container bcs.

We use two instances of the class Constant to define the source B and the traction T.

// Define source and boundary traction functions
auto B = std::make_shared<Constant>(0.0, -0.5, 0.0);
auto T = std::make_shared<Constant>(0.1,  0.0, 0.0);

The solution for the displacement will be an instance of the class Function, living in the function space V; we define it here:

// Define solution function
auto u = std::make_shared<Function>(V);

Next, we set the material parameters

// Set material parameters
const double E  = 10.0;
const double nu = 0.3;
auto mu = std::make_shared<Constant>(E/(2*(1 + nu)));
auto lambda = std::make_shared<Constant>(E*nu/((1 + nu)*(1 - 2*nu)));

Now, we can initialize the bilinear and linear forms (a, L) using the previously defined FunctionSpace V. We attach the material parameters and previously initialized functions to the forms.

// Create (linear) form defining (nonlinear) variational problem
HyperElasticity::ResidualForm F(V);
F.mu = mu; F.lmbda = lambda; F.u = u;
F.B = B; F.T = T;

// Create Jacobian dF = F' (for use in nonlinear solver).
HyperElasticity::JacobianForm J(V, V);
J.mu = mu; J.lmbda = lambda; J.u = u;

Now, we have specified the variational forms and can consider the solution of the variational problem.

// Solve nonlinear variational problem F(u; v) = 0
solve(F == 0, *u, bcs, J);

Finally, the solution u is saved to a file named displacement.pvd in VTK format.

  // Save solution in VTK format
  File file("displacement.pvd");
  file << *u;

  return 0;
}