From c5570c8b83fa48bb740cef413116d66a581c5990 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 21 May 2025 23:54:35 +0200 Subject: [PATCH] Enhance the stability checks (#147) * Check stability in all directions * Update __about__.py --- src/matadi/__about__.py | 2 +- src/matadi/_lab_compressible.py | 18 ++++++++-------- src/matadi/_lab_incompressible.py | 35 ++++++++++++++----------------- 3 files changed, 26 insertions(+), 29 deletions(-) diff --git a/src/matadi/__about__.py b/src/matadi/__about__.py index 493f741..260c070 100644 --- a/src/matadi/__about__.py +++ b/src/matadi/__about__.py @@ -1 +1 @@ -__version__ = "0.3.0" +__version__ = "0.3.1" diff --git a/src/matadi/_lab_compressible.py b/src/matadi/_lab_compressible.py index 373a0c9..75c8037 100644 --- a/src/matadi/_lab_compressible.py +++ b/src/matadi/_lab_compressible.py @@ -37,11 +37,11 @@ def stability(stretch, stretches_23): for j, b in enumerate(c): B[i, j] = A[(*a, *b)] - # unit force in direction 1 + # unit forces in all directions # calculate linear solution of stretch 1 resulting from unit load - dl = np.linalg.inv(B)[0, 0] + dl = np.diag(np.linalg.inv(B)) - return dl > 0 + return np.all(dl > 0) def stress_free(stretches_23, stretch): s = stress(stretch, *stretches_23) @@ -77,11 +77,11 @@ def stability(stretch, stretch_3): for j, b in enumerate(c): B[i, j] = A[(*a, *b)] - # unit force in direction 1 + # unit forces in all directions # calculate linear solution of stretch 1 resulting from unit load - dl = np.linalg.inv(B)[0, 0] + dl = np.diag(np.linalg.inv(B)) - return dl > 0 + return np.all(dl > 0) def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] @@ -116,11 +116,11 @@ def stability(stretch, stretch_3): for j, b in enumerate(c): B[i, j] = A[(*a, *b)] - # init unit force in direction 1 + # init unit forces in all directions # calculate linear solution of stretch 1 resulting from unit load - dl = np.linalg.inv(B)[0, 0] + dl = np.diag(np.linalg.inv(B)) - return dl > 0 + return np.all(dl > 0) def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] diff --git a/src/matadi/_lab_incompressible.py b/src/matadi/_lab_incompressible.py index ef37a50..6ca354b 100644 --- a/src/matadi/_lab_incompressible.py +++ b/src/matadi/_lab_incompressible.py @@ -3,6 +3,9 @@ import matplotlib.pyplot as plt import numpy as np +from ._templates import MaterialHyperelastic +from .math import det + class LabIncompressible: def __init__(self, material, title=None): @@ -47,30 +50,24 @@ def stability(stretch): P = self.material.gradient([F])[0] A = self.material.hessian([F])[0] + d2JdFdF = MaterialHyperelastic(lambda F: det(F)).hessian([F])[0] + + # hydrostatic stress + p = -P[-1, -1] * kinematics(stretch)[-1] + + # convert hessian to (3, 3) matrix + B = np.zeros((3, 3)) + c = [(0, 0), (1, 1), (2, 2)] - # convert hessian to (3, 3) matrix and take (2, 2) submatrix - B = np.zeros((2, 2)) - delta = np.eye(2) - - for i in range(2): - for j in range(2): - B[i, j] = ( - A[i, i, j, j] - - delta[-1, j] / kinematics(stretch)[i] * P[-1, -1] - - kinematics(stretch)[-1] - / kinematics(stretch)[i] - * A[-1, -1, j, j] - + kinematics(stretch)[-1] - / kinematics(stretch)[i] ** 2 - * P[-1, -1] - * delta[i, j] - ) + for i, a in enumerate(c): + for j, b in enumerate(c): + B[i, j] = (A + p * d2JdFdF)[(*a, *b)] # init unit force in direction 1 # calculate linear solution of stretch 1 resulting from unit load - dl = np.linalg.inv(B)[0, 0] + dl = np.diag(np.linalg.inv(B)) - return dl > 0 + return np.all(dl > 0) return ( stress_free(stretch),