Finite elasticity - gravity loading of cantilever beam

Introduction

This tutorial demonstrates how to setup and solve a nonlinear continuum mechanics problem using OpenCMISS-Iron in Python. For the purpose of this tutorial, we will be solving a 3D finite elasticity problem on an isotropic, homogenous rectangular beam under gravity load with the bottom face fixed.

OpenCMISS-Iron info: - Python script file - /opencmiss-build/opencmiss/install/virtual_environments/oclibs_venv_py38_release/lib/python_._/site-packages/opencmiss/iron/iron.py - API Documentation - http://opencmiss.org/documentation/apidoc/iron/latest/python/classiron_1_1_nodes.html

Learning outcomes

At the end of this tutorial you will:

• Know the steps involved in setting up and solving a nonlinear finite elasticity simulation with OpenCMISS-Iron.
• Know how different boundary conditions and mechanical loads can be applied to a geometry.
• Know how information about the deformation can be extracted for analysis.

Problem summary

The finite elasticity stress equilibrium equation

In this example we are solving the stress equilibrium equation for nonlinear finite elasticity. The stress equilibrium equation represents the principle of linear momentum as follows:

\[\displaystyle \nabla \sigma + \rho \mathbf{b}-\rho \mathbf {a}=0. \qquad \text{in} \qquad \Omega\]

where \(\sigma\) is the Cauchy stress tensor, \(\mathbf{b}\) is the body force vector and \(\mathbf{a}\) is the vector representing the acceleration due to any unbalanced forces.

The boundary equations for the stress equilibrium equation are partitioned into the Dirichlet boundary conditions representing a fixed displacement \(u_d\) over the boundary \(\Gamma_d\), and the Neumann boundary conditions representing the traction forces \(\mathbf{t}\) applied on the boundary \(\Gamma_t\) along the normal \(\mathbf{n}\) as follows:

\[\begin{split}\begin{aligned} \displaystyle u &= u_d \quad &\text{on} \quad \Gamma_d \\ \displaystyle \sigma^T \mathbf{n} &= \mathbf{t} \quad &\text{on} \quad \Gamma_t \end{aligned}\end{split}\]

In this example we will solve the stress equilibrium equation for a incompressible material, which we shall perform using the incompressibility constraint:

\[\displaystyle I_3-1 = 0\]

We will also use a Mooney Rivlin constitutive equation to describe the behaviour of the material.

Solution variables

The problem we are about to solve includes 4 dependent variables. The first three variables represent each of the 3D coordinates \((x,y,z)\) of the mesh nodes in the deformed state. The incompressibility constraint is satisfied using lagrange multipliers, which are represented as a scalar variable, \(p\). This variable is often referred to as the hydrostatic pressure.

Constitutive relationship

Another important set of equations that need to be included when solving a mechanics problem are the stress-strain relationships or constitutive relationships. There are a set of constitutive equations that have already been included within the OpenCMISS-Iron library as EquationSet subtypes. We will demonstrate how these constants are used to incorporate the constitutive equation in the simulation when setting up the equation set.

Setup

Loading the OpenCMISS-Iron library

In order to use OpenCMISS-Iron we have to first import the opencmiss.iron module from the OpenCMISS-Iron package.

import numpy as np

# Initialise OpenCMISS-Iron.
from opencmiss.iron import iron

Set up parameters

We first specify a set of variables that will be used throughout the tutorial. This is a common practice among experienced users to set these variables up at the start because we can then change these easily and re-run the rest of the code as required.

# Set constants
X, Y, Z = (1, 2, 3)

# Specify the number of local element directions.
dimensions = 3
number_of_xi = dimensions

# Specify the number of elements along each element direction.
# Note that the number of elements in each direction were set to create a mesh with uniform square elements.
number_global_x_elements = 5
number_global_y_elements = 5
number_global_z_elements = 5

# Set dimensions of the cube (in mm).
width = 40
length = 40
height = 90

interpolation_type = 1 # Linear Lagrange interpolation.
number_of_gauss_per_xi = 2 # Gauss points along each local element coordinate direction (xi).
use_pressure_basis = False
number_of_load_increments = 12

NOTE: If we have more than 2000 elements in each of the three directions, it exceeds computer memory and crashes. Once the RAM storage is exceeded, it uses the Swap space located on hdd (very slow) and puts us in an endless loop of copying to and receiving data from the RAM. Don’t do this on HPC, it will obliterate other people’s work.

# Version and derivative numbers for this specific problem as constants.
VersionNumber = 1
DerivativeNumber = 1

• interpolation_type is an integer variable that can be changed to choose one of the nine basis interpolation types defined in OpenCMISS-Iron Basis Interpolation Specifications Constants
• number_of_gauss_per_xi is the number of Gauss points used along a given local element direction for integrating the equilibrium equations. Note that there is a minimum number of gauss points required for a given interpolation scheme (e.g. linear Lagrange interpolation requires at least 2 Gauss points along each local element direction). Using less than the minimum number will result in the solver not converging.
• use_pressure_basis is a boolean variable that we can set to true to include the incompressibility constraint or set to false if the material we want to represent in the model is nearly incompressible or compressible.
• number_of_load_increments sets the number of load steps to be taken to solve the nonlinear mechanics problem. This concept is briefly explained below.

The weak form finite element matrix equations for the stress equilibrium equations above have been implemented in the OpenCMISS-Iron libraries. Numerical treatment of the equations will show that the equations are nonlinear and, specifically, the determination of the deformed configuration coordinates turn out to be like a root finding exercise. We need to find the deformed coordinates \(\mathbf{x}\) such that \(f(x)=0\). These equations are therefore solved using a nonlinear iterative solver, Newton Raphson technique to be exact. Nonlinear solvers require an initial guess at the root of the equation and if we are far away from the actual solution then the solvers can diverge and give spurious \(x\) values. Therefore, in most simulations it is best to split an entire mechanical load into lots of smaller loads. This is what number_of_load_increments allows us to do. You will see it come to use in the control loop section of the tutorial below.

The next section describes how we can interact with the OpenCMISS-Iron library through it’s object-oriented API to create and solve the mechanics problem.

NOTE: If the cell pops up with an unexplained CMFE error in Jupyter, please see the running Powershell terminal for more info.

Step by step guide

1. Creating a coordinate system

First we construct a coordinate system that will be used to describe the geometry in our problem. The 3D geometry will exist in a 3D space, so we need a 3D coordinate system.

# Create a 3D rectangular cartesian coordinate system.
coordinate_system_user_number = 1
coordinate_system = iron.CoordinateSystem()
coordinate_system.CreateStart(coordinate_system_user_number)
coordinate_system.DimensionSet(dimensions)
coordinate_system.CreateFinish()

2. Creating basis functions

The finite element description of our fields requires a basis function to interpolate field values over elements. In OpenCMISS-Iron, we start by initalising a basis object on which can specific an interpolation scheme. In the following section, you can choose from a number of basis function interpolation types. These are defined in the OpenCMISS-Iron Basis Interpolation Specifications Constants. The interpolation_type variable that was set at the top of the tutorial will be used to determine which basis interpolation will be used in this simulation. Note that you will need to ensure that an appropriate number_of_gauss_per_xi have been set if you change the basis interpolation.

We have also initialised a second pressure_basis object for the hydrostatic pressure. In mechanics theory, it is generally well understood that the interpolation scheme for the hydrostatic pressure basis should be one order lower than the geometric basis. use_pressure_basis, which is set at the top of the tutorial, determines whether the incompressibility constraint will be used in the simulation or not. If the pressure basis is not defined then the simulation is set up for either nearly incompressible constitutive equations or for constitutive equations that describe the mechanical properties of compressible solid materials.

basis_user_number = 1
pressure_basis_user_number = 2

# Define basis parameters.
if dimensions == 1:
    xi_interpolation = [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]
    number_of_guass_xi = [3]
elif dimensions == 2:
    xi_interpolation = [
        iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE,
        iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]
    number_of_guass_xi = [3, 3]
elif dimensions == 3:
    xi_interpolation = [
        iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE,
        iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE,
        iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]
    number_of_guass_xi = [3, 3, 3]

# Define geometric basis.
basis = iron.Basis()
basis.CreateStart(basis_user_number)
if interpolation_type in (1,2,3,4):
    basis.TypeSet(iron.BasisTypes.LAGRANGE_HERMITE_TP)
elif interpolation_type in (7,8,9):
    basis.TypeSet(iron.BasisTypes.SIMPLEX)
basis.NumberOfXiSet(number_of_xi)
basis.InterpolationXiSet(
    [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]*number_of_xi)
if number_of_gauss_per_xi>0:
    basis.QuadratureNumberOfGaussXiSet( [number_of_gauss_per_xi]*number_of_xi)
basis.CreateFinish()

if use_pressure_basis:
    # Define hydrostatic pressure basis.
    pressure_basis = iron.Basis()
    pressure_basis.CreateStart(pressure_basis_user_number)
    if interpolation_type in (1,2,3,4):
        pressure_basis.TypeSet(iron.BasisTypes.LAGRANGE_HERMITE_TP)
    elif interpolation_type in (7,8,9):
        pressure_basis.TypeSet(iron.BasisTypes.SIMPLEX)
    pressure_basis.NumberOfXiSet(number_of_xi)
    pressure_basis.InterpolationXiSet(
        [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]*number_of_xi)
    if number_of_gauss_per_xi > 0:
        pressure_basis.QuadratureNumberOfGaussXiSet(
          [number_of_gauss_per_xi]*number_of_xi)
    pressure_basis.CreateFinish()

3. Creating a region

Next we create a region that our fields will be defined on, and tell it to use the 3D coordinate system we created previously.

# Create a region and assign the coordinate system to the region.
region_user_number = 1
region = iron.Region()
region.CreateStart(region_user_number,iron.WorldRegion)
region.LabelSet("Region")
region.CoordinateSystemSet(coordinate_system)
region.CreateFinish()

4. Setting up a simple cuboid mesh

In this example we will use the iron.GeneratedMesh() class of OpenCMISS to automatically create a 3D geometric mesh on which to solve the mechanics problem. We will create a regular mesh of size width (defined along x), height (defined along y) and depth (defined along z) and divide the mesh into number_global_x_elements in the X direction, number_global_y_elements in the Y direction and number_global_z_elements in the Z direction. We will then tell it to use the basis we created previously.

# Start the creation of a generated mesh in the region.
generated_mesh_user_number = 1
generated_mesh = iron.GeneratedMesh()
generated_mesh.CreateStart(generated_mesh_user_number, region)
generated_mesh.TypeSet(iron.GeneratedMeshTypes.REGULAR)
if use_pressure_basis:
    generated_mesh.BasisSet([basis, pressure_basis])
else:
    generated_mesh.BasisSet([basis])
    generated_mesh.ExtentSet([width, length, height])
    generated_mesh.NumberOfElementsSet(
        [number_global_x_elements, number_global_y_elements, number_global_z_elements])

# Finish the creation of a generated mesh in the region.
mesh_user_number = 1
mesh = iron.Mesh()
generated_mesh.CreateFinish(mesh_user_number,mesh)

5. Decomposing the mesh

Once the mesh has been created we can decompose it into a number of domains in order to allow for parallelism. We choose the options to let OpenCMISS-Iron to calculate the best way to break up the mesh. We also set the number of domains to be equal to the number of computational nodes this example is running on. In this example, we will only be using a single domain. Look for our parallelisation example for a demonstration of how to execute simulations using parallel processing techniques.

# Get the number of computational nodes.
number_of_computational_nodes = iron.ComputationalNumberOfNodesGet()

# Create a decomposition for the mesh.
decomposition_user_number = 1
decomposition = iron.Decomposition()
decomposition.CreateStart(decomposition_user_number,mesh)
decomposition.TypeSet(iron.DecompositionTypes.CALCULATED)
decomposition.NumberOfDomainsSet(number_of_computational_nodes)
decomposition.CreateFinish()

6. Creating a geometric field

Now that the mesh has been decomposed we are in a position to create fields. The first field we need to create is the geometric field. Once we have finished creating the field, we can change the field degrees of freedom (DOFs) to give us our geometry. Since the mesh was constructed using the OpenCMISS-Iron GenerateMesh class, we can use it’s GeometricParametersCalculate method to automatically calculate and populate the geometric field parameters of the regular mesh.

# Create a field for the geometry.
geometric_field_user_number = 1

geometric_field = iron.Field()
geometric_field.CreateStart(geometric_field_user_number,region)
geometric_field.MeshDecompositionSet(decomposition)
geometric_field.TypeSet(iron.FieldTypes.GEOMETRIC)
geometric_field.VariableLabelSet(iron.FieldVariableTypes.U,"Geometry")
geometric_field.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,1)
geometric_field.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,1)
geometric_field.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,1)
if interpolation_type == 4:
    # Set arc length scaling for cubic-Hermite elements.
    geometric_field.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN)
geometric_field.CreateFinish()

# Update the geometric field parameters from generated mesh.
generated_mesh.GeometricParametersCalculate(geometric_field)

Visualising the geometry

We now visualise the geometry using pythreejs. Note that this visualisation currently only supports elements with linear Lagrange interpolation. This includes the node numbers for all elements.

import sys
sys.path.insert(1, '../../../tools/')
import threejs_visualiser
vertices, faces = threejs_visualiser.visualise(
    mesh, geometric_field, dimensions, xi_interpolation, variable=iron.FieldVariableTypes.U, node_labels=True)

Hint: Why are the nodes arranged the way they are? Is there a certain pattern that node ID allocations follow?

Try out several meshes of different element sizes and see if you can answer those questions.

7. Creating fields

Dependent field

When solving the mechanics equations set, we require somewhere to store the deformed geometry (our solution). In OpenCMISS-Iron, we store the solutions to equations sets in a dependent field that contains our dependent variables. Note that the dependent field has been pre-defined in the OpenCMISS-Iron library to contain four components when solving ProblemTypes.FINITE_ELASTICITY with a EquationsSetSubTypes.MOONEY_RIVLIN subtype. The first three components store the deformed coordinates and the fourth stores the hydrostatic pressure.

Remember that use_pressure_basis can be set to true or false to switch between compressible and incompressible elasticity.

One can either make the hydrostatic pressure constant within an element by calling the dependent_field.ComponentInterpolationSet method with the iron.FieldInterpolationTypes.ELEMENT_BASED option. Alternatively, we can make the hydrostatic pressure variable across different nodes by calling the dependent_field.ComponentInterpolationSet method with the iron.FieldInterpolationTypes.NODE_BASED option. In this tutorial, we have chosen the, hydrostatic pressure to be nodally interpolated if use_pressure_basis is true, or use element based interpolation if use_pressure_basis is false.

dependent_field_user_number = 2

dependent_field = iron.Field()
dependent_field.CreateStart(dependent_field_user_number, region)
dependent_field.MeshDecompositionSet(decomposition)
dependent_field.TypeSet(iron.FieldTypes.GEOMETRIC_GENERAL)
dependent_field.GeometricFieldSet(geometric_field)
dependent_field.DependentTypeSet(iron.FieldDependentTypes.DEPENDENT)
dependent_field.VariableLabelSet(iron.FieldVariableTypes.U, "Dependent")
dependent_field.NumberOfVariablesSet(2)

# Set the number of components for the U variable (position) and the DELUDELN
# (forces).
dependent_field.NumberOfComponentsSet(iron.FieldVariableTypes.U, 4)
dependent_field.NumberOfComponentsSet(iron.FieldVariableTypes.DELUDELN, 4)

if use_pressure_basis:
    # Set the hydrostatic pressure to be nodally based and use the second mesh component.
    # U variable (position)
    dependent_field.ComponentInterpolationSet(
        iron.FieldVariableTypes.U, 4, iron.FieldInterpolationTypes.NODE_BASED)
    dependent_field.ComponentMeshComponentSet(
        iron.FieldVariableTypes.U, 4, 2)

    # DELUDELN variable (forces)
    dependent_field.ComponentInterpolationSet(
        iron.FieldVariableTypes.DELUDELN, 4,
        iron.FieldInterpolationTypes.NODE_BASED)
    dependent_field.ComponentMeshComponentSet(
        iron.FieldVariableTypes.DELUDELN, 4, 2)

    if interpolation_type == 4:
        # Set arc length scaling for cubic-Hermite elements.
        dependent_field.FieldScalingTypeSet(
            iron.FieldScalingTypes.ARITHMETIC_MEAN)
else:
    # Set the hydrostatic pressure to be constant within each element.
    dependent_field.ComponentInterpolationSet(
        iron.FieldVariableTypes.U, 4,
        iron.FieldInterpolationTypes.ELEMENT_BASED)
    dependent_field.ComponentInterpolationSet(
        iron.FieldVariableTypes.DELUDELN, 4,
        iron.FieldInterpolationTypes.ELEMENT_BASED)
dependent_field.CreateFinish()

This dependent field needs to be initialised before the simulation is run. To this end, we copy the values of the coordinates from the geometric field into the dependent field in the below code snippet. The hydrostatic pressure field is set to 0.0.

# Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure.
iron.Field.ParametersToFieldParametersComponentCopy(
    geometric_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1,
    dependent_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1)
iron.Field.ParametersToFieldParametersComponentCopy(
    geometric_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 2,
    dependent_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 2)
iron.Field.ParametersToFieldParametersComponentCopy(
    geometric_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 3,
    dependent_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 3)
iron.Field.ComponentValuesInitialiseDP(
    dependent_field, iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 4, 0.0)

Material field

We now set up a new field called the material field, which will store the constitutive equation parameters of the Mooney Rivlin equation. This field can be set to have the same values throughout the mesh to represent a homogeneous material as shown below. If you want to describe heterogeneous materials, you can set the values of these parameters differently across the mesh (e.g. either using an nodally varying variation, or a element constant variation). This is not shown here but we’ll give you a little example of this at the end of this tutorial. Below, the ComponentValuesInitialiseDP function sets all the nodal values to be the same: 2.0 for c10 and 0.0 for c01 for variable type u and the density value for variable type v.

Note that since we have two variables of interest for this field: displacement and density, we must define the two variables in the field using the method: VariableTypesSet(). See iron.py for further documentation.

# Density parameters
#density=0.0009 #in g mm^-3
density = 0.00102 #9345

# Create the material field.
material_field_user_number = 3
material_field = iron.Field()
material_field.CreateStart(material_field_user_number, region)
material_field.TypeSet(iron.FieldTypes.MATERIAL)
material_field.MeshDecompositionSet(decomposition)
material_field.GeometricFieldSet(geometric_field)
material_field.NumberOfVariablesSet(2)
material_field.VariableTypesSet([iron.FieldVariableTypes.U, iron.FieldVariableTypes.V])
material_field.VariableLabelSet(iron.FieldVariableTypes.U, "Material")
material_field.VariableLabelSet(iron.FieldVariableTypes.V, "Density")

# Set the number of components for the Mooney Rivlin constitutive equation (2).
material_field.NumberOfComponentsSet(iron.FieldVariableTypes.U, 2)

# Set a new 1-component variable for density.
material_field.NumberOfComponentsSet(iron.FieldVariableTypes.V, 1)

for component in [1, 2]:
    material_field.ComponentInterpolationSet(
      iron.FieldVariableTypes.U, component,
      iron.FieldInterpolationTypes.ELEMENT_BASED)
    if interpolation_type == 4:
        # Set arc length scaling for cubic-Hermite elements.
        material_field.FieldScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN)
material_field.CreateFinish()

# Set Mooney-Rivlin constants c10 and c01 respectively.
material_field.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, 2.5)
material_field.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 2, 0.0)
material_field.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.V,iron.FieldParameterSetTypes.VALUES, 1, density)

Equation set field

# Equation set field.
equations_set_field_user_number = 4
equations_set_field = iron.Field()
equations_set = iron.EquationsSet()

8. Defining the finite elasticity equations set

We now specify that we want to solve a finite elasticity equation, and identify the specific material constitutive equation that we wish to use to describe the mechanical behaviour of the cube.

The key constants used to define the equation set are: - ProblemClasses.ELASTICITY defines that the equation set is of the elasticity class. - ProblemTypes.FINITE_ELASTICITY defines that the finite elasticity equations set will be used. - EquationsSetSubTypes.MOONEY_RIVLIN defines that the Mooney Rivlin constitutive equation that should be used from the range of constitutive equations that are implemented within the OpenCMISS-Iron library. You can find more information on these by browsing the OpenCMISS-Iron Equations Set Subtypes Constants. Future tutorials will demonstrate how you can dynamically specify constitutive relations using CellML.

Setting up equations set

equations_set_user_number = 1

# Finite elasticity equation specification.
equations_set_specification = [iron.ProblemClasses.ELASTICITY,
    iron.ProblemTypes.FINITE_ELASTICITY,
    iron.EquationsSetSubtypes.MOONEY_RIVLIN]

# Add the geometric field and equations set field that we created earlier (note
# that while we defined the geometric field above, we only initialised an empty
# field for the equations_set_field. When an empty field is provided to the
# equations_set, it will automatically populate it with default values).
equations_set.CreateStart(
    equations_set_user_number, region, geometric_field, equations_set_specification,
    equations_set_field_user_number, equations_set_field)

# Add the dependent field that we created earlier.
equations_set.DependentCreateStart(dependent_field_user_number, dependent_field)
equations_set.DependentCreateFinish()

# Add the material field that we created earlier.
equations_set.MaterialsCreateStart(material_field_user_number, material_field)
equations_set.MaterialsCreateFinish()

Making a source field

We wish to apply the gravity load as a body force, which acts on the entire cube volume. To add this term. another field is created for this problem. In this example, it is called a ‘source’ field. All relevant parameters related to the gravity load are also defined here.

This field uses the equations_set, which hasn’t been completed yet. Hence, remember to finish this process.

Note when defining the gravity vector, the angle is with respect to the positive z axis.

# Set up source field for gravity.
sourceFieldUserNumber = 5

# Gravity parameters.
angle = np.deg2rad(60)
gravity= [0.0, -9.81*np.sin(angle), -9.81*np.cos(angle)] #in ms^-2.
print("Gravity:", gravity)

#Create the source field with the gravity vector
sourceField = iron.Field()
equations_set.SourceCreateStart(sourceFieldUserNumber,sourceField)
if interpolation_type == 4:
    sourceField.fieldScalingType = iron.FieldScalingTypes.ARITHMETIC_MEAN
else:
    sourceField.fieldScalingType = iron.FieldScalingTypes.UNIT
equations_set.SourceCreateFinish()

#Set the gravity vector component values
sourceField.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,gravity[0])
sourceField.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2,gravity[1])
sourceField.ComponentValuesInitialiseDP(
    iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3,gravity[2])

# Finish creating equations set.
equations_set.CreateFinish() # once this is called, cant add anything else to equations_set.

Gravity: [0.0, -8.495709211125343, -4.905000000000001]

Once the equations set is defined, we create the equations that use our fields to construct equations matrices and vectors.

# Create equations.
equations = iron.Equations()
equations_set.EquationsCreateStart(equations)
equations.SparsityTypeSet(iron.EquationsSparsityTypes.SPARSE)
equations.OutputTypeSet(iron.EquationsOutputTypes.NONE)
equations_set.EquationsCreateFinish()

9. Defining the problem

Now that we have defined the equations we can now create our problem to be solved by OpenCMISS-Iron. We create a standard finite elasticity problem, which is a member of the elasticity field problem class. The problem control loop uses the default load increment loop and hence does not require a subtype.

# Define the problem.
problem_user_number = 1

problem = iron.Problem()
problem_specification = (
  [iron.ProblemClasses.ELASTICITY,
   iron.ProblemTypes.FINITE_ELASTICITY,
   iron.ProblemSubtypes.NONE])
problem.CreateStart(problem_user_number, problem_specification)
problem.CreateFinish()

10. Defining control loops

The problem type defines a control loop structure that is used when solving the problem. The OpenCMISS-Iron control loop is a “supervisor” for the computational process. We may have multiple control loops with nested sub loops, and control loops can have different types, for example load incremented loops or time loops for dynamic problems. These control loops have been defined in the OpenCMISS-Iron library for the finite elasticity type of equations as a load increment loop. If we wanted to access the control loop and modify it we would use the problem.ControlLoopGet method before finishing the creation of the control loops. In the below code snippet we get the control loop to set the number of load increments to be used to solve the problem using the variable number_of_load_increments.

# Create the problem control loop.
problem.ControlLoopCreateStart()
control_loop = iron.ControlLoop()
problem.ControlLoopGet([iron.ControlLoopIdentifiers.NODE], control_loop)
control_loop.MaximumIterationsSet(number_of_load_increments)
problem.ControlLoopCreateFinish()

11. Defining solvers

After defining the problem structure we can create the solvers that will be run to actually solve our problem. As the finite elasticity equations are nonlinear, we require a nonlinear solver. Nonlinear solvers typically involve a linearisation step, and therefore a linear solver is also required. In OpenCMISS-Iron, we can start the creation of solvers by calling the problem.SolversCreateStart() method, whose properties can be specified, before they are finalised using a call to the problem.SolversCreateFinish() method. Once finalised, only solver parameter (.e.g tolerances) can be changed, however, fundamental properties (e.g. which solver library to use) cannot be changed. If an additional solver is required, the existing solver can be destroyed and recreated, or another solver can be constructed.

# Number of iterations the solver runs for.
iteration_num = 1000

# Solving process.
nonlinear_solver = iron.Solver()
linear_solver = iron.Solver()
problem.SolversCreateStart()
problem.SolverGet([iron.ControlLoopIdentifiers.NODE], 1, nonlinear_solver)
nonlinear_solver.OutputTypeSet(iron.SolverOutputTypes.NONE)
nonlinear_solver.NewtonJacobianCalculationTypeSet(
  iron.JacobianCalculationTypes.EQUATIONS)
nonlinear_solver.NewtonLinearSolverGet(linear_solver)
nonlinear_solver.NewtonMaximumIterationsSet(iteration_num)
linear_solver.LinearTypeSet(iron.LinearSolverTypes.DIRECT)
problem.SolversCreateFinish()

12. Defining solver equations

After defining our solver we can create the equations for the solver to solve. This is achieved by adding our equations sets to an OpenCMISS-Iron solver_equations object. In this example, we have just one equations set to add, however, for coupled problems, we may have multiple equations sets in the solver equations.

solver = iron.Solver()
solver_equations = iron.SolverEquations()
problem.SolverEquationsCreateStart()
problem.SolverGet([iron.ControlLoopIdentifiers.NODE], 1, solver)
solver.SolverEquationsGet(solver_equations)
solver_equations.SparsityTypeSet(iron.SolverEquationsSparsityTypes.SPARSE)
_ = solver_equations.EquationsSetAdd(equations_set)
#solver_equations.Finalise()
problem.SolverEquationsCreateFinish()
#problem.Finalise()

13. Defining the boundary conditions

The final step in configuring the problem is to define the boundary conditions for the problem. Each line of code below sets a dirichlet boundary conditions that prescribe a nodal coordinate value. The constant iron.BoundaryConditionsTypes.FIXED indicates that the value needs to be fixed to a certain value specified in the final argument of the boundary_conditions.AddNode method.

Note that for the cantilever beam example, the bottom face is fixed in the x, y, and z directions. Thus. a boundary condition value of 0 will be added to whichever node has a z-coordinate value of 0 (see the visualisation above).

To understand how to automatically prescribe these, please run the following blocks of code associated with the nodes and elements in a mesh.

Number of nodes in a mesh.

nodes = iron.MeshNodes()
mesh.NodesGet(1, nodes)
num_nodes = nodes.NumberOfNodesGet()
node_nums = (np.arange(num_nodes)+1).tolist()

A recommended approach would be to print the outputs of the commands and review the associated documentation in iron.py to understand what the code is doing, as it relates to how the boundary conditions are being prescribed.

Number of elements in a mesh.

elements = iron.MeshElements()
mesh.ElementsGet(1, elements)
num_elements = mesh.NumberOfElementsGet()
element_nums = (np.arange(num_elements)+1).tolist()

Prescribe boundary conditions.

boundary_conditions = iron.BoundaryConditions()
solver_equations.BoundaryConditionsCreateStart(boundary_conditions)

# Applying zero-displacement boundary condition on the nodes of the bottom face.
for node_id, node in enumerate(range(1,num_nodes+1)):
    value = dependent_field.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, VersionNumber, DerivativeNumber, int(node), 3)
    if np.allclose(value, 0.0, atol = 1e-06):
        for component_idx, component in enumerate(range(X,Z+1)):
            boundary_conditions.AddNode(dependent_field, iron.FieldVariableTypes.U, VersionNumber, DerivativeNumber, int(node), component, iron.BoundaryConditionsTypes.FIXED, 0.0)
            fixed_nodes = int(node)

# Number of fixed nodes.
print(fixed_nodes)

36

Record reference configuration information.

For evaluating nodal displacement, we wish to record the node’s coordinates before and after deformation. The following code records the nodal coordinates in the undeformed state for each respective node.

# Initialise nodes and reference coordinate lists.
node_value = []
reference = []

# Reference configuration coordinates for all nodes.
for node_id, node in enumerate(range(1,num_nodes+1)):
    for component_idx, component in enumerate(range(X,Z+1)):
        coordinate = dependent_field.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, VersionNumber, DerivativeNumber, int(node), component)
        if component == 1:
            x_value = coordinate
        elif component == 2:
            y_value = coordinate
        elif component == 3:
            z_value = coordinate
        node_value.append(node)
        reference.append(coordinate)

    #print(x_value, y_value, z_value)

    # Let's pick a corner node at the free end, with nodal coordinates (0,0,height).
    if (x_value, y_value, z_value) == (0, 0, float(height)):
        node_want = int(node)

# Convert list to array.
node_value = np.array(node_value)
reference = np.array(reference)

print(node_want)

181

We then construct the solver matrices and vectors by making a call to the solver_equations.BoundaryConditionsCreateFinish() method.

solver_equations.BoundaryConditionsCreateFinish()

14. Solving the problem

After our problem solver equations have been fully defined, we are now ready to solve our problem. When we call the Solve method of the problem it will loop over the control loops and control loop solvers to solve our problem.

# Solve the problem.
try:
    problem.Solve()
    print("Problem solved.")

except:
    print("Too few iterations.")

Problem solved.

Visualising results

We can now visualise the resulting deformation as animation using pythreejs.

vertices, faces = threejs_visualiser.visualise(
    mesh, geometric_field, dimensions, xi_interpolation, dependent_field=dependent_field,
    variable=iron.FieldVariableTypes.U, resolution=8, mechanics_animation=True)

If the simulation doesn’t show any obvious deformation, consider the following: - Check the nodal displacement outputs for the corner nodes. - Are the boundary conditions prescribed correctly? - How is the gravity load implemented? - Is the density value large enough? Is it within the range for the problem you are solving? - Decrease the stiffness value?

Writing vertices and cells of triangular mesh components to an obj file for external visualiser.

cells = [("triangle", np.array([[0, 1, 2]]))]

# Creating triangular cells on each face of geometry.
cells = []
for face in faces:
    cells.append(("triangle", np.array([face])))

# Writing vertices and cell data points to an external file.
import meshio
meshio.write_points_cells(
    "foo.obj",
    vertices,
    cells)

Once an obj file has been created, view the output using online 3D viewers or softwares such as ParaView.

Extracting nodal displacements (corner nodes) for the model.

The following code records the nodal coordinates in the deformed state for each respective node.

# Initialise nodes and deformed coordinates lists.
deformed = []

for node_id, node in enumerate(range(1,num_nodes+1)):
    for component_idx, component in enumerate(range(X,Z+1)):
        coordinate = dependent_field.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, VersionNumber, DerivativeNumber, int(node), component)
        deformed.append(coordinate)

# Convert list to array.
deformed = np.array(deformed)

# Nodal displacement info.
displacement = deformed - reference
node_displacement = displacement[3*(node_want - 1):3*node_want]

# Outputting the relevant info.
print("Number of elements:", num_elements)
print("Node:", node_want)
print('Displacements (mm):', abs(node_displacement))

Number of elements: 125
Node: 181
Displacements (mm): [ 0.25217352 37.8664746  19.51400273]

QUESTION: From inspection, what do these values mean physically? Is the model behaving as you expect?

Degrees of Freedom in the mesh

Degrees of Freedom (DOF) refers to the number of degrees the nodes in a mesh can move freely. These are constrained by the nodal displacements in different directions and the pressure field in each element.

For linear lagragian interpolation, each node in the mesh can move in three different directions (DOF = 3) whilst the fixed nodes on the boundary have no degrees of freedom (DOF = 0). There is one pressure field for each element, check out the computation below.

DoFs = 3*(num_nodes - fixed_nodes) + num_elements
print("Degrees of Freedom:", DoFs)

Degrees of Freedom: 665

15. Exporting results

Before we export the results in Cmgui format, we will first create a new deformed_field and pressure_field to separately hold the solution for the deformed geometry and hydrostatic pressure for visualisation in Cmgui (this enables simplified access to these fields in Cmgui visualisation scripts).

deformed_field_user_number = 6
deformed_field = iron.Field()
deformed_field.CreateStart(deformed_field_user_number, region)
deformed_field.MeshDecompositionSet(decomposition)
deformed_field.TypeSet(iron.FieldTypes.GEOMETRIC)
deformed_field.VariableLabelSet(iron.FieldVariableTypes.U, "DeformedGeometry")
for component in [1, 2, 3]:
    deformed_field.ComponentMeshComponentSet(
            iron.FieldVariableTypes.U, component, 1)
if interpolation_type == 4:
    # Set arc length scaling for cubic-Hermite elements.
    deformed_field.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN)
deformed_field.CreateFinish()

pressure_field_user_number = 7
pressure_field = iron.Field()
pressure_field.CreateStart(pressure_field_user_number, region)
pressure_field.MeshDecompositionSet(decomposition)
pressure_field.VariableLabelSet(iron.FieldVariableTypes.U, "Pressure")
pressure_field.ComponentMeshComponentSet(iron.FieldVariableTypes.U, 1, 1)
pressure_field.ComponentInterpolationSet(
  iron.FieldVariableTypes.U, 1, iron.FieldInterpolationTypes.ELEMENT_BASED)
pressure_field.NumberOfComponentsSet(iron.FieldVariableTypes.U, 1)
pressure_field.CreateFinish()

# Copy deformed geometry into deformed field.
for component in [1, 2, 3]:
    dependent_field.ParametersToFieldParametersComponentCopy(
        iron.FieldVariableTypes.U,
        iron.FieldParameterSetTypes.VALUES, component,
        deformed_field, iron.FieldVariableTypes.U,
        iron.FieldParameterSetTypes.VALUES, component)

# Copy the hydrostatic pressure solutions from the dependent field into the
# pressure field.
dependent_field.ParametersToFieldParametersComponentCopy(
  iron.FieldVariableTypes.U,
  iron.FieldParameterSetTypes.VALUES, 4,
  pressure_field, iron.FieldVariableTypes.U,
  iron.FieldParameterSetTypes.VALUES, 1)

# Export results to exnode and exelem format.
fields = iron.Fields()
fields.CreateRegion(region)
fields.NodesExport("cube", "FORTRAN")
fields.ElementsExport("cube", "FORTRAN")
fields.Finalise()

Finalising session

problem.Destroy()
coordinate_system.Destroy()
region.Destroy()
basis.Destroy()
iron.Finalise()

Additional exercises

Run the Jupyter notebook with a larger number of elements, how does the nodal displacement change? At what number of elements will nodal displacements converge?

Contributors

Tutorial created by Max Dang Vu, with input and guidance from Dr. Prasad Babarenda Gamage and the OpenCMISS team.