Post-Processing¶
Strain¶
Remember that the strain \(\boldsymbol{\varepsilon}\) is defined as the symmetric gradient of the displacements \(\boldsymbol{u}\):
or in index notation
In our finite element implementation we used the Voigt notation for the strains and stresses, hence this formula has to be adapted according to this notation. This is will be “automatically” taken care of when we use the shape function gradient matrix \(\boldsymbol{B}\). When you are using the finite element method you typically compute the strains on the element level using the symmetric shape function gradients:
Implement a method in the Element
class that computes the
strain vector at the position \(xi\), \(eta\) according to the
following specifications:
-
Element.
strain
(xi, eta, u)¶ Return the strain vector of the element subject to displacements u at the reference space position (xi, eta).
Parameters: - xi (float) – Reference space coordinate at which the strain should be computed.
- eta (float) – Reference space coordinate at which the strain should be computed.
- u (numpy.ndarray) – Global displacement vector that results from
solve()
.
Returns: The strain vector of the element subject to displacements u at position (xi, eta) given as \(\begin{bmatrix} \varepsilon_{xx} & \varepsilon_{yy} & \varepsilon_{xy} \end{bmatrix}^{\mathrm{T}}\).
Return type:
Stress¶
Using the strains computing the stresses is straightforward via the popular Hooke’s law
which in terms of a finite element is usually handled via
Implement a method in the Element
class that computes the
stress vector at the position \(xi\), \(eta\) according to the
following specifications:
-
Element.
stress
(xi, eta, u)¶ Return the stress vector of the element subject to displacements u at the reference space position (xi, eta).
Parameters: - xi (float) – Reference space coordinate at which the stress should be computed.
- eta (float) – Reference space coordinate at which the stress should be computed.
- u (numpy.ndarray) – Global displacement vector that results from
solve()
.
Returns: The stress vector of the element subject to displacements u at position (xi, eta) given as \(\begin{bmatrix} \sigma_{xx} & \sigma_{yy} & \sigma_{xy} \end{bmatrix}^{\mathrm{T}}\).
Return type:
Plotting Fields¶
Now that you have your finite element code up and running it is time to see the results of your hard work. For this a function is provided:
-
plot_mesh
(ax, elements, u=None, u_amp=1, field=None, field_label=None, field_resolution=4, padding=0.1, show_axes=True, verbose=True)¶ Plot a finite element mesh and – if given – its solution and fields.
Parameters: - ax (matplotlib.axes.Axes) – The axes which the mesh should be plotted in.
- elements (list[fem.element.Element]) – List of elements that make up the mesh that should be plotted.
- u (numpy.ndarray) – Solution of the FEM system returned by
solve()
. - u_amp (float) – Amplification of the displacements u.
- field (function) – Function with the interface
f(element, xi, eta, u)
where element is anElement
, xi and eta are reference coordinates and u is the solution of the FEM system returned bysolve()
. - field_label (str) – Name of the field that should be plotted.
- field_resolution (int) – Amount of partitions the field interpolation should use.
- padding (float) – The padding of the axis figure.
- show_axes (bool) – Whether to show the axes or not.
- verbose (bool) – Whether to show the element and node numbers or not.
Using this you can plot the results:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from fem.node import Node
from fem.element import Element
from fem.system import System
from fem.utilities import plot_mesh, refine_elements
# Make sure no nodes are stored internally.
Node.clear()
# Initialize an element defining the geometry
initial_elements = [Element([Node(0, 1), Node(0, 0), Node(2, 0.5), Node(2, 1)],
3e7, 0.3)]
# Refine the initial elements
elements = refine_elements(initial_elements, 3)
system = System(elements)
# Define boundary conditions
for node in system.nodes:
x, y = node.coordinates
if np.isclose(x, 0):
system.dirichlet_bc[2*node.id] = 0
system.dirichlet_bc[2*node.id+1] = 0
if np.isclose(y, 1):
if np.isclose(x, 0) or np.isclose(x, 2):
system.neumann_bc[2*node.id+1] = -20
else:
system.neumann_bc[2*node.id+1] = -40
u = system.solve()
# Define the function of the field that should be plotted.
def von_mises_stress(element, xi, eta, u):
s = element.stress(xi, eta, u)
return np.sqrt(s[0]**2 - s[0]*s[1] + s[1]**2 + 3*s[2]**2)
# Plot the mesh
fig, ax = plt.subplots()
plot_mesh(ax, elements, u=u, u_amp=500, field=von_mises_stress,
field_label='von Mises stress', field_resolution=11, show_axes=True,
verbose=False)
plt.show()
