auxiliary_logarithmic_stress_patch.py

Go to the documentation of this file.

"""Affine patch with three unequal principal stretches and a finite rotation."""
 
import numpy as np
 
 
LOG_STRETCH = np.array([0.04, -0.015, 0.005])
ROTATION = np.array([0.02, -0.015, 0.01])
C10 = 1.7
K = 8.5
 
 
def deformation_gradient():
    angle = np.linalg.norm(ROTATION)
    x, y, z = ROTATION
    # MoFEM uses hat(omega)_ij = epsilon_ijk omega_k.
    skew = np.array([[0, z, -y], [-z, 0, x], [y, -x, 0]])
    rotation = (np.eye(3) + np.sin(angle) / angle * skew
                + (1 - np.cos(angle)) / angle**2 * skew @ skew)
    return rotation @ np.diag(np.exp(LOG_STRETCH))
 
 
def analytical_disp(delta_t, t, x, y, z, tx, ty, tz, block_name):
    points = np.column_stack((x, y, z))
    return np.ascontiguousarray(t * points @ (deformation_gradient() - np.eye(3)).T)
 
 
def analytical_traction(delta_t, t, x, y, z, tx, ty, tz, block_name):
    gradient = np.eye(3) + t * (deformation_gradient() - np.eye(3))
    jacobian = np.linalg.det(gradient)
    inverse_transpose = np.linalg.inv(gradient).T
    stress = (2 * C10 * jacobian**(-2 / 3)
              * (gradient - np.sum(gradient**2) / 3 * inverse_transpose)
              + K * jacobian * (jacobian - 1) * inverse_transpose)
    normals = np.column_stack((tx, ty, tz))
    normals /= np.linalg.norm(normals, axis=1)[:, None]
    return np.ascontiguousarray(normals @ stress.T)
 
 
def analytical_external_strain(delta_t, t, x, y, z, block_name):
    return np.zeros((len(x), 1))
 
 
def analytical_elastic(delta_t, t, x, y, z, block_name):
    return np.ones((len(x), 1))