Regularization Mode

This example will explain the basic usage of a regularisation with saenopy.

Regularization mode solves the unconstrained problem of measured displacements, which should be fitted by deforming the mesh and thereby generating forces. To avoid spurious forces all over the mesh, a regularization term is used to suppress “noise” forces.

The input is:
  • the material parameters (see Material)

  • the mesh (nodes and connectivity) (see Mesh)

  • the measured displacements \(U^\mathrm{measured}_i\) for all nodes

The output is:
  • the forces \(f_i(U)\) for all nodes

  • the displacements \(U_i\) for all nodes (should be similar to the measured displacements)

This mode is used to calculate from a measured noisy displacement field, the forces that generated this displacement field.

First the basic theoretical background is explained, followed by the code to use it.

Theory

The inverse problem, to fit forces for measured displacements is under determined. Therefore, a regularization scheme is needed, to circumvent this problem. The target function \(L(U)\) is extended with a regularisation term. We use here the Thikonov regularization, \(|f_i(U)|^2_A\), the weighted sum over all forces (other regularization schemes can also be used).

\[L(U) = | U_i - U_i^\mathrm{measured}|^2 + |f_i(U)|^2_A\]

The subscript \(A\) denotes here a weighted sum over the forces. If \(A\) would be a constant, forces at every point of the volume would be penalized the same. But we normally expect strong forces at a few nodes, caused by the cell, and small forces at other nodes, caused by noise.

To not smooth the strongest forces, a lower weight is assigned to nodes which obtained a high force in the iteration process (Huber, 2004).

\[\begin{split}A_{ii}(f_i) = \begin{cases} \alpha, & \text{if}\ |f_i| < 1.345 \cdot \mathrm{median}(|f|) \\ \alpha \frac{1.345 \cdot \mathrm{median}(|f|)}{|f_i|}, & \text{otherwise} \end{cases}\end{split}\]

The target function \(L(U)\) is minimized if \(U\) fulfills the following condition:

\[\underbrace{(\mathbf{I} + \mathbf{K}(U) \cdot \mathbf{A} \cdot \mathbf{K}(U))_{ij}}_{A_{ij}} \cdot \underbrace{\Delta u_j}_{x_j} = \underbrace{U_i^\mathrm{measured} + (\mathbf{K}(U) \cdot \mathbf{A} \cdot f(U))_i}_{b_i}\]

This linear equation (of the form \(A_{ij}\cdot x_j = b_i\)) is solved using the conjugate gradient method to obtain a value for \(\Delta u_i\). A tiny fraction (stepper) of this is applied to the displacement \(U_i\):

\[U^\prime_i = U_i + \mathrm{stepper} \cdot \Delta u_i.\]

With the new displacement \(U^\prime\), the stiffness matrix \(K_{ij}(U^\prime)\), the nodal forces \(f_i(U^\prime)\) and the total energy \(E(U^\prime)\) are updated.

From these nodal forces the weight matrix \(A_{ii}(f_i(U^\prime))\) is updated. And the linear equation is again solved for the new stiffness matrix \(K_{ij}(U^\prime)\) and weight matrix \(A_{ii}(f_i(U^\prime))\).

This procedure is iterated until the total energy \(E\) of the system converges. The convergence criterion is:

\[\frac{\mathrm{std}(\{E^{t-6}, E^{t-5}, ..., E^{t-1}, E^t\})}{\mathrm{mean}(\{E^{t-6}, E^{t-5}, ..., E^{t-1}, E^t\})} \leq \mathrm{rel\_conv\_crit}\]

So when the standard deviation of \(E\) divided by the mean of \(E\) for the last 6 iterations is lower than the threshold rel_conv_crit.

Example

The following code example (the code of all boxed joined in a single file) can be downloaded at regularization.py.

First, import the Solver class.

[1]:
from saenopy import Solver

# initialize the object
M = Solver()

Set the material model (see Material):

[2]:
from saenopy.materials import SemiAffineFiberMaterial

# provide a material model
material = SemiAffineFiberMaterial(1645, 0.0008, 1.0075, 0.033)
M.setMaterialModel(material)

Define the mesh (see Mesh).

[3]:
import numpy as np

# define the coordinates of the nodes of the mesh
# the array has to have the shape N_v x 3
R = np.array([[0., 0., 0.],  # 0
              [0., 1., 0.],  # 1
              [1., 1., 0.],  # 2
              [1., 0., 0.],  # 3
              [0., 0., 1.],  # 4
              [1., 0., 1.],  # 5
              [1., 1., 1.],  # 6
              [0., 1., 1.]]) # 7

# define the tetrahedra of the mesh
# the array has to have the shape N_t x 4
# every entry is an index referencing a verces in R (indices start with 0)
T = np.array([[0, 1, 7, 2],
              [0, 2, 5, 3],
              [0, 4, 5, 7],
              [2, 5, 6, 7],
              [0, 7, 5, 2]])

And hand the data over to the Solver object.

[4]:
# provide the node data
M.setNodes(R)
# and the tetrahedron data
M.setTetrahedra(T)

Now we have to specify which displacements to fit (see Measured displacement).

[5]:
# the displacements of the nodes which shall be fitted
# during the solving
U = np.array([[0   , 0, 0],  # 0
              [0   , 0, 0],  # 1
              [0.01, 0, 0],  # 2
              [0.01, 0, 0],  # 3
              [0   , 0, 0],  # 4
              [0.01, 0, 0],  # 5
              [0.01, 0, 0],  # 6
              [0   , 0, 0]]) # 7

# hand the displacements over to the class instance
M.setTargetDisplacements(U)

Now we can start the regularisation process (see parameters see solve_regularized()).

[6]:
# call the regularisation
M.solve_regularized(stepper=0.1, alpha=0.001);
e:\saenopy\saenopy\solver.py:339: NumbaWarning: Cannot cache compiled function "_get_applied_epsilon" as it uses dynamic globals (such as ctypes pointers and large global arrays)
  @staticmethod

We can now view the result with the integrated mesh viewer. The nodes (yellow points) are connected by tetrahedra (yellow lines). The displacements (green vectors) cause the reaction forces (red vectors).

[7]:
M.viewMesh(50, 1)