Documentation for Release 2024.2

# Solving More Elasticity Problems

This tutorial is automatically generated from TestSolvingMoreElasticityProblemsTutorial.hpp at revision 3c544f98da9c. Note that the code is given in full at the bottom of the page.

### Introduction

In this second solid mechanics tutorial, we illustrate some other possibilities: using tractions that are defined with a function, or tractions that depend on the deformed body (eg normal pressure boundary conditions), specifying non-zero displacement boundary conditions, and then displacement boundary conditions only in some directions, and doing compressible solves.

These includes are the same as before

The incompressible solver

An incompressible material law

These two are specific to compressible problems

This include should generally go last to avoid issues on old library versions

This function is used in the first test

### Incompressible deformation: non-zero displacement boundary conditions, functional tractions

We now consider a more complicated example. We prescribe particular new locations for the nodes on the Dirichlet boundary, and also show how to prescribe a traction that is given in functional form rather than prescribed for each boundary element.

Create a mesh

Use a different material law this time, an exponential material law.
The material law needs to inherit from `AbstractIncompressibleMaterialLaw`

,
and there are a few implemented, see `continuum_mechanics/src/problem/material_laws`

Now specify the fixed nodes, and their new locations. Create `std::vector`

s for each.

Loop over the mesh nodes

If the node is on the Y=0 surface (the LHS)

Add it to the list of fixed nodes

and define a new position x=(X,0.1*X^2^)

Now collect all the boundary elements on the top surface, as before, except here we don’t create the tractions for each element

If Y=1, have found a boundary element

Create a problem definition object, and this time calling `SetFixedNodes`

which takes in the new locations of the fixed nodes.

Now call `SetTractionBoundaryConditions`

, which takes in a vector of
boundary elements as in the previous test. However this time the second argument
is a *function pointer* (just the name of the function) to a
function returning traction in terms of position (and time [see below]).
This function is defined above, before the tests. It has to take in a `c_vector`

(position)
and a double (time), and returns a `c_vector`

(traction), and will only be called
using points in the boundary elements being passed in. The function `MyTraction`

(defined above, before the tests) above defines a horizontal traction (ie a shear stress, since it is
applied to the top surface) which increases in magnitude across the object.

Note: You can also call `problem_defn.SetBodyForce(MyBodyForce)`

, passing in a function
instead of a vector, although isn’t really physically useful, it is only really useful
for constructing problems with exact solutions.

Create the solver as before

Call `Solve()`

Another quick check

Visualise as before.

**Advanced:** Note that the function `MyTraction`

takes in time, which it didn’t use. In the above it would have been called
with t=0. The current time can be set using `SetCurrentTime()`

. The idea is that the user may want to solve a
sequence of static problems with time-dependent tractions (say), for which they should allow `MyTraction`

to
depend on time, and put the solve inside a time-loop, for example:

In this the current time would be passed through to `MyTraction`

Create Cmgui output

This is just to check that nothing has been accidentally changed in this test

### Sliding boundary conditions

It is common to require a Dirichlet boundary condition where the displacement/position in one dimension is fixed, but the displacement/position in the others are free. This can be easily done when collecting the new locations for the fixed nodes, as shown in the following example. Here, we take a unit square, apply gravity downward, and suppose the Y=0 surface is like a frictionless boundary, so that, for the nodes on Y=0, we specify y=0 but leave x free (Here (X,Y)=old position, (x,y)=new position). (Note though that this wouldn’t be enough to uniquely specify the final solution - an arbitrary translation in the Y direction could be added a solution to obtain another valid solution, so we fully fix the node at the origin.)

Create fixed nodes and locations…

Fix node 0 (the node at the origin)

For the rest, if the node is on the Y=0 surface..

..add it to the list of fixed nodes..

..and define y to be 0 but x is fixed

Set the material law and fixed nodes, add some gravity, and solve

Check the node at (1,0) has moved but has stayed on Y=0

### Compressible deformation, and other bits and pieces

In this test, we will show the (very minor) changes required to solve a compressible nonlinear
elasticity problem, we will describe and show how to specify ‘pressure on deformed body’
boundary conditions, we illustrate how a quadratic mesh can be generated using a linear mesh
input files, and we also illustrate how `Solve()`

can be called repeatedly, with loading
changing between the solves.

Note: for other examples of compressible solves, including problems with an exact solution, see the
file `continuum_mechanics/test/TestCompressibleNonlinearElasticitySolver.hpp`

All mechanics problems must take in quadratic meshes, but the mesh files for (standard) linear meshes in Triangles/Tetgen can be automatically converted to quadratic meshes, by simply doing the following. (The mesh loaded here is a disk centred at the origin with radius 1).

Compressible problems require a compressible material law, ie one that
inherits from `AbstractCompressibleMaterialLaw`

. The `CompressibleMooneyRivlinMaterialLaw`

is one such example; instantiate one of these

For this problem, we fix the nodes on the surface for which Y < -0.9

We will (later) apply Neumann boundary conditions to surface elements which lie below Y=0, and these are collected here. (Minor, subtle, comment: we don’t bother here checking Y>-0.9, so the surface elements collected here include the ones on the Dirichlet boundary. This doesn’t matter as the Dirichlet boundary conditions to the nonlinear system essentially overwrite an Neumann-related effects).

Create the problem definition class, and set the law again, this time stating that the law is compressible

Set the fixed nodes and gravity

The elasticity solvers have two nonlinear solvers implemented, one hand-coded and one which uses PETSc’s SNES solver. The latter is not the default but can be more robust (and will probably be the default in later versions). This is how it can be used. (This option can also be called if the compiled binary is run from the command line (see Running Binaries From Command Line) using the option “-mech_use_snes”).

This line tells the solver to output info about the nonlinear solve as it progresses, and can be used with or without the SNES option above. The corresponding command line option is “-mech_verbose”

Declare the compressible solver, which has the same interface as the incompressible
one, and call `Solve()`

Now we call add additional boundary conditions, and call `Solve() again. Firstly: these
Neumann conditions here are not specified traction boundary conditions (such BCs are specified
on the undeformed body), but instead, the (more natural) specification of a pressure
exactly in the *normal direction on the deformed body*. We have to provide a set of boundary
elements of the mesh, and a pressure to act on those elements. The solver will automatically
compute the deformed normal directions on which the pressure acts. Note: with this type of
BC, the ordering of the nodes on the boundary elements needs to be consistent, otherwise some
normals will be computed to be inward and others outward. The solver will check this on the
original mesh and throw an exception if this is not the case. (The required ordering is such that:
in 2D, surface element nodes are ordered anticlockwise (looking at the whole mesh); in 3D, looking
at any surface element from outside the mesh, the three nodes are ordered anticlockwise. (Triangle/tetgen
automatically create boundary elements like this)).

Call `Solve()`

again, so now solving with fixed nodes, gravity, and pressure. The solution from
the previous solve will be used as the initial guess. Although at the moment the solution from the
previous call to `Solve()`

will be over-written, calling `Solve()`

repeatedly may be useful for
some problems: sometimes, Newton’s method will fail to converge for given force/pressures etc, and it can
be (very) helpful to *increment* the loading. For example, set the gravity to be (0,-9.81/3), solve,
then set it to be (0,-2*9.81/3), solve again, and finally set it to be (0,-9.81) and solve for a
final time