Linear elastic#

Note

For readability the examples are only presented on the Python API. Converting them to C++ is trivial: the API is (almost) identical, and also the API of xtensor is almost identical to that of NumPy.

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.

sketch.svg

Below, an example is described line-by-line. For simplicity the material model is taken from GMatElastic, as this example is about FEM, not about constitutive modelling. The full example can be downloaded run as follows:

example.py

In this example we illustrate some basic plotting.

Proceed for example as follows:

Get the prerequisites:

conda install -c conda-forge python-goosefem python-gmatelastic goosempl

Run (and plot):

python fixed-displacement_elastic.py --plot

Include library#

import argparse
import sys

import GMatElastic
import GooseFEM
import numpy as np

Define mesh#

# define mesh
mesh = GooseFEM.Mesh.Quad4.Regular(5, 5)

# mesh dimensions
nelem = mesh.nelem
nne = mesh.nne
ndim = mesh.ndim

# mesh definition, displacement, external forces
coor = mesh.coor
conn = mesh.conn
dofs = mesh.dofs
disp = np.zeros_like(coor)
fext = np.zeros_like(coor)

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 (nodesRightEdge, nodesTopEdge, nodesLeftEdge, nodesBottomEdge).

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#

iip = np.concatenate(
    (
        dofs[mesh.nodesRightEdge, 0],
        dofs[mesh.nodesTopEdge, 1],
        dofs[mesh.nodesLeftEdge, 0],
        dofs[mesh.nodesBottomEdge, 1],
    )
)

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

\[\begin{split}\texttt{u} = \begin{bmatrix} \texttt{u}_u \\ \texttt{u}_p \end{bmatrix}\end{split}\]

where the subscript \(u\) and \(p\) 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#

vector = GooseFEM.VectorPartitioned(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’ \(\leftrightarrow\) ‘dofval’ using dofs and iip,
‘nodevec’ \(\leftrightarrow\) ‘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

System matrix#

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).

See also

Matrix representation
Details: Matrix

Element definition#

elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor))
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

Material definition#

mat = GMatElastic.Cartesian3d.Elastic2d(K=np.ones([nelem, nip]), G=np.ones([nelem, nip]))

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:

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

Compute strain#

ue = vector.AsElement(disp)
elem.symGradN_vector(ue, mat.Eps)
mat.refresh()

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

The displacement per element ue eliminate the connectivity from the equation for a very significant part of the API. This makes things much easier to control and much more modular.

Warning

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

Note

A function like

ue = vector.AsElement(disp)

starting with a capital letter, allocates the output data. However, if that data was already allocated, it can be modified in-place with the corresponding function starting with a lowercase letter, e.g.

vector.asElement(disp, ue)

Tip

To allocate nodal vectors or tensors per element, use the convenience functions:

# nodal vectors ("fint", "fext", "fext", "disp", or "coor")
variable = np.zeros(vector.shape_nodevec())

# vector per element ("ue" or "fe")
variable = np.zeros(vector.shape_elemvec())

# matrix per element ("Ke")
variable = np.zeros(vector.shape_elemmat())

Assemble system#

fe = elem.Int_gradN_dot_tensor2_dV(mat.Sig)
fint = vector.AssembleNode(fe)

# stiffness matrix
Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(mat.C)
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 \(\rightarrow\) fint and Ke \(\rightarrow\) 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.

Solve#

# residual
fres = fext - fint

# set fixed displacements
disp[mesh.nodesRightEdge, 0] = +0.2
disp[mesh.nodesTopEdge, 1] = -0.2
disp[mesh.nodesLeftEdge, 0] = 0.0  # not strictly needed: default == 0
disp[mesh.nodesBottomEdge, 1] = 0.0  # not strictly needed: default == 0

# solve
Solver.solve(K, fres, disp)

We now compute the residual forces (trivially zero in this case, but kept for pedagogical reasons). We then 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. As an example, the same operation with manual book-keeping is included in

example.py

Post-process#

Strain and stress#

vector.asElement(disp, ue)
elem.symGradN_vector(ue, mat.Eps)
mat.refresh()

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

Residual force#

# internal force
elem.int_gradN_dot_tensor2_dV(mat.Sig, fe)
vector.assembleNode(fe, fint)

# apply reaction force
vector.copy_p(fint, fext)

# residual
fres = fext - fint

# print residual
assert np.isclose(np.sum(np.abs(fres)) / np.sum(np.abs(fext)), 0)

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

Plot#

Let’s extract the average (equivalent) stress per element:

    Sigav = np.average(mat.Sig, weights=dV, axis=1)
    sigeq = GMatElastic.Cartesian3d.Sigeq(Sigav)

And plot an equivalent stress on a deformed mesh:

    gplt.patch(coor=coor, conn=conn, linestyle="--", axis=ax)