text:
Reference documentation for deal.II version Git e7c9a24 2019-02-08 08:26:19 +0100
The 'Quasi-Static Finite-Strain Compressible Elasticity' code gallery program

This program was contributed by Jean-Paul Pelteret <jppelteret@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

# Pictures from this code gallery program

## Overview

The Cook membrane (or cantilever) problem is a classic benchmark test for finite element formulations for solid mechanics. It is typically used to test for and demonstrate the shear-locking (or locking-free) behaviour of a finite element ansatz under quasi-incompressible conditions. As it is so widely referred to in the literature on finite-strain elasticity, we reproduce the example here. However, we consider on the compressible case to avoid many of the complexities that arise in @ref step_44 "step-44", which provides an efficient approach to deal with the quasi-incompressible case.

### A classical approach to solving the cook membrane problem.

In this work we take a classical approach to solving the equations governing quasi-static finite-strain compressible elasticity, with code based on @ref step_44 "step-44". The formulation adopted here is that seen in many texts on solid mechanics and can be used as the starting point for extension into many topics such as material anisotropy, rate dependence or plasticity, or even as a component of multi-physics problems.

The basic problem configuration is summarised in the following image. A beam of specific dimensions is fixed at one end and a uniform traction load is applied at the other end such that the total force acting on this surface totals 1 Newton. Displacement in the third coordinate direction (out of plane) is prevented in order to impose plane strain conditions.

Note that we perform a three-dimensional computation as, for this particular formulation, the two-dimensional case corresponds to neither plane-strain nor plane-stress conditions.

## Requirements

Version 8.5 or greater of deal.II

## Compiling and running

Similar to the example programs, run

1 cmake -DDEAL_II_DIR=/path/to/deal.II .

in this directory to configure the problem. You can switch between debug and release mode by calling either

1 make debug

or

1 make release

The problem may then be run with

1 make run

## Reference for this work

If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: J-P. V. Pelteret and A. McBride, The deal.II code gallery: Quasi-Static Finite-Strain Compressible Elasticity, 2016. DOI: 10.5281/zenodo.1228964

### Acknowledgements

The support of this work by the European Research Council (ERC) through the Advanced Grant 289049 MOCOPOLY is gratefully acknowledged by the first author.

## Recommended Literature

C. Miehe (1994), Aspects of the formulation and finite element implementation of large strain isotropic elasticity International Journal for Numerical Methods in Engineering 37 , 12, 1981-2004. DOI: 10.1002/nme.1620371202; G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: 978-0-471-82319-3; P. Wriggers (2008), Nonlinear finite element methods, Springer. DOI: 10.1007/978-3-540-71001-1; T.J.R. Hughes (2000), The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover. ISBN: 978-0486411811

The derivation of the finite-element problem, namely the definition and linearisation of the residual and their subsequent discretisation are quite lengthy and involved. Thankfully, the classical approach adopted in this work is well documented and therefore does not need to be reproduced here. We refer the reader to, among many other possible texts, Holzapfel (2001) and Wriggers (2008) for a detailed description of the approach applied in this work. It amounts to a reduction and slight reworking of @ref step_44 "step-44" (accounting for the removal of the two additional fields used therein). We also refer the reader to @ref step_44 "step-44" for a brief overview of the continuum mechanics and kinematics related to solid mechanics. We also provide two alternative assembly mechanisms, which can be selected within the parameter file, that use automatic differentiation to assemble the linear system. These not only demonstrate how one might employ such automated techniques to linearise the residual, but also serve as a verification of the correctness of the hand-developed implementation of the tangent stiffness matrix.

## Results

These results were produced using the following material properties: Shear modulus is 422.5kPa Poisson ratio is 0.3

The 32x32x1 discretised reference geometry looks as follows:

And an example of the displaced solution is given in the next image.

Below we briefly document the tip displacement as predicted for different discretisation levels and ansatz for the displacement field. A direct and, by visual inspection, favourable comparison of the following results can be made with those found in Miehe (1994). Since the material is compressible, shear-locking is not exhibited by the beam for low-order elements.

#### Number of degrees of freedom

Elements per edge Q1 Q2
1 24 81
2 54 225
4 150 729
8 486 2601
16 1734 9801
32 6534 38025
64 25350 149769

#### Tip y-displacement (in mm)

Elements per edge Q1 Q2
1 5.15 12.19
2 8.72 13.83
4 12.02 14.22
8 13.61 14.30
16 14.13 14.32
32 14.28 14.33
64 14.32 14.33

# Annotated version of cook_membrane.cc

/ * ---------------------------------------------------------------------
* Copyright (C) 2010 - 2015 by the deal.II authors and
* Jean-Paul Pelteret and Andrew McBride
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
* /
/ *
* Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg,
* Andrew McBride, University of Cape Town, 2015, 2017
* /

We start by including all the necessary deal.II header files and some C++ related ones. They have been discussed in detail in previous tutorial programs, so you need only refer to past tutorials for details.

#include <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition_selector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/config.h>
#if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS)
#endif

These must be included below the AD headers so that their math functions are available for use in the definition of tensors and kinematic quantities

#include <deal.II/physics/elasticity/kinematics.h>
#include <deal.II/physics/elasticity/standard_tensors.h>
#include <iostream>
#include <fstream>

We then stick everything that relates to this tutorial program into a namespace of its own, and import all the deal.II function and class names into it:

namespace Cook_Membrane
{
using namespace dealii;

### Run-time parameters

There are several parameters that can be set in the code so we set up a ParameterHandler object to read in the choices at run-time.

namespace Parameters
{

#### Assembly method

Here we specify whether automatic differentiation is to be used to assemble the linear system, and if so then what order of differentiation is to be employed.

struct AssemblyMethod
{
unsigned int automatic_differentiation_order;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void AssemblyMethod::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Assembly method");
{
prm.declare_entry("Automatic differentiation order", "0",
"The automatic differentiation order to be used in the assembly of the linear system.\n"
"# Order = 0: Both the residual and linearisation are computed manually.\n"
"# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
"# Order = 2: Both the residual and linearisation are computed using AD.");
}
}
void AssemblyMethod::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Assembly method");
{
automatic_differentiation_order = prm.get_integer("Automatic differentiation order");
}
}

#### Finite Element system

Here we specify the polynomial order used to approximate the solution. The quadrature order should be adjusted accordingly.

struct FESystem
{
unsigned int poly_degree;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void FESystem::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Finite element system");
{
prm.declare_entry("Polynomial degree", "2",
"Displacement system polynomial order");
}
}
void FESystem::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Finite element system");
{
poly_degree = prm.get_integer("Polynomial degree");
}
}

#### Geometry

Make adjustments to the problem geometry and its discretisation.

struct Geometry
{
unsigned int elements_per_edge;
double scale;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Geometry::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Geometry");
{
prm.declare_entry("Elements per edge", "32",
"Number of elements per long edge of the beam");
prm.declare_entry("Grid scale", "1e-3",
"Global grid scaling factor");
}
}
void Geometry::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Geometry");
{
elements_per_edge = prm.get_integer("Elements per edge");
scale = prm.get_double("Grid scale");
}
}

#### Materials

We also need the shear modulus $$\mu$$ and Poisson ration $$\nu$$ for the neo-Hookean material.

struct Materials
{
double nu;
double mu;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Materials::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
prm.declare_entry("Poisson's ratio", "0.3",
Patterns::Double(-1.0,0.5),
"Poisson's ratio");
prm.declare_entry("Shear modulus", "0.4225e6",
"Shear modulus");
}
}
void Materials::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
nu = prm.get_double("Poisson's ratio");
mu = prm.get_double("Shear modulus");
}
}

#### Linear solver

Next, we choose both solver and preconditioner settings. The use of an effective preconditioner is critical to ensure convergence when a large nonlinear motion occurs within a Newton increment.

struct LinearSolver
{
std::string type_lin;
double tol_lin;
double max_iterations_lin;
std::string preconditioner_type;
double preconditioner_relaxation;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void LinearSolver::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Linear solver");
{
prm.declare_entry("Solver type", "CG",
Patterns::Selection("CG|Direct"),
"Type of solver used to solve the linear system");
prm.declare_entry("Residual", "1e-6",
"Linear solver residual (scaled by residual norm)");
prm.declare_entry("Max iteration multiplier", "1",
"Linear solver iterations (multiples of the system matrix size)");
prm.declare_entry("Preconditioner type", "ssor",
Patterns::Selection("jacobi|ssor"),
"Type of preconditioner");
prm.declare_entry("Preconditioner relaxation", "0.65",
"Preconditioner relaxation value");
}
}
void LinearSolver::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Linear solver");
{
type_lin = prm.get("Solver type");
tol_lin = prm.get_double("Residual");
max_iterations_lin = prm.get_double("Max iteration multiplier");
preconditioner_type = prm.get("Preconditioner type");
preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
}
}

#### Nonlinear solver

A Newton-Raphson scheme is used to solve the nonlinear system of governing equations. We now define the tolerances and the maximum number of iterations for the Newton-Raphson nonlinear solver.

struct NonlinearSolver
{
unsigned int max_iterations_NR;
double tol_f;
double tol_u;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void NonlinearSolver::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Nonlinear solver");
{
prm.declare_entry("Max iterations Newton-Raphson", "10",
"Number of Newton-Raphson iterations allowed");
prm.declare_entry("Tolerance force", "1.0e-9",
"Force residual tolerance");
prm.declare_entry("Tolerance displacement", "1.0e-6",
"Displacement error tolerance");
}
}
void NonlinearSolver::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Nonlinear solver");
{
max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
tol_f = prm.get_double("Tolerance force");
tol_u = prm.get_double("Tolerance displacement");
}
}

#### Time

Set the timestep size $$\varDelta t$$ and the simulation end-time.

struct Time
{
double delta_t;
double end_time;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Time::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Time");
{
prm.declare_entry("End time", "1",
"End time");
prm.declare_entry("Time step size", "0.1",
"Time step size");
}
}
void Time::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Time");
{
end_time = prm.get_double("End time");
delta_t = prm.get_double("Time step size");
}
}

#### All parameters

Finally we consolidate all of the above structures into a single container that holds all of our run-time selections.

struct AllParameters :
public AssemblyMethod,
public FESystem,
public Geometry,
public Materials,
public LinearSolver,
public NonlinearSolver,
public Time
{
AllParameters(const std::string &input_file);
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
AllParameters::AllParameters(const std::string &input_file)
{
declare_parameters(prm);
prm.parse_input(input_file);
parse_parameters(prm);
}
void AllParameters::declare_parameters(ParameterHandler &prm)
{
AssemblyMethod::declare_parameters(prm);
FESystem::declare_parameters(prm);
Geometry::declare_parameters(prm);
Materials::declare_parameters(prm);
LinearSolver::declare_parameters(prm);
NonlinearSolver::declare_parameters(prm);
Time::declare_parameters(prm);
}
void AllParameters::parse_parameters(ParameterHandler &prm)
{
AssemblyMethod::parse_parameters(prm);
FESystem::parse_parameters(prm);
Geometry::parse_parameters(prm);
Materials::parse_parameters(prm);
LinearSolver::parse_parameters(prm);
NonlinearSolver::parse_parameters(prm);
Time::parse_parameters(prm);
}
}

### Time class

A simple class to store time data. Its functioning is transparent so no discussion is necessary. For simplicity we assume a constant time step size.

class Time
{
public:
Time (const double time_end,
const double delta_t)
:
timestep(0),
time_current(0.0),
time_end(time_end),
delta_t(delta_t)
{}
virtual ~Time()
{}
double current() const
{
return time_current;
}
double end() const
{
return time_end;
}
double get_delta_t() const
{
return delta_t;
}
unsigned int get_timestep() const
{
return timestep;
}
void increment()
{
time_current += delta_t;
++timestep;
}
private:
unsigned int timestep;
double time_current;
const double time_end;
const double delta_t;
};

### Compressible neo-Hookean material within a one-field formulation

As discussed in the literature and step-44, Neo-Hookean materials are a type of hyperelastic materials. The entire domain is assumed to be composed of a compressible neo-Hookean material. This class defines the behaviour of this material within a one-field formulation. Compressible neo-Hookean materials can be described by a strain-energy function (SEF) $$\Psi = \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J)$$.

The isochoric response is given by $$\Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3]$$ where $$c_{1} = \frac{\mu}{2}$$ and $$\overline{I}_{1}$$ is the first invariant of the left- or right-isochoric Cauchy-Green deformation tensors. That is $$\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})$$. In this example the SEF that governs the volumetric response is defined as $$\Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1 - 2\textrm{ln}\; J ]$$, where $$\kappa:= \lambda + 2/3 \mu$$ is the bulk modulus and $$\lambda$$ is Lame's first parameter.

The following class will be used to characterize the material we work with, and provides a central point that one would need to modify if one were to implement a different material model. For it to work, we will store one object of this type per quadrature point, and in each of these objects store the current state (characterized by the values or measures of the displacement field) so that we can compute the elastic coefficients linearized around the current state.

template <int dim,typename NumberType>
class Material_Compressible_Neo_Hook_One_Field
{
public:
Material_Compressible_Neo_Hook_One_Field(const double mu,
const double nu)
:
kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
c_1(mu / 2.0)
{
Assert(kappa > 0, ExcInternalError());
}
~Material_Compressible_Neo_Hook_One_Field()
{}

The first function is the total energy $$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}$$.

NumberType
get_Psi(const NumberType &det_F,
{
return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
}

The second function determines the Kirchhoff stress $$\boldsymbol{\tau} = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}$$

get_tau(const NumberType &det_F,
{

See Holzapfel p231 eq6.98 onwards

return get_tau_vol(det_F) + get_tau_iso(b_bar);
}

The fourth-order elasticity tensor in the spatial setting $$\mathfrak{c}$$ is calculated from the SEF $$\Psi$$ as $$J \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}$$ where $$\mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial \mathbf{C} \partial \mathbf{C}}$$

get_Jc(const NumberType &det_F,
{
return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
}
private:

Define constitutive model parameters $$\kappa$$ (bulk modulus) and the neo-Hookean model parameter $$c_1$$:

const double kappa;
const double c_1;

Value of the volumetric free energy

NumberType
get_Psi_vol(const NumberType &det_F) const
{
return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
}

Value of the isochoric free energy

NumberType
get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
{
return c_1 * (trace(b_bar) - dim);
}

Derivative of the volumetric free energy with respect to $$J$$ return $$\frac{\partial \Psi_{\text{vol}}(J)}{\partial J}$$

NumberType
get_dPsi_vol_dJ(const NumberType &det_F) const
{
return (kappa / 2.0) * (det_F - 1.0 / det_F);
}

The following functions are used internally in determining the result of some of the public functions above. The first one determines the volumetric Kirchhoff stress $$\boldsymbol{\tau}_{\textrm{vol}}$$. Note the difference in its definition when compared to step-44.

get_tau_vol(const NumberType &det_F) const
{
return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
}

Next, determine the isochoric Kirchhoff stress $$\boldsymbol{\tau}_{\textrm{iso}} = \mathcal{P}:\overline{\boldsymbol{\tau}}$$:

Then, determine the fictitious Kirchhoff stress $$\overline{\boldsymbol{\tau}}$$:

get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
{
return 2.0 * c_1 * b_bar;
}

Second derivative of the volumetric free energy wrt $$J$$. We need the following computation explicitly in the tangent so we make it public. We calculate $$\frac{\partial^2 \Psi_{\textrm{vol}}(J)}{\partial J \partial J}$$

NumberType
get_d2Psi_vol_dJ2(const NumberType &det_F) const
{
return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
}

Calculate the volumetric part of the tangent $$J \mathfrak{c}_\textrm{vol}$$. Again, note the difference in its definition when compared to step-44. The extra terms result from two quantities in $$\boldsymbol{\tau}_{\textrm{vol}}$$ being dependent on $$\boldsymbol{F}$$.

get_Jc_vol(const NumberType &det_F) const
{

See Holzapfel p265

return det_F
* ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
- (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
}

Calculate the isochoric part of the tangent $$J \mathfrak{c}_\textrm{iso}$$:

get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
{
const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
const SymmetricTensor<4, dim> tau_iso_x_I
= outer_product(tau_iso,
const SymmetricTensor<4, dim> I_x_tau_iso
tau_iso);
const SymmetricTensor<4, dim> c_bar = get_c_bar();
return (2.0 / dim) * trace(tau_bar)
- (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
}

Calculate the fictitious elasticity tensor $$\overline{\mathfrak{c}}$$. For the material model chosen this is simply zero:

get_c_bar() const
{
}
};

As seen in step-18, the  PointHistory  class offers a method for storing data at the quadrature points. Here each quadrature point holds a pointer to a material description. Thus, different material models can be used in different regions of the domain. Among other data, we choose to store the Kirchhoff stress $$\boldsymbol{\tau}$$ and the tangent $$J\mathfrak{c}$$ for the quadrature points.

template <int dim,typename NumberType>
class PointHistory
{
public:
PointHistory()
{}
virtual ~PointHistory()
{}

The first function is used to create a material object and to initialize all tensors correctly: The second one updates the stored values and stresses based on the current deformation measure $$\textrm{Grad}\mathbf{u}_{\textrm{n}}$$.

void setup_lqp (const Parameters::AllParameters &parameters)
{
material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
parameters.nu));
}

We offer an interface to retrieve certain data. This is the strain energy:

NumberType
get_Psi(const NumberType &det_F,
{
return material->get_Psi(det_F,b_bar);
}

Here are the kinetic variables. These are used in the material and global tangent matrix and residual assembly operations: First is the Kirchhoff stress:

get_tau(const NumberType &det_F,
{
return material->get_tau(det_F,b_bar);
}

And the tangent:

get_Jc(const NumberType &det_F,
{
return material->get_Jc(det_F,b_bar);
}

In terms of member functions, this class stores for the quadrature point it represents a copy of a material type in case different materials are used in different regions of the domain, as well as the inverse of the deformation gradient...

private:
std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
};

### Quasi-static compressible finite-strain solid

Forward declarations for classes that will perform assembly of the linear system.

template <int dim,typename NumberType>
struct Assembler_Base;
template <int dim,typename NumberType>
struct Assembler;

The Solid class is the central class in that it represents the problem at hand. It follows the usual scheme in that all it really has is a constructor, destructor and a run() function that dispatches all the work to private functions of this class:

template <int dim,typename NumberType>
class Solid
{
public:
Solid(const Parameters::AllParameters &parameters);
virtual
~Solid();
void
run();
private:

We start the collection of member functions with one that builds the grid:

void
make_grid();

Set up the finite element system to be solved:

void
system_setup();

Several functions to assemble the system and right hand side matrices using multithreading. Each of them comes as a wrapper function, one that is executed to do the work in the WorkStream model on one cell, and one that copies the work done on this one cell into the global object that represents it:

void
assemble_system(const BlockVector<double> &solution_delta);

We use a separate data structure to perform the assembly. It needs access to some low-level data, so we simply befriend the class instead of creating a complex interface to provide access as necessary.

friend struct Assembler_Base<dim,NumberType>;
friend struct Assembler<dim,NumberType>;

Apply Dirichlet boundary conditions on the displacement field

void
make_constraints(const int &it_nr);

Create and update the quadrature points. Here, no data needs to be copied into a global object, so the copy_local_to_global function is empty:

void
setup_qph();

Solve for the displacement using a Newton-Raphson method. We break this function into the nonlinear loop and the function that solves the linearized Newton-Raphson step:

void
solve_nonlinear_timestep(BlockVector<double> &solution_delta);
std::pair<unsigned int, double>
solve_linear_system(BlockVector<double> &newton_update);

Solution retrieval as well as post-processing and writing data to file :

get_total_solution(const BlockVector<double> &solution_delta) const;
void
output_results() const;

Finally, some member variables that describe the current state: A collection of the parameters used to describe the problem setup...

const Parameters::AllParameters &parameters;

...the volume of the reference and current configurations...

double vol_reference;
double vol_current;

...and description of the geometry on which the problem is solved:

Triangulation<dim> triangulation;

Also, keep track of the current time and the time spent evaluating certain functions

Time time;

A storage object for quadrature point information. As opposed to step-18, deal.II's native quadrature point data manager is employed here.

A description of the finite-element system including the displacement polynomial degree, the degree-of-freedom handler, number of DoFs per cell and the extractor objects used to retrieve information from the solution vectors:

const unsigned int degree;
const FESystem<dim> fe;
DoFHandler<dim> dof_handler_ref;
const unsigned int dofs_per_cell;

Description of how the block-system is arranged. There is just 1 block, that contains a vector DOF $$\mathbf{u}$$. There are two reasons that we retain the block system in this problem. The first is pure laziness to perform further modifications to the code from which this work originated. The second is that a block system would typically necessary when extending this code to multiphysics problems.

static const unsigned int n_blocks = 1;
static const unsigned int n_components = dim;
static const unsigned int first_u_component = 0;
enum
{
u_dof = 0
};
std::vector<types::global_dof_index> dofs_per_block;

Rules for Gauss-quadrature on both the cell and faces. The number of quadrature points on both cells and faces is recorded.

const QGauss<dim> qf_cell;
const QGauss<dim - 1> qf_face;
const unsigned int n_q_points;
const unsigned int n_q_points_f;

Objects that store the converged solution and right-hand side vectors, as well as the tangent matrix. There is a ConstraintMatrix object used to keep track of constraints. We make use of a sparsity pattern designed for a block system.

ConstraintMatrix constraints;
BlockSparsityPattern sparsity_pattern;

Then define a number of variables to store norms and update norms and normalisation factors.

struct Errors
{
Errors()
:
norm(1.0), u(1.0)
{}
void reset()
{
norm = 1.0;
u = 1.0;
}
void normalise(const Errors &rhs)
{
if (rhs.norm != 0.0)
norm /= rhs.norm;
if (rhs.u != 0.0)
u /= rhs.u;
}
double norm, u;
};
Errors error_residual, error_residual_0, error_residual_norm, error_update,
error_update_0, error_update_norm;

Methods to calculate error measures

void
get_error_residual(Errors &error_residual);
void
get_error_update(const BlockVector<double> &newton_update,
Errors &error_update);

Print information to screen in a pleasing way...

static
void
void
print_conv_footer();
void
print_vertical_tip_displacement();
};

### Implementation of the Solid class

#### Public interface

We initialise the Solid class using data extracted from the parameter file.

template <int dim,typename NumberType>
Solid<dim,NumberType>::Solid(const Parameters::AllParameters &parameters)
:
parameters(parameters),
vol_reference (0.0),
vol_current (0.0),
triangulation(Triangulation<dim>::maximum_smoothing),
time(parameters.end_time, parameters.delta_t),
timer(std::cout,
TimerOutput::summary,
TimerOutput::wall_times),
degree(parameters.poly_degree),

The Finite Element System is composed of dim continuous displacement DOFs.

fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
dof_handler_ref(triangulation),
dofs_per_cell (fe.dofs_per_cell),
u_fe(first_u_component),
dofs_per_block(n_blocks),
n_q_points (qf_cell.size()),
n_q_points_f (qf_face.size())
{
}

The class destructor simply clears the data held by the DOFHandler

template <int dim,typename NumberType>
Solid<dim,NumberType>::~Solid()
{
dof_handler_ref.clear();
}

In solving the quasi-static problem, the time becomes a loading parameter, i.e. we increasing the loading linearly with time, making the two concepts interchangeable. We choose to increment time linearly using a constant time step size.

We start the function with preprocessing, and then output the initial grid before starting the simulation proper with the first time (and loading) increment.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::run()
{
make_grid();
system_setup();
output_results();
time.increment();

We then declare the incremental solution update $$\varDelta \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}$$ and start the loop over the time domain.

At the beginning, we reset the solution update for this time step...

BlockVector<double> solution_delta(dofs_per_block);
while (time.current() <= time.end())
{
solution_delta = 0.0;

...solve the current time step and update total solution vector $$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} + \varDelta \mathbf{\Xi}$$...

solve_nonlinear_timestep(solution_delta);
solution_n += solution_delta;

...and plot the results before moving on happily to the next time step:

output_results();
time.increment();
}

Lastly, we print the vertical tip displacement of the Cook cantilever after the full load is applied

print_vertical_tip_displacement();
}

### Private interface

#### Solid::make_grid

On to the first of the private member functions. Here we create the triangulation of the domain, for which we choose a scaled an anisotripically discretised rectangle which is subsequently transformed into the correct of the Cook cantilever. Each relevant boundary face is then given a boundary ID number.

We then determine the volume of the reference configuration and print it for comparison.

template <int dim>
Point<dim> grid_y_transform (const Point<dim> &pt_in)
{
const double &x = pt_in[0];
const double &y = pt_in[1];
const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam
const double y_lower = 0.0 + (44.0/48.0)*x; // Line defining lower edge of beam
const double theta = y/44.0; // Fraction of height along left side of beam
const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation
Point<dim> pt_out = pt_in;
pt_out[1] = y_transform;
return pt_out;
}
template <int dim,typename NumberType>
void Solid<dim,NumberType>::make_grid()
{

Divide the beam, but only along the x- and y-coordinate directions

std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);

Only allow one element through the thickness (modelling a plane strain condition)

if (dim == 3)
repetitions[dim-1] = 1;
const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
const Point<dim> top_right = (dim == 3 ? Point<dim>(48.0, 44.0, 0.5) : Point<dim>(48.0, 44.0));
repetitions,
bottom_left,
top_right);

Since we wish to apply a Neumann BC to the right-hand surface, we must find the cell faces in this part of the domain and mark them with a distinct boundary ID number. The faces we are looking for are on the +x surface and will get boundary ID 11. Dirichlet boundaries exist on the left-hand face of the beam (this fixed boundary will get ID 1) and on the +Z and -Z faces (which correspond to ID 2 and we will use to impose the plane strain condition)

const double tol_boundary = 1e-6;
triangulation.begin_active(), endc = triangulation.end();
for (; cell != endc; ++cell)
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() == true)
{
if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
cell->face(face)->set_boundary_id(1); // -X faces
else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
cell->face(face)->set_boundary_id(11); // +X faces
else if (std::abs(std::abs(cell->face(face)->center()[0]) - 0.5) < tol_boundary)
cell->face(face)->set_boundary_id(2); // +Z and -Z faces
}

Transform the hyper-rectangle into the beam shape

GridTools::transform(&grid_y_transform<dim>, triangulation);
GridTools::scale(parameters.scale, triangulation);
vol_reference = GridTools::volume(triangulation);
vol_current = vol_reference;
std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
}

#### Solid::system_setup

Next we describe how the FE system is setup. We first determine the number of components per block. Since the displacement is a vector component, the first dim components belong to it.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::system_setup()
{
timer.enter_subsection("Setup system");
std::vector<unsigned int> block_component(n_components, u_dof); // Displacement

The DOF handler is then initialised and we renumber the grid in an efficient manner. We also record the number of DOFs per block.

dof_handler_ref.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler_ref, block_component);
DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
block_component);
std::cout << "Triangulation:"
<< "\n\t Number of active cells: " << triangulation.n_active_cells()
<< "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
<< std::endl;

Setup the sparsity pattern and tangent matrix

tangent_matrix.clear();
{
const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
BlockDynamicSparsityPattern csp(n_blocks, n_blocks);
csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
csp.collect_sizes();

Naturally, for a one-field vector-valued problem, all of the components of the system are coupled.

Table<2, DoFTools::Coupling> coupling(n_components, n_components);
for (unsigned int ii = 0; ii < n_components; ++ii)
for (unsigned int jj = 0; jj < n_components; ++jj)
coupling[ii][jj] = DoFTools::always;
coupling,
csp,
constraints,
false);
sparsity_pattern.copy_from(csp);
}
tangent_matrix.reinit(sparsity_pattern);

We then set up storage vectors

system_rhs.reinit(dofs_per_block);
system_rhs.collect_sizes();
solution_n.reinit(dofs_per_block);
solution_n.collect_sizes();

...and finally set up the quadrature point history:

setup_qph();
}

#### Solid::setup_qph

The method used to store quadrature information is already described in step-18 and step-44. Here we implement a similar setup for a SMP machine.

Firstly the actual QPH data objects are created. This must be done only once the grid is refined to its finest level.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::setup_qph()
{
std::cout << " Setting up quadrature point data..." << std::endl;
triangulation.end(),
n_q_points);

Next we setup the initial quadrature point data. Note that when the quadrature point data is retrieved, it is returned as a vector of smart pointers.

triangulation.begin_active(); cell != triangulation.end(); ++cell)
{
const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
Assert(lqph.size() == n_q_points, ExcInternalError());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
lqph[q_point]->setup_lqp(parameters);
}
}

#### Solid::solve_nonlinear_timestep

The next function is the driver method for the Newton-Raphson scheme. At its top we create a new vector to store the current Newton update step, reset the error storage objects and print solver header.

template <int dim,typename NumberType>
void
Solid<dim,NumberType>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
{
std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
<< time.current() << "s" << std::endl;
BlockVector<double> newton_update(dofs_per_block);
error_residual.reset();
error_residual_0.reset();
error_residual_norm.reset();
error_update.reset();
error_update_0.reset();
error_update_norm.reset();

We now perform a number of Newton iterations to iteratively solve the nonlinear problem. Since the problem is fully nonlinear and we are using a full Newton method, the data stored in the tangent matrix and right-hand side vector is not reusable and must be cleared at each Newton step. We then initially build the right-hand side vector to check for convergence (and store this value in the first iteration). The unconstrained DOFs of the rhs vector hold the out-of-balance forces. The building is done before assembling the system matrix as the latter is an expensive operation and we can potentially avoid an extra assembly process by not assembling the tangent matrix when convergence is attained.

unsigned int newton_iteration = 0;
for (; newton_iteration < parameters.max_iterations_NR;
++newton_iteration)
{
std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;

If we have decided that we want to continue with the iteration, we assemble the tangent, make and impose the Dirichlet constraints, and do the solve of the linearized system:

make_constraints(newton_iteration);
assemble_system(solution_delta);
get_error_residual(error_residual);
if (newton_iteration == 0)
error_residual_0 = error_residual;

We can now determine the normalised residual error and check for solution convergence:

error_residual_norm = error_residual;
error_residual_norm.normalise(error_residual_0);
if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
&& error_residual_norm.u <= parameters.tol_f)
{
std::cout << " CONVERGED! " << std::endl;
print_conv_footer();
break;
}
const std::pair<unsigned int, double>
lin_solver_output = solve_linear_system(newton_update);
get_error_update(newton_update, error_update);
if (newton_iteration == 0)
error_update_0 = error_update;

We can now determine the normalised Newton update error, and perform the actual update of the solution increment for the current time step, update all quadrature point information pertaining to this new displacement and stress state and continue iterating:

error_update_norm = error_update;
error_update_norm.normalise(error_update_0);
solution_delta += newton_update;
std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
<< std::scientific << lin_solver_output.first << " "
<< lin_solver_output.second << " " << error_residual_norm.norm
<< " " << error_residual_norm.u << " "
<< " " << error_update_norm.norm << " " << error_update_norm.u
<< " " << std::endl;
}

At the end, if it turns out that we have in fact done more iterations than the parameter file allowed, we raise an exception that can be caught in the main() function. The call AssertThrow(condition, exc_object) is in essence equivalent to if (!cond) throw exc_object; but the former form fills certain fields in the exception object that identify the location (filename and line number) where the exception was raised to make it simpler to identify where the problem happened.

AssertThrow (newton_iteration <= parameters.max_iterations_NR,
ExcMessage("No convergence in nonlinear solver!"));
}

This program prints out data in a nice table that is updated on a per-iteration basis. The next two functions set up the table header and footer:

template <int dim,typename NumberType>
{
static const unsigned int l_width = 87;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
std::cout << " SOLVER STEP "
<< " | LIN_IT LIN_RES RES_NORM "
<< " RES_U NU_NORM "
<< " NU_U " << std::endl;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
}
template <int dim,typename NumberType>
void Solid<dim,NumberType>::print_conv_footer()
{
static const unsigned int l_width = 87;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
std::cout << "Relative errors:" << std::endl
<< "Displacement:\t" << error_update.u / error_update_0.u << std::endl
<< "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
<< "v / V_0:\t" << vol_current << " / " << vol_reference
<< std::endl;
}

At the end we also output the result that can be compared to that found in the literature, namely the displacement at the upper right corner of the beam.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::print_vertical_tip_displacement()
{
static const unsigned int l_width = 87;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
Point<dim> soln_pt (48.0*parameters.scale,60.0*parameters.scale);
if (dim == 3)
soln_pt[2] = 0.5*parameters.scale;
double vertical_tip_displacement = 0.0;
double vertical_tip_displacement_check = 0.0;
dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
for (; cell != endc; ++cell)
{

if (cell->point_inside(soln_pt) == true)

for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
if (cell->vertex(v).distance(soln_pt) < 1e-6)
{

Extract y-component of solution at the given point This point is coindicent with a vertex, so we can extract it directly as we're using FE_Q finite elements that have support at the vertices

vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));

Sanity check using alternate method to extract the solution at the given point. To do this, we must create an FEValues instance to help us extract the solution value at the desired point

const MappingQ<dim> mapping (parameters.poly_degree);
const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
AssertThrow(soln_qrule.size() == 1, ExcInternalError());
FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
fe_values_soln.reinit(cell);

Extract y-component of solution at given point

std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
fe_values_soln[u_fe].get_function_values(solution_n,
soln_values);
vertical_tip_displacement_check = soln_values[0][u_dof+1];
break;
}
}
AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"))
std::cout << "Vertical tip displacement: " << vertical_tip_displacement
<< "\t Check: " << vertical_tip_displacement_check
<< std::endl;
}

#### Solid::get_error_residual

Determine the true residual error for the problem. That is, determine the error in the residual for the unconstrained degrees of freedom. Note that to do so, we need to ignore constrained DOFs by setting the residual in these vector components to zero.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
{
BlockVector<double> error_res(dofs_per_block);
for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
if (!constraints.is_constrained(i))
error_res(i) = system_rhs(i);
error_residual.norm = error_res.l2_norm();
error_residual.u = error_res.block(u_dof).l2_norm();
}

#### Solid::get_error_udpate

Determine the true Newton update error for the problem

template <int dim,typename NumberType>
void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
Errors &error_update)
{
BlockVector<double> error_ud(dofs_per_block);
for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
if (!constraints.is_constrained(i))
error_ud(i) = newton_update(i);
error_update.norm = error_ud.l2_norm();
error_update.u = error_ud.block(u_dof).l2_norm();
}

#### Solid::get_total_solution

This function provides the total solution, which is valid at any Newton step. This is required as, to reduce computational error, the total solution is only updated at the end of the timestep.

template <int dim,typename NumberType>
Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
{
BlockVector<double> solution_total(solution_n);
solution_total += solution_delta;
return solution_total;
}

#### Solid::assemble_system

template <int dim,typename NumberType>
struct Assembler_Base
{
virtual ~Assembler_Base() {}

Here we deal with the tangent matrix assembly structures. The PerTaskData object stores local contributions.

{
const Solid<dim,NumberType> *solid;
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
:
solid (solid),
cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
cell_rhs(solid->dofs_per_cell),
local_dof_indices(solid->dofs_per_cell)
{}
void reset()
{
cell_matrix = 0.0;
cell_rhs = 0.0;
}
};

On the other hand, the ScratchData object stores the larger objects such as the shape-function values array (Nx) and a shape function gradient and symmetric gradient vector which we will use during the assembly.

struct ScratchData_ASM
{
const BlockVector<double> &solution_total;
FEValues<dim> fe_values_ref;
FEFaceValues<dim> fe_face_values_ref;
ScratchData_ASM(const FiniteElement<dim> &fe_cell,
const QGauss<dim> &qf_cell,
const UpdateFlags uf_cell,
const QGauss<dim-1> & qf_face,
const UpdateFlags uf_face,
const BlockVector<double> &solution_total)
:
solution_total(solution_total),
fe_values_ref(fe_cell, qf_cell, uf_cell),
fe_face_values_ref(fe_cell, qf_face, uf_face),
std::vector<Tensor<2,dim,NumberType> >(fe_cell.dofs_per_cell)),
std::vector<SymmetricTensor<2,dim,NumberType> >
(fe_cell.dofs_per_cell))
{}
ScratchData_ASM(const ScratchData_ASM &rhs)
:
solution_total (rhs.solution_total),
fe_values_ref(rhs.fe_values_ref.get_fe(),
rhs.fe_values_ref.get_update_flags()),
fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
rhs.fe_face_values_ref.get_update_flags()),
{}
void reset()
{
const unsigned int n_q_points = fe_values_ref.get_quadrature().size();
const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell;
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
{
}
}
}
};

Of course, we still have to define how we assemble the tangent matrix contribution for a single cell.

void
assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM &scratch,
{

Due to the C++ specialization rules, we need one more level of indirection in order to define the assembly routine for all different number. The next function call is specialized for each NumberType, but to prevent having to specialize the whole class along with it we have inlined the definition of the other functions that are common to all implementations.

assemble_system_tangent_residual_one_cell(cell, scratch, data);
assemble_neumann_contribution_one_cell(cell, scratch, data);
}

This function adds the local contribution to the system matrix.

void
{
const ConstraintMatrix &constraints = data.solid->constraints;
BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;
BlockVector<double> &system_rhs = const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;
data.cell_matrix, data.cell_rhs,
data.local_dof_indices,
tangent_matrix, system_rhs);
}
protected:

This function needs to exist in the base class for Workstream to work with a reference to the base class.

virtual void
assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/ *cell* /,
ScratchData_ASM &/ *scratch* /,
{
}
void
assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM &scratch,
{

Aliases for data referenced from the Solid class

const unsigned int &n_q_points_f = data.solid->n_q_points_f;
const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
const Parameters::AllParameters &parameters = data.solid->parameters;
const Time &time = data.solid->time;
const FESystem<dim> &fe = data.solid->fe;
const unsigned int &u_dof = data.solid->u_dof;

Next we assemble the Neumann contribution. We first check to see it the cell face exists on a boundary on which a traction is applied and add the contribution if this is the case.

for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() == true
&& cell->face(face)->boundary_id() == 11)
{
scratch.fe_face_values_ref.reinit(cell, face);
for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
++f_q_point)
{

We specify the traction in reference configuration. For this problem, a defined total vertical force is applied in the reference configuration. The direction of the applied traction is assumed not to evolve with the deformation of the domain.

Note that the contributions to the right hand side vector we compute here only exist in the displacement components of the vector.

const double time_ramp = (time.current() / time.end());
const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area)
dir[1] = 1.0;
const Tensor<1, dim> traction = magnitude*dir;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int i_group =
fe.system_to_base_index(i).first.first;
if (i_group == u_dof)
{
const unsigned int component_i =
const double Ni =
scratch.fe_face_values_ref.shape_value(i,
f_q_point);
const double JxW = scratch.fe_face_values_ref.JxW(
f_q_point);
data.cell_rhs(i) += (Ni * traction[component_i])
* JxW;
}
}
}
}
}
};
template <int dim>
struct Assembler<dim,double> : Assembler_Base<dim,double>
{
typedef double NumberType;
using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
virtual ~Assembler() {}
virtual void
assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM &scratch,
{

Aliases for data referenced from the Solid class

const unsigned int &n_q_points = data.solid->n_q_points;
const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
const FESystem<dim> &fe = data.solid->fe;
const unsigned int &u_dof = data.solid->u_dof;
const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
data.reset();
scratch.reset();
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
const std::vector<std::shared_ptr<const PointHistory<dim,NumberType> > > lqph =
Assert(lqph.size() == n_q_points, ExcInternalError());

We first need to find the solution gradients at quadrature points inside the current cell and then we update each local QP using the displacement gradient:

Now we build the local cell stiffness matrix. Since the global and local system matrices are symmetric, we can exploit this property by building only the lower half of the local matrix and copying the values to the upper half.

In doing so, we first extract some configuration dependent variables from our QPH history objects for the current quadrature point.

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const NumberType det_F = determinant(F);
const Tensor<2,dim,NumberType> F_inv = invert(F);
Assert(det_F > NumberType(0.0), ExcInternalError());
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
const unsigned int k_group = fe.system_to_base_index(k).first.first;
if (k_group == u_dof)
{
}
else
Assert(k_group <= u_dof, ExcInternalError());
}
const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
const SymmetricTensor<4,dim,NumberType> Jc = lqph[q_point]->get_Jc(det_F,b_bar);
const Tensor<2,dim,NumberType> tau_ns (tau);

Next we define some aliases to make the assembly process easier to follow

const double JxW = scratch.fe_values_ref.JxW(q_point);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int component_i = fe.system_to_component_index(i).first;
const unsigned int i_group = fe.system_to_base_index(i).first.first;
if (i_group == u_dof)
data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
else
Assert(i_group <= u_dof, ExcInternalError());
for (unsigned int j = 0; j <= i; ++j)
{
const unsigned int component_j = fe.system_to_component_index(j).first;
const unsigned int j_group = fe.system_to_base_index(j).first.first;

This is the $$\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}$$ contribution. It comprises a material contribution, and a geometrical stress contribution which is only added along the local matrix diagonals:

if ((i_group == j_group) && (i_group == u_dof))
{
data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
if (component_i == component_j) // geometrical stress contribution
data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
}
else
Assert((i_group <= u_dof) && (j_group <= u_dof),
}
}
}

Finally, we need to copy the lower half of the local matrix into the upper half:

for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
data.cell_matrix(i, j) = data.cell_matrix(j, i);
}
};
template <int dim>
{
virtual ~Assembler() {}
virtual void
assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM &scratch,
{

Aliases for data referenced from the Solid class

const unsigned int &n_q_points = data.solid->n_q_points;
const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
const FESystem<dim> &fe = data.solid->fe;
const unsigned int &u_dof = data.solid->u_dof;
const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
data.reset();
scratch.reset();
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
Assert(lqph.size() == n_q_points, ExcInternalError());
const unsigned int n_independent_variables = data.local_dof_indices.size();
std::vector<double> local_dof_values(n_independent_variables);
cell->get_dof_values(scratch.solution_total,
local_dof_values.begin(),
local_dof_values.end());

We now retrieve a set of degree-of-freedom values that have the operations that are performed with them tracked.

for (unsigned int i=0; i<n_independent_variables; ++i)

Accumulate the residual value for each degree of freedom. Note: Its important that the vectors is initialised (zero'd) correctly.

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
const unsigned int k_group = fe.system_to_base_index(k).first.first;
if (k_group == u_dof)
{
}
else
Assert(k_group <= u_dof, ExcInternalError());
}

Next we define some position-dependent aliases, again to make the assembly process easier to follow.

const double JxW = scratch.fe_values_ref.JxW(q_point);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int i_group = fe.system_to_base_index(i).first.first;
if (i_group == u_dof)
else
Assert(i_group <= u_dof, ExcInternalError());
}
}
for (unsigned int I=0; I<n_independent_variables; ++I)
{
data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
for (unsigned int J=0; J<n_independent_variables; ++J)
{

Compute the gradients of the residual entry [forward-mode]

data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
}
}
}
};
template <int dim>
{
virtual ~Assembler() {}
virtual void
assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM &scratch,
{

Aliases for data referenced from the Solid class

const unsigned int &n_q_points = data.solid->n_q_points;
const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
data.reset();
scratch.reset();
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
Assert(lqph.size() == n_q_points, ExcInternalError());
const unsigned int n_independent_variables = data.local_dof_indices.size();
std::vector<double> local_dof_values(n_independent_variables);
cell->get_dof_values(scratch.solution_total,
local_dof_values.begin(),
local_dof_values.end());

We now retrieve a set of degree-of-freedom values that have the operations that are performed with them tracked.

for (unsigned int i=0; i<n_independent_variables; ++i)

Next we compute the total potential energy of the element. This is defined as follows: Total energy = (internal - external) energies Note: Its important that this value is initialised (zero'd) correctly.

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{

Next we define some position-dependent aliases, again to make the assembly process easier to follow.

const double JxW = scratch.fe_values_ref.JxW(q_point);

We extract the configuration-dependent material energy from our QPH history objects for the current quadrature point and integrate its contribution to increment the total cell energy.

}

Compute derivatives of reverse-mode AD variables

for (unsigned int I=0; I<n_independent_variables; ++I)
{

This computes the adjoint df/dX_{i} [reverse-mode]

data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
for (unsigned int J=0; J<n_independent_variables; ++J)
{

Compute the gradients of the residual entry [forward-mode]

data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
}
}
}
};
#endif

Since we use TBB for assembly, we simply setup a copy of the data structures required for the process and pass them, along with the memory addresses of the assembly functions to the WorkStream object for processing. Note that we must ensure that the matrix is reset before any assembly operations can occur.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
{
timer.enter_subsection("Assemble linear system");
std::cout << " ASM " << std::flush;
tangent_matrix = 0.0;
system_rhs = 0.0;
const UpdateFlags uf_face(update_values |
const BlockVector<double> solution_total(get_total_solution(solution_delta));
typename Assembler_Base<dim,NumberType>::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
Assembler<dim,NumberType> assembler;
WorkStream::run(dof_handler_ref.begin_active(),
dof_handler_ref.end(),
static_cast<Assembler_Base<dim,NumberType>&>(assembler),
&Assembler_Base<dim,NumberType>::assemble_system_one_cell,
&Assembler_Base<dim,NumberType>::copy_local_to_global_ASM,
scratch_data,
}

#### Solid::make_constraints

The constraints for this problem are simple to describe. However, since we are dealing with an iterative Newton method, it should be noted that any displacement constraints should only be specified at the zeroth iteration and subsequently no additional contributions are to be made since the constraints are already exactly satisfied.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::make_constraints(const int &it_nr)
{
std::cout << " CST " << std::flush;

Since the constraints are different at different Newton iterations, we need to clear the constraints matrix and completely rebuild it. However, after the first iteration, the constraints remain the same and we can simply skip the rebuilding step if we do not clear it.

if (it_nr > 1)
return;
constraints.clear();
const bool apply_dirichlet_bc = (it_nr == 0);

The boundary conditions for the indentation problem are as follows: On the -x, -y and -z faces (ID's 0,2,4) we set up a symmetry condition to allow only planar movement while the +x and +y faces (ID's 1,3) are traction free. In this contrived problem, part of the +z face (ID 5) is set to have no motion in the x- and y-component. Finally, as described earlier, the other part of the +z face has an the applied pressure but is also constrained in the x- and y-directions.

In the following, we will have to tell the function interpolation boundary values which components of the solution vector should be constrained (i.e., whether it's the x-, y-, z-displacements or combinations thereof). This is done using ComponentMask objects (see GlossComponentMask) which we can get from the finite element if we provide it with an extractor object for the component we wish to select. To this end we first set up such extractor objects and later use it when generating the relevant component masks:

Fixed left hand side of the beam

{
const int boundary_id = 1;
if (apply_dirichlet_bc == true)
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
else
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
}

Zero Z-displacement through thickness direction This corresponds to a plane strain condition being imposed on the beam

if (dim == 3)
{
const int boundary_id = 2;
const FEValuesExtractors::Scalar z_displacement(2);
if (apply_dirichlet_bc == true)
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
else
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
}
constraints.close();
}

#### Solid::solve_linear_system

As the system is composed of a single block, defining a solution scheme for the linear problem is straight-forward.

template <int dim,typename NumberType>
std::pair<unsigned int, double>
Solid<dim,NumberType>::solve_linear_system(BlockVector<double> &newton_update)
{
BlockVector<double> A(dofs_per_block);
BlockVector<double> B(dofs_per_block);
unsigned int lin_it = 0;
double lin_res = 0.0;

We solve for the incremental displacement $$d\mathbf{u}$$.

{
timer.enter_subsection("Linear solver");
std::cout << " SLV " << std::flush;
if (parameters.type_lin == "CG")
{
const int solver_its = tangent_matrix.block(u_dof, u_dof).m()
* parameters.max_iterations_lin;
const double tol_sol = parameters.tol_lin
* system_rhs.block(u_dof).l2_norm();
SolverControl solver_control(solver_its, tol_sol);
SolverCG<Vector<double> > solver_CG(solver_control, GVM);

We've chosen by default a SSOR preconditioner as it appears to provide the fastest solver convergence characteristics for this problem on a single-thread machine. However, for multicore computing, the Jacobi preconditioner which is multithreaded may converge quicker for larger linear systems.

PreconditionSelector<SparseMatrix<double>, Vector<double> >
preconditioner (parameters.preconditioner_type,
parameters.preconditioner_relaxation);
preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
newton_update.block(u_dof),
system_rhs.block(u_dof),
preconditioner);
lin_it = solver_control.last_step();
lin_res = solver_control.last_value();
}
else if (parameters.type_lin == "Direct")
{

Otherwise if the problem is small enough, a direct solver can be utilised.

A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
lin_it = 1;
lin_res = 0.0;
}
else
Assert (false, ExcMessage("Linear solver type not implemented"));
}

Now that we have the displacement update, distribute the constraints back to the Newton update:

constraints.distribute(newton_update);
return std::make_pair(lin_it, lin_res);
}

#### Solid::output_results

Here we present how the results are written to file to be viewed using ParaView or Visit. The method is similar to that shown in the tutorials so will not be discussed in detail.

template <int dim,typename NumberType>
void Solid<dim,NumberType>::output_results() const
{
DataOut<dim> data_out;
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(dim,
std::vector<std::string> solution_name(dim, "displacement");
data_out.attach_dof_handler(dof_handler_ref);
solution_name,
data_component_interpretation);

Since we are dealing with a large deformation problem, it would be nice to display the result on a displaced grid! The MappingQEulerian class linked with the DataOut class provides an interface through which this can be achieved without physically moving the grid points in the Triangulation object ourselves. We first need to copy the solution to a temporary vector and then create the Eulerian mapping. We also specify the polynomial degree to the DataOut object in order to produce a more refined output data set when higher order polynomials are used.

Vector<double> soln(solution_n.size());
for (unsigned int i = 0; i < soln.size(); ++i)
soln(i) = solution_n(i);
MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
data_out.build_patches(q_mapping, degree);
std::ostringstream filename;
filename << "solution-" << time.get_timestep() << ".vtk";
std::ofstream output(filename.str().c_str());
data_out.write_vtk(output);
}
}

### Main function

Lastly we provide the main driver function which appears no different to the other tutorials.

int main (int argc, char *argv[])
{
using namespace dealii;
using namespace Cook_Membrane;
const unsigned int dim = 2;
try
{
deallog.depth_console(0);
Parameters::AllParameters parameters("parameters.prm");
if (parameters.automatic_differentiation_order == 0)
{
std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;

Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
typedef double NumberType;
Solid<dim,NumberType> solid_3d(parameters);
solid_3d.run();
}
else if (parameters.automatic_differentiation_order == 1)
{
std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;

Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
Solid<dim,NumberType> solid_3d(parameters);
solid_3d.run();
}
else if (parameters.automatic_differentiation_order == 2)
{
std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;

Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
1);
Solid<dim,NumberType> solid_3d(parameters);
solid_3d.run();
}
#endif
else
{
AssertThrow(false,
ExcMessage("The selected assembly method is not supported. "
"You need deal.II 9.0 and Trilinos with the Sacado package "
"to enable assembly using automatic differentiation."));
}
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!"
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}


text:
# Modules

ModuleSource FileDescription
mod_arpack_typemod_arpack_type.f08

Module containing a dedicated type for the ARPACK solvers. All needed variables are defined and allocated here, along with sanity checks when setting certain values.

mod_assertmod_assert.f08

Module defining the assertion routine used by assert.fpp.

mod_atmosphere_curvesmod_atmosphere_curves.f08

Contains data for a realistic solar atmosphere model taken from Avrett, E. H., & Loeser, R. (2008). Models of the solar chromosphere and transition region from SUMER and HRTS observations: formation of the extreme-ultraviolet spectrum of hydrogen, carbon, and oxygen. ApJS, 175(1), 229. link.

mod_boundary_managermod_boundary_manager.f08
smod_essential_boundariessmod_essential_boundaries.f08
smod_natural_boundariessmod_natural_boundaries.f08
smod_natural_bounds_conductionsmod_natural_bounds_conduction.f08
smod_natural_bounds_flowsmod_natural_bounds_flow.f08
smod_natural_bounds_hallsmod_natural_bounds_hall.f08
smod_natural_bounds_regularsmod_natural_bounds_regular.f08
smod_natural_bounds_resistivesmod_natural_bounds_resistive.f08
smod_natural_bounds_viscositysmod_natural_bounds_viscosity.f08
mod_check_valuesmod_check_values.f08

This module contains various methods to check for small, NaN or negative values, equal values or inf values. Interfaces are provided for functionality with real and complex variables.

mod_cooling_curvesmod_cooling_curves.f08

All data of the different cooling curves is contained in this module, along with handling piecewise cooling curves.

mod_eigenfunctionsmod_eigenfunctions.f08

Main module responsible for eigenfunction management. Contains routines and interfaces to initialise and calculate the eigenfunctions and derived quantities.

smod_base_efssmod_base_efs.f08

This submodule initialises and calculates the base eigenfunctions, i.e. the ones corresponding to the basic state vector variables (rho, v1, v2, etc.)

smod_ef_operationssmod_ef_operations.f08

This submodule contains procedures and functions to be used when we are either assembling eigenfunctions based on the eigenvectors, or retransforming them.

smod_derived_efssmod_derived_efs.f08

This submodule calculates various quantities derived from the base eigenfunctions: - entropy S - $\nabla \cdot \mathbf{v}_1$ - all 3 components of $\nabla \times \mathbf{v}_1$ - $\nabla \cdot \mathbf{B}_1$ - all 3 components of $\mathbf{B}_1$ - all 3 components of $\nabla \times \mathbf{B}_1$

smod_derived_pp_efssmod_derived_pp_efs.f08

This submodule calculates quantities derived from the base eigenfunctions that are parallel or perpendicular to the background magnetic field: - parallel and perpendicular components of $\mathbf{B}_1$ - parallel and perpendicular components of $\nabla \times \mathbf{B}_1$ - parallel and perpendicular components of $\mathbf{v}_1$ - parallel and perpendicular components of $\nabla \times \mathbf{v}_1$

mod_equilibriummod_equilibrium.f08

Parent module governing all equilibrium types and submodules. This module contains all equilibrium types and the initial declarations of the module subroutines. Every equilibrium submodule extends this module, implementing one of the module subroutines declared here. All "main" equilibrium configurations are set in the submodules. The ones that depend on "main" arrays, like radiative cooling, are set here through calls to their respective modules.

This submodule defines a simple, adiabatic homogeneous medium in Cartesian geometry. The geometry can be overridden using the parfile.

smod_equil_constant_currentsmod_equil_constant_current.f08

This submodule defines an equilibrium in cylindrical geometry with a constant axial current. The geometry can be overridden using the parfile.

smod_equil_coronal_flux_tubesmod_equil_coronal_flux_tube.f08

This submodule defines a magnetic flux tube embedded in a uniform magnetic environment. In this case the flux tube is under coronal conditions $c_s < c_A < c_{Ae}$ where the subscript e denotes the outer region. More specifically the equilibrium is defined as $c_{Ae} = 5c_s, c_{se} = c_s/2, c_A = 2c_s$. The geometry can be overridden in the parfile, and is cylindrical by default for $r \in [0, 10]$.

smod_equil_couette_flowsmod_equil_couette_flow.f08

This submodule defines a steady plane Couette flow in a Cartesian geometry.

smod_equil_discrete_alfvensmod_equil_discrete_alfven.f08

This submodule defines an equilibrium in cylindrical geometry with an axial current profile ($\nu = 2$), modelling a solar coronal loop in which discrete Alfvén waves are present. The geometry can be overridden in the parfile.

smod_equil_flow_driven_instabilitiessmod_equil_flow_driven_instabilities.f08

This submodule defines flow driven instabilities in a Cartesian geometry. This equilibrium can not be called explicitly from the parfile, but rather acts as a "parent setup" for the Rayleigh-Taylor and Kelvin-Helmholtz submodules which use this specific kind of equilibrium but with different parameters. This submodule is called within its implicit children.

smod_equil_gold_hoylesmod_equil_gold_hoyle.f08

This submodule defines a Gold-Hoyle equilibrium in cylindrical geometry. This equilibrium configuration models a filament with a uniform twist such that all fieldlines perform an equal amount of turns around the cylinder axis. The geometry can be overridden in the parfile.

smod_equil_gravito_acousticsmod_equil_gravito_acoustic.f08

This submodule defines an equilibrium in Cartesian geometry with a stratified equilibrium profile, giving rise to gravito-acoustic waves. No magnetic fields are included, such that this treats the hydrodynamic regime. The geometry can be overridden using the parfile.

smod_equil_gravito_mhdsmod_equil_gravito_mhd.f08

This submodule defines an equilibrium in Cartesian geometry with a stratified equilibrium profile, giving rise to gravito-MHD waves. The geometry can be overridden using the parfile.

smod_equil_harris_sheetsmod_equil_harris_sheet.f08

This submodule defines a resistive equilibrium with tearing modes created by a Harris sheet.

smod_equil_interchange_modessmod_equil_interchange_modes.f08

This submodule defines an exponentially stratified medium in Cartesian geometry with a constant gravity term and magnetic shear. The geometry can be overridden in the parfile.

smod_equil_internal_kink_instabilitysmod_equil_internal_kink_instability.f08

This submodule defines internal kink modes in force-free magnetic fields. The geometry is cylindrical with parabolic density and velocity profiles, The geometry can be overridden in the parfile.

smod_equil_isothermal_atmospheresmod_equil_isothermal_atmosphere.f08

This submodule defines an equilibrium in Cartesian geometry with a stratified equilibrium profile, representing a solar magnetic atmosphere. The equilibrium is isothermal with a constant magnetic field. The scale height is given by The geometry is fixed to Cartesian, boundaries can be overridden using the parfile.

smod_equil_kelvin_helmholtz_cdsmod_equil_kelvin_helmholtz_cd.f08

This submodule defines an unperturbed magnetised jet model in cylindrical geometry, giving rise to Kelvin-Helmholtz and current-driven instabilities. The geometry is fixed for this problem; the cylinder wall is dependent on the equilibrium parameters and is given by 2rj.

smod_equil_KHIsmod_equil_KHI.f08

This submodule defines Kelvin-Helmholtz instabilities in Cartesian geometry. This equilibrium is a specific case of the flow driven instabilities. No magnetic fields are considered in this case (pure HD).

smod_equil_magnetothermal_instabilitiessmod_equil_magnetothermal_instabilities.f08

This submodule defines an equilibrium containing magnetothermal instabilities in a cylindrical geometry. The geometry can be overridden in the parfile.

smod_equil_MRIsmod_equil_MRI_accretion.f08

This submodule defines magneto-rotational instabilities in an accretion disk. Due to the special nature of this equilibrium x_start is hardcoded to one and can not be overridden in the parfile, the same goes for the geometry which is hardcoded to 'cylindrical'. The outer edge can be chosen freely. This equilibrium is chosen in such a way that the angular rotation is of order unity, implying Keplerian rotation. The thin-disk approximation is valid with small magnetic fields, but still large enough to yield magneto-rotational instabilities. Gravity is assumed to go like $g \thicksim 1/r^2$.

smod_equil_photospheric_flux_tubesmod_equil_photospheric_flux_tube.f08

This submodule defines a magnetic flux tube embedded in a uniform magnetic environment. In this case the flux tube is under photospheric conditions $c_{Ae} < c_s < c_{se} < c_A$ where the subscript e denotes the outer region. More specifically the equilibrium is defined as $c_{Ae} = c_s/2, c_{se} = 3c_s/2, c_A = 2c_s$. The geometry can be overridden in the parfile, and is cylindrical by default for $r \in [0, 10]$.

smod_equil_resistive_homosmod_equil_resistive_homo.f08

This submodule defines a simple, homogeneous medium in Cartesian geometry with a constant resistivity value. The geometry can be overridden using the parfile.

smod_equil_resistive_tearingsmod_equil_resistive_tearing.f08

This submodule defines an equilibrium in Cartesian geometry with a constant resistivity value. Parameters are taken in such a way as to allow for resistive tearing modes. The geometry can be overridden using the parfile.

smod_equil_resistive_tearing_flowsmod_equil_resistive_tearing_flow.f08

This submodule defines an equilibrium in Cartesian geometry with a flow profile and a constant resistivity value. Parameters are taken in such a way as to allow for resistive tearing modes. The geometry can be overridden using the parfile.

smod_equil_resonant_absorptionsmod_equil_resonant_absorption.f08

This submodule defines an inhomogeneous medium in Cartesian geometry with a constant resistivity value. Two (constant) density profiles are defined which are connected by a sine profile, representing the interface inbetween. This density profile allows for resonant absorption, the geometry can be overridden in the parfile.

smod_equil_rotating_plasma_cylindersmod_equil_rotating_plasma_cylinder.f08

This submodule defines a cylindrical equilibrium, resembling a rotating plasma cylinder. The geometry can be overridden using the parfile.

smod_equil_RTIsmod_equil_RTI.f08

This submodule defines Rayleigh-Taylor instabilities in Cartesian geometry. This equilibrium is a specific case of the flow driven instabilities.

smod_equil_RTI_KHIsmod_equil_RTI_KHI.f08

This submodule defines Rayleigh-Taylor and Kelvin-Helmholtz instabilities in Cartesian geometry. This equilibrium is a specific case of the flow driven instabilities.

smod_equil_RTI_theta_pinchsmod_equil_RTI_theta_pinch.f08

This submodule defines Rayleigh-Taylor instabilities in rotating theta pinches. The straight cylinder approximation is used with a constant angular frequency. Density and pressure profiles decrease over the domain, with a uni-directional increasing magnetic field profile. Mode numbers $k = 0$ correspond to HD Rayleigh-Taylor instabilities, while $k \neq 0$ represent MHD RTIs. The geometry is hardcoded to 'cylindrical', the domain is forced to $0 <= r <= 1$ through division by x_end.

smod_equil_suydam_clustersmod_equil_suydam_cluster.f08

This submodule defines a cylindrical equilibrium in which a Suydam surface is present, such that this gives rises to Suydam cluster modes. The geometry can be overriden using the parfile.

smod_equil_taylor_couettesmod_equil_taylor_couette.f08

This submodule defines a steady Taylor-Couette flow in a cylindrical geometry where a fluid is confined between two (rotating) coaxial cylinders (without a magnetic field).

smod_equil_tc_pinchsmod_equil_tc_pinch.f08

This submodule defines a steady Taylor-Couette flow in a cylindrical geometry where a plasma is confined between two (rotating) coaxial cylinders with an azimuthal magnetic field and constant resistivity.

smod_user_definedsmod_user_defined.f08

Submodule for user-defined equilibria. Look at the examples in the equilibria subdirectory or consult the website for more information.

mod_equilibrium_paramsmod_equilibrium_params.f08

Module containing all equilibrium parameters. All parameters used in the equilibrium submodules are defined here for convenience, including the wave numbers $k_2$ and $k_3$. All of these values are NaN initially, such that variables that are not properly set propagate their value and are easy to spot in follow-up checks.

mod_exceptionsmod_exceptions.f08

Module to explicitly handle exceptions. Depending on the application at hand we override what happens when an exception is raised, which is useful for testing purposes (no error stop if we expect something to fail). Loosely based on an example given in https://github.com/Goddard-Fortran-Ecosystem/pFUnit_demos/blob/main/Basic/src/throw.F90

mod_global_variablesmod_global_variables.f08

All global variables are defined in this module, and can be accessed throughout the entire code. Most variables are first initialised to a default value, and are properly set either in the parfile or in an equilibrium submodule.

mod_gridmod_grid.f08

Module containing all grid-related things. Contains subroutines to create the base grid, Gaussian grid and scale factors. Also handles mesh accumulation. An integral of $f(x)$ in $[a, b]$ can be approximated with where $w_i$ and $x_i$ are the weights and nodes of the Gaussian quadrature. The Gaussian grid is hence set up in every interval $[a, b]$ across the four nodes $j$ as

mod_hallmod_hall.f08

!> Module containing Hall-related routines, calculates and sets the Hall and inertia term factors based on the specified profiles.

mod_inputmod_input.f08

Module to handle parfile reading. Contains subroutines to retrieve the parfile based on the commandline arguments and to read the parfile, setting the global variables.

mod_inspectionsmod_inspections.f08

Module to inspect if certain conditions are fulfilled by doing additional sanity checks on the equilibrium configuration. For cylindrical geometries we check if $k_2$ is an integer and if the on-axis values obey regularity conditions. Equilibrium balance for both the Cartesian and cylindrical cases is checked.

mod_integrationmod_integration.f08

Module responsible for integration of differential equations, useful when setting equilibria or integrating the equilibrium equation. Contains subroutines to numerically solve the following systems of differential equations: These are solved using a fifth-order Runge-Kutta method.

mod_interpolationmod_interpolation.f08

Module responsible for table interpolations and array lookups. Contains subroutines for table interpolations, numerical derivatives of arrays and lookup functions. Subroutines are loosely based on routines implemented in the MPI-AMRVAC code.

mod_loggingmod_logging.f08

Main handler for console print statements. The level of information printed to the console depends on the corresponding global variable called logging_level defined in the parfile.

mod_make_subblockmod_make_subblock.f08

Module to build a single quadblock by creating the four comprising subblocks. This is done for every grid interval, for specific integral elements and spline functions.

mod_matrix_creationmod_matrix_creation.f08

Assembly of the finite element matrices $\mathcal{A}$ and $\mathcal{B}$ in the matrix eigenvalue problem $\mathcal{A}\textbf{X} = \omega\mathcal{B}\textbf{X}$. Matrix creation proceeds by first assembling a quadblock at every grid interval using the integrals and finite element basis functions, after which this quadblock is shifted along the main diagonal.

mod_matrix_managermod_matrix_manager.f08
smod_conduction_matrixsmod_conduction_matrix.f08
smod_cooling_matrixsmod_cooling_matrix.f08
smod_flow_matrixsmod_flow_matrix.f08
smod_hall_matrixsmod_hall_matrix.f08
smod_regular_matrixsmod_regular_matrix.f08
smod_resistive_matrixsmod_resistive_matrix.f08
smod_viscosity_matrixsmod_viscosity_matrix.f08
mod_matrix_operationsmod_matrix_operations.f08

Module containing various matrix operations (inverting, multiplying, etc). Does calls to relevant BLAS or LAPACK routines.

mod_matrix_shortcutsmod_matrix_shortcuts.f08
mod_outputmod_output.f08

This module contains all routines for file opening and file writing.

mod_paintingmod_painting.f08

This module handles formatting of terminal-printed strings. Contains subroutines to colourise strings.

mod_physical_constantsmod_physical_constants.f08

All physical constants used in the code are defined in this module. We include values both in SI units and in cgs units for convenience. All values are taken from the NRL Plasma Formulary.

Module containing radiative cooling-related routines. This module is responsible for initialising the radiative cooling variables and a correct handling of the cooling curves. If an interpolated cooling curve is selected this module calls the interpolation module to create one.

mod_resistivitymod_resistivity.f08

Module containing resistivity-related routines, calculates and sets the resistivity values based on the equilibrium configuration.

mod_solar_atmospheremod_solar_atmosphere.f08

Module to set a realistic solar atmosphere, using tabulated density and temperature profiles (see mod_atmosphere_curves), in Cartesian geometries only.

mod_solversmod_solvers.f08

Parent module for everything solver-related. Interfaces to the different submodules are defined here, and the solve_evp routine calls the correct solver based on parfile settings.

smod_arnoldismod_arnoldi.f08

Submodule containing the implementation of the various ARPACK solvers. Preprocessor directives are used here, if ARPACK was not found during compilation most of this submodule will not be compiled to avoid unknown calls to methods. We first cycle through the possible modes, set up the relevant matrix operators and call the corresponding methods. Then the reverse communication interface znaupd is called, which will keep iterating until either everything is converged or the maximum number of iterations maxiter is reached. Then eigenvalues are extracted using zneupd, eigenvectors are calculated and the residual of each eigenvalue is computed which should be small for converged values.

smod_qr_invertsmod_qr_invert.f08

Submodule containing the implementation of the QR-invert algorithm. The original problem is written as a standard eigenvalue problem through where $\mathcal{B}$ is always properly conditioned. Various calls to mod_matrix_operations are done to invert the B-matrix and to do matrix multiplications. Eventually a call to LAPACK's zgeev routine is done to obtain all eigenvalues and eigenvectors.

smod_qz_directsmod_qz_direct.f08

Submodule containing the implementation of the QZ-direct algorithm. We keep the general form of the eigenvalue problem and solve this directly by calling LAPACK's zggev routine.

mod_spline_functionsmod_spline_functions.f08

Module that calculates the finite element basis functions. The different basis functions and their derivatives used throughout the code are calculated here. All routines defined in this module simply return the basis functions for a specific point r in the interval (rj_lo, rj_hi).

mod_thermal_conductionmod_thermal_conduction.f08

This module is responsible for calculating and setting the thermal conduction values based on the equilibrium configuration.

mod_timingmod_timing.f08

Module to provide timing facilities.

mod_typesmod_types.f08

Module containing the different types used in the code. All types are defined here, this includes all types where the equilibrium arrays are set as attributes, as well as a type to handle eigenfunctions.

mod_unitsmod_units.f08

This module handles checking and setting all unit normalisations. When setting the normalisations, one should always specify a unit length and a unit magneticfield. The normalisations use a reference value of 2 for plasma beta as is common practice, such that we can write $\beta/2 = \mu * p / B^2 = 1$. The pressure unit is then calculated as Then, if a unit density is specified this fixes the unit temperature and vice-versa, through use of the ideal gas law The other normalisations are then fixed, and are set through Where the latter two denote the luminosity and thermal conduction normalisations.

format: plain


text:
# DFT+U+V

This tutorial aims at explaining how to perform DFT+U+V calculations in Octopus. This correspond to adding not only an on-site Hubbard interaction U, but also an intersite interaction V. As a prototypical example, we will consider bulk silicon.

The DFT+U+V method, as well as its performances and implementation details are discussed in Ref.1

### Input

The input file we will use is the following one:

CalculationMode = gs
PeriodicDimensions = 3

Spacing = 0.5

%LatticeVectors
0.0 | 0.5 | 0.5
0.5 | 0.0 | 0.5
0.5 | 0.5 | 0.0
%

a = 10.18
%LatticeParameters
a | a | a
%

%ReducedCoordinates
"Si" | 0.0 | 0.0 | 0.0
"Si" | 1/4 | 1/4 | 1/4
%

nk = 2
%KPointsGrid
nk |  nk |  nk
0.5 | 0.5 | 0.5
0.5 | 0.0 | 0.0
0.0 | 0.5 | 0.0
0.0 | 0.0 | 0.5
%

ExtraStates = 4

%KPointsPath
10 |  10 |  15
0.5 | 0.0 | 0.0  # L point
0.0 | 0.0 | 0.0  # Gamma point
0.0 | 0.5 | 0.5  # X point
1.0 | 1.0 | 1.0  # Another Gamma point
%

DFTULevel = dft_u_acbn0
AOLoewdin = yes
UseAllAtomicOrbitals = yes
SkipSOrbitals = no

ACBN0IntersiteInteraction = yes
ACBN0IntersiteCutoff = 7

ExperimentalFeatures = yes

Most of the variable are the same as in the introduction tutorial on solids [Getting started with periodic systems](../Getting started with periodic systems), and the DFT+U related variables are explained in the tutorial DFT+U. For convenience, we are doing here a ground-state calculation with a k-point grid and a k-point path. This allows to get the bandstructure as a direct output of a ground-state calculation. However, the unoccupied states might not be converged with this approach. To perform a proper band-structure calculation, see the tutorial [Getting started with periodic systems](../Getting started with periodic systems).

Compared to a more conventional DFT+U calculation, we note few differences here:

• AOLoewdin = yes: We want to perform a Löwdin orthonormalization of the atomic orbitals. This is very important to avoid double counting between different orbitals.
• UseAllAtomicOrbitals = yes and SkipSOrbitals = no: Here we want to use all the atomic orbitals available in the pseudopotentials. In this case, the s orbitals are important to obtain a meaningful orthonormalization. Note that in Si, we indeed want to use all atomic orbitals. But in some strongly correlated materials, one might as well only use one set of atomic orbital using the Species block, as done for NiO in the DFT+U tutorial.
• ACBN0IntersiteInteraction = yes: We are asking for the intersite interaction (+V). The number of neighbors for which we are including the intersite interaction is given by a cutoff in real space, specified by the variable ACBN0IntersiteCutoff. Here we are using 7 Bohr, which gives us the first nearest neighbor only. In practice, this value needs to be converged, or needs to be chosen if one wants to include only certain shells of neighbors.

### Output

After running Octopus using the above input, we can look at the output. In the static/info file, we find the direct and indirect bandgap

Direct gap at ik=    2 of  0.0304 H
Indirect gap between ik=    8 and ik=    2 of  0.0177 H

Direct gap at ik=    6 of  0.1346 H
Indirect gap between ik=    6 and ik=   24 of  0.1148 H

As a comparison, without adding the +U+V, one finds at the LDA level a much smaller bandgap

Direct gap at ik=   14 of  0.1346 H
Indirect gap between ik=    2 and ik=   11 of  0.1148 H

Direct gap at ik=    1 of  0.1041 H
Indirect gap between ik=    3 and ik=    9 of  0.0810 H

The effect of bandgap opening can be understood more precisely from the bandstructure. As one can see, the conduction bands have been rigidly shifted toward higher energies, as expected from more advance calculations such as hybrid functionals or GW calculations.

Band structure of bulk silicon from DFT+U+V. The zero of energy has been shifted to the maximum of the occupied bands.

## References

1. N. Tancogne-Dejean, and A. Rubio, Parameter-free hybridlike functional based on an extended Hubbard model: DFT+U+V, Physical Review B 102 155117 (2020); ↩