# Linear elastic¶

Consider a uniform linear elastic bar that is extended by a uniform fixed displacement on both sides. This problem can be modelled and discretised using symmetry as show below. In this example we will furthermore assume that the bar is sufficiently thick in the out-of-plane direction to be modelled using two-dimensional plane strain.

Below an example is described line-by-line. The full example can be downloaded:

Note

This example is also available using the Python interface
(`example.py`

).
Compared to the C++ API, the Python API requires more data-allocation,
in particular for the functions *AsElement* and *AssembleNode*.
See: Data allocation.

## Include library¶

```
#include <GMatElastic/Cartesian3d.h>
#include <GooseFEM/GooseFEM.h>
#include <GooseFEM/MatrixPartitioned.h>
#include <highfive/H5Easy.hpp>
```

The first step is to include the header-only library. Note that for this example we also make use of a material model (GMatElastic) and a library to write (and read) HDF5 files (HighFive).

## Define mesh¶

```
// define mesh
GooseFEM::Mesh::Quad4::Regular mesh(5, 5);
// mesh dimensions
size_t nelem = mesh.nelem();
size_t nne = mesh.nne();
size_t ndim = mesh.ndim();
// mesh definitions
xt::xtensor<double, 2> coor = mesh.coor();
xt::xtensor<size_t, 2> conn = mesh.conn();
xt::xtensor<size_t, 2> dofs = mesh.dofs();
// node sets
xt::xtensor<size_t, 1> nodesLft = mesh.nodesLeftEdge();
xt::xtensor<size_t, 1> nodesRgt = mesh.nodesRightEdge();
xt::xtensor<size_t, 1> nodesTop = mesh.nodesTopEdge();
xt::xtensor<size_t, 1> nodesBot = mesh.nodesBottomEdge();
```

A mesh is defined using GooseFEM.
As observed the *mesh* is a class that has methods to extract the relevant information
such as the nodal coordinates (*coor*), the connectivity (*conn*),
the degrees-of-freedom per node (*dofs*) and several node-sets that
will be used to impose the sketched boundary conditions
(*nodesLeft*, *nodesRight*, *nodesTop*, *nodesBottom*).

Note that:

- The connectivity (
*conn*) contains information of which nodes, in which order, belong to which element. - The degrees-of-freedom per node (
*dofs*) contains information of how a nodal vector (a vector stored per node) can be transformed to a list of degrees-of-freedom as used in the linear system (although this can be mostly done automatically as we will see below).

See also

- Terminology
- Details: Mesh::Quad4

## Define partitioning¶

```
xt::xtensor<size_t, 1> iip = xt::concatenate(xt::xtuple(
xt::view(dofs, xt::keep(nodesRgt), 0),
xt::view(dofs, xt::keep(nodesTop), 1),
xt::view(dofs, xt::keep(nodesLft), 0),
xt::view(dofs, xt::keep(nodesBot), 1)));
```

We will reorder such that degrees-of-freedom are ordered such that

where the subscript and respectively denote
*Unknown* and *Prescribed*
degrees-of-freedom.
To achieve this we start by collecting all prescribed degrees-of-freedom in *iip*.

## (Avoid) Book-keeping¶

```
GooseFEM::VectorPartitioned vector(conn, dofs, iip);
```

To switch between the three of GooseFEM’s data-representations,
an instance of the *Vector* class is used.
This instance, *vector*, will enable us to switch between a vector field (e.g. the displacement)

- collected per node,
- collected per degree-of-freedom, or
- collected per element.

Note

The *Vector* class collects most, if not all, the burden of book-keeping.
It is thus here that *conn*, *dofs*, and *iip* are used. In particular,

- ‘nodevec’ ‘dofval’ using
*dofs*and*iip*, - ‘nodevec’ ‘elemvec’ using
*conn*.

By contrast, most of GooseFEM’s other methods receive the relevant representation,
and consequently require no problem specific knowledge.
They thus do not have to supplied with *conn*, *dofs*, or *iip*.

See also

- Vector representation
- Data storage
- Details: Vector

## System matrix¶

```
GooseFEM::MatrixPartitioned K(conn, dofs, iip);
GooseFEM::MatrixPartitionedSolver<> Solver;
```

We now also allocate the system/stiffness system (stored as sparse matrix). Like vector, it can accept and return different vector representations, in addition to the ability to assemble from element system matrices.

In addition we allocate the accompanying sparse solver, that we will use to solve a linear system of equations. Note that the solver-class takes care of factorising only when needed (when the matrix has been changed).

Note

Here, the default solver is used (which is the default template, hence the “<>”). To use other solvers see: Linear solver.

See also

- Matrix representation
- Details: Matrix

## Allocate nodal vectors¶

```
xt::xtensor<double, 2> disp = xt::zeros<double>(coor.shape());
xt::xtensor<double, 2> fint = xt::zeros<double>(coor.shape());
xt::xtensor<double, 2> fext = xt::zeros<double>(coor.shape());
xt::xtensor<double, 2> fres = xt::zeros<double>(coor.shape());
```

*disp*: nodal displacements*fint*: nodal internal forces*fext*: nodal external forces*fres*: nodal residual forces

Note

To allocate nodal vectors the convenience function:

```
vector.AllocateNodevec(); // allocate
vector.AllocateElemvec(0.0); // allocate & (zero-)initialise
```

is available, which takes care of getting the right shape. E.g.

```
auto disp = vector.AllocateNodevec(0.0);
```

## Allocate element vectors¶

```
xt::xtensor<double, 3> ue = xt::empty<double>({nelem, nne, ndim});
xt::xtensor<double, 3> fe = xt::empty<double>({nelem, nne, ndim});
xt::xtensor<double, 3> Ke = xt::empty<double>({nelem, nne * ndim, nne * ndim});
```

*ue*: displacement*fe*: force*Ke*: tangent matrix

Warning

Upsizing (e.g. *disp* *ue*) can be done uniquely,
but downsizing (e.g. *fe* *fint*) can be done in more than one way,
see Conversion.
We will get back to this point below.

Note

To allocate element vectors the convenience function:

```
vector.AllocateElemvec(); // allocate
vector.AllocateElemvec(0.0); // allocate & (zero-)initialise
```

is available, which takes care of getting the right shape. E.g.

```
auto ue = vector.AllocateElemvec(0.0);
```

Note

To allocate element matrices the convenience function:

```
vector.AllocateElemmat(); // allocate
vector.AllocateElemmat(0.0); // allocate & (zero-)initialise
```

is available, which takes care of getting the right shape. E.g.

```
auto Ke = vector.AllocateElemmat(0.0);
```

## Element definition¶

```
GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor));
size_t nip = elem.nip();
```

At this moment the interpolation and quadrature is allocated.
The shape functions and integration points (that can be customised) are stored in this class.
As observed, no further information is needed than the number of elements and
the nodal coordinates per element.
Both are contained in the output of `vector.AsElement(coor)`

, which is an ‘elemvec’ of
shape “[nelem, nne, ndim]”.
This illustrates that problem specific book-keeping is isolated to the main program,
using *Vector* as tool.

Note

The shape-functions are computed when constructing this class, they are not recomputed when evaluating them. One can recompute them if the nodal coordinates change using “.update_x(…)”, however, this is only relevant in a large deformation setting.

See also

- Vector representation
- Data storage
- Details: Vector
- Details: Element::Quad4

## Material definition¶

```
GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1.0, 1.0);
```

We now define a uniform linear elastic material, using an external library that is tuned to GooseFEM. This material library will translate a strain tensor per integration point to a stress tensor per integration point and a stiffness tensor per integration point.

See also

Material libraries tuned to GooseFEM include:

- GMatElastic
- GMatElastoPlastic
- GMatElastoPlasticFiniteStrainSimo
- GMatElastoPlasticQPot
- GMatElastoPlasticQPot3d
- GMatNonLinearElastic

But other libraries can also be easily used with (simple) wrappers.

## Integration point tensors¶

```
xt::xtensor<double, 4> Eps = xt::empty<double>({nelem, nip, 3ul, 3ul});
xt::xtensor<double, 4> Sig = xt::empty<double>({nelem, nip, 3ul, 3ul});
xt::xtensor<double, 6> C = xt::empty<double>({nelem, nip, 3ul, 3ul, 3ul, 3ul});
```

These strain, stress, and stiffness tensors per integration point are now allocated. Note that these tensors are 3-d while our problem was 2-d. This is thanks to the plane strain assumption, and the element definition that ignores all third-dimension components.

Note

To allocate integration point the convenience function:

```
quad.AllocateQtensor<rank>(); // allocate
quad.AllocateQtensor<rank>(0.0); // allocate & (zero-)initialise
```

is available, which takes care of getting the right shape. E.g.

```
auto Eps = quad.AllocateQtensor<2>();
auto C = quad.AllocateQtensor<4>();
```

From Python simply use `rank`

as the first function argument.
Furthermore, for scalar you could use `AllocateQscalar()`

which is equivalent to `AllocateQtensor<0>()`

.

## Compute strain¶

```
vector.asElement(disp, ue);
elem.symGradN_vector(ue, Eps);
```

The strain per integration point is now computed using the current nodal displacements
(stored as ‘elemvec’ in *ue*) and the gradient of the shape functions.

Note

*ue* is the output of `vector.asElement(disp, ue)`

.
Using this syntax re-allocation of *ue* is avoided.
If this optimisation is irrelevant for you problem (or if you are using the Python interface),
please use the same function, but starting with a capital:

```
ue = vector.AsElement(disp);
```

Note that this allows the one-liner

```
Eps = elem.SymGradN_vector(vector.AElement(disp));
```

## Compute stress and tangent¶

```
mat.tangent(Eps, Sig, C);
```

The stress and stiffness tensors are now computed for each integration point (completely independently) using the external material model.

Note

*Sig* and *C* are the output variables that were preallocated in the main.

## Assemble system¶

```
// internal force
elem.int_gradN_dot_tensor2_dV(Sig, fe);
vector.assembleNode(fe, fint);
// stiffness matrix
elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke);
K.assemble(Ke);
```

The stress stored per integration point (*Sig*) is now converted to
nodal internal force vectors stored per element (*fe*).
Using *vector* this ‘elemvec’ representation is then converted of a
‘nodevec’ representation in *fint*.
Likewise, the stiffness tensor stored for integration point (*C*) are converted
to system matrices stored per element (‘elemmat’) and finally assembled to
the global stiffness matrix.

Warning

Please note that downsizing
(*fe* *fint* and *Ke* *K*) can be done in two ways,
and that “assemble…” is the right function here as it adds entries that occur
more than once.
In contrast “as…” would not result in what we want here.

Note

Once more, *fe*, *fint*, and *Ke* are output variables. Less efficient, but shorter, is:

```
// internal force
fint = vector.AssembleNode(elem.Int_gradN_dot_tensor2_dV(Sig));
// stiffness matrix
K.assemble(elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C));
```

## Solve¶

```
// set fixed displacements
xt::view(disp, xt::keep(nodesRgt), 0) = +0.1;
xt::view(disp, xt::keep(nodesTop), 1) = -0.1;
xt::view(disp, xt::keep(nodesLft), 0) = 0.0;
xt::view(disp, xt::keep(nodesBot), 1) = 0.0;
// residual
xt::noalias(fres) = fext - fint;
// solve
Solver.solve(K, fres, disp);
```

We now prescribe the displacement of the Prescribed degrees-of-freedom directly
in the nodal displacements *disp* and compute the residual force.
This is follows by partitioning and solving, all done internally in the *MatrixPartitioned* class.

## Post-process¶

### Strain and stress¶

```
vector.asElement(disp, ue);
elem.symGradN_vector(ue, Eps);
mat.stress(Eps, Sig);
```

The strain and stress per integration point are recomputed for post-processing.

### Residual force¶

```
// internal force
elem.int_gradN_dot_tensor2_dV(Sig, fe);
vector.assembleNode(fe, fint);
// apply reaction force
vector.copy_p(fint, fext);
// residual
xt::noalias(fres) = fext - fint;
// print residual
std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl;
```

We convince ourselves that the solution is indeed in mechanical equilibrium.

### Store & plot¶

```
// average stress per node
xt::xtensor<double, 4> dV = elem.AsTensor<2>(elem.dV());
xt::xtensor<double, 3> SigAv = xt::average(Sig, dV, {1});
// write output
H5Easy::File file("output.h5", H5Easy::File::Overwrite);
H5Easy::dump(file, "/coor", coor);
H5Easy::dump(file, "/conn", conn);
H5Easy::dump(file, "/disp", disp);
H5Easy::dump(file, "/Sig", SigAv);
```

Finally we store some fields for plotting using
`plot.py`

.

## Manual partitioning¶

To verify how partitioning and solving is done internally using the *MatrixPartitioned* class,
the same example is provided where partitioning is done manually: