From 1776cdd792ffb5e1530b2b53b485d9f8a93d7a63 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 19 May 2025 23:15:58 +0200 Subject: [PATCH 1/6] Simplify the stability check in `LabCompressible` also, don't check the volume ratio anymore --- README.md | 5 +-- src/matadi/_lab_compressible.py | 75 ++++++--------------------------- 2 files changed, 15 insertions(+), 65 deletions(-) diff --git a/README.md b/README.md index 060f979..5b2f83d 100644 --- a/README.md +++ b/README.md @@ -186,10 +186,7 @@ fig2, ax2 = lab.plot_shear(data) ![Lab experiments shear(Microsphere)](https://raw.githubusercontent.com/adtzlr/matadi/main/docs/images/plot_shear_lab-microsphere.svg) -Unstable states of deformation can be indicated as dashed lines with the stability argument `lab.plot(data, stability=True)`. This checks whether if -a) the volume ratio is greater zero, -b) the monotonic increasing slope of force per undeformed area vs. stretch and -c) the sign of the resulting stretch from a small superposed force in one direction. +Unstable states of deformation can be indicated as dashed lines with the stability argument `lab.plot(data, stability=True)`. This checks whether all incremental stretches due to a small superposed normal force in one direction are positive. ## Hints and usage in FEM modules For tensor-valued material definitions use `MaterialTensor` (e.g. any stress-strain relation). Also, please have a look at [casADi's documentation](https://web.casadi.org/). It is very powerful but unfortunately does not support all the Python stuff you would expect. For example Python's default if-else-statements can't be used in combination with symbolic conditions (use `math.if_else(cond, if_true, if_false)` instead). Contrary to [casADi](https://web.casadi.org/), the gradient of the eigenvalue function is stabilized by a perturbation of the diagonal components. diff --git a/src/matadi/_lab_compressible.py b/src/matadi/_lab_compressible.py index 779e417..7756615 100644 --- a/src/matadi/_lab_compressible.py +++ b/src/matadi/_lab_compressible.py @@ -25,9 +25,8 @@ def stress(stretch, stretch_2, stretch_3): F = np.diag([stretch, stretch_2, stretch_3]) return self.material.gradient([F])[0] - def stability(stretch, stretches_23, stretches_23_eps): + def stability(stretch, stretches_23): F = np.diag([stretch, *stretches_23]) - G = np.diag([stretch + 1e-6, *stretches_23_eps]) A = self.material.hessian([F])[0] # convert hessian to (3, 3) matrix @@ -38,21 +37,11 @@ def stability(stretch, stretches_23, stretches_23_eps): for j, b in enumerate(c): B[i, j] = A[(*a, *b)] - # init unit force in direction 1 - df = np.zeros(3) - df[0] = 1 - + # unit force in direction 1 # calculate linear solution of stretch 1 resulting from unit load - dl = (np.linalg.inv(B) @ df)[0] - - # check volume ratio - J = stretch * stretches_23[0] * stretches_23[1] - - # check slope of force - Q = self.material.gradient([G])[0][0, 0] - P = self.material.gradient([F])[0][0, 0] + dl = np.linalg.inv(B)[0, 0] - return True if dl > 0 and J > 0 and (Q - P) > 0 else False + return np.all(dl > 0) def stress_free(stretches_23, stretch): s = stress(stretch, *stretches_23) @@ -62,18 +51,10 @@ def stress_free(stretches_23, stretch): if not res.success: res = root(stress_free, np.ones(2) * 1 / np.sqrt(stretch), args=(stretch,)) - res_eps = root(stress_free, np.ones(2), args=(stretch + 1e-6,)) - if not res_eps.success: - res_eps = root( - stress_free, - np.ones(2) * 1 / np.sqrt(stretch + 1e-6), - args=(stretch + 1e-6,), - ) - return ( stress(stretch, *res.x)[0, 0], *res.x, - stability(stretch, res.x, res_eps.x), + stability(stretch, res.x), ) def _biaxial(self, stretch): @@ -84,9 +65,8 @@ def stress(stretch, stretch_3): F = np.diag([stretch, stretch, stretch_3]) return self.material.gradient([F])[0] - def stability(stretch, stretch_3, stretch_3_eps): + def stability(stretch, stretch_3): F = np.diag([stretch, stretch, stretch_3]) - G = np.diag([stretch + 1e-6, stretch + 1e-6, stretch_3_eps]) A = self.material.hessian([F])[0] # convert hessian to (3, 3) matrix @@ -97,21 +77,11 @@ def stability(stretch, stretch_3, stretch_3_eps): for j, b in enumerate(c): B[i, j] = A[(*a, *b)] - # init unit force in direction 1 - df = np.zeros(3) - df[0] = 1 - + # unit force in direction 1 # calculate linear solution of stretch 1 resulting from unit load - dl = (np.linalg.inv(B) @ df)[0] + dl = np.linalg.inv(B)[0, 0] - # check volume ratio - J = stretch**2 * stretch_3 - - # check slope of force - Q = self.material.gradient([G])[0][0, 0] - P = self.material.gradient([F])[0][0, 0] - - return True if dl > 0 and J > 0 and (Q - P) > 0 else False + return np.all(dl > 0) def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] @@ -119,14 +89,11 @@ def stress_free(stretch_3, stretch): res = root(stress_free, np.ones(1), args=(stretch,)) stretch_3 = res.x[0] - res_eps = root(stress_free, np.ones(1), args=(stretch + 1e-6,)) - stretch_3_eps = res_eps.x[0] - return ( stress(stretch, stretch_3)[0, 0], stretch, stretch_3, - stability(stretch, stretch_3, stretch_3_eps), + stability(stretch, stretch_3), ) def _planar(self, stretch): @@ -137,9 +104,8 @@ def stress(stretch, stretch_3): F = np.diag([stretch, 1, stretch_3]) return self.material.gradient([F])[0] - def stability(stretch, stretch_3, stretch_3_eps): + def stability(stretch, stretch_3): F = np.diag([stretch, 1, stretch_3]) - G = np.diag([stretch + 1e-6, 1, stretch_3_eps]) A = self.material.hessian([F])[0] # convert hessian to (3, 3) matrix @@ -151,20 +117,10 @@ def stability(stretch, stretch_3, stretch_3_eps): B[i, j] = A[(*a, *b)] # init unit force in direction 1 - df = np.zeros(3) - df[0] = 1 - # calculate linear solution of stretch 1 resulting from unit load - dl = (np.linalg.inv(B) @ df)[0] + dl = np.linalg.inv(B)[0, 0] - # check volume ratio - J = stretch * stretch_3 - - # check slope of force - Q = self.material.gradient([G])[0][0, 0] - P = self.material.gradient([F])[0][0, 0] - - return True if dl > 0 and J > 0 and (Q - P) > 0 else False + return np.all(dl > 0) def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] @@ -172,14 +128,11 @@ def stress_free(stretch_3, stretch): res = root(stress_free, np.ones(1), args=(stretch,)) stretch_3 = res.x[0] - res_eps = root(stress_free, np.ones(1), args=(stretch + 1e-6,)) - stretch_3_eps = res_eps.x[0] - return ( stress(stretch, stretch_3)[0, 0], 1, stretch_3, - stability(stretch, stretch_3, stretch_3_eps), + stability(stretch, stretch_3), ) def _shear(self, shear): From fcdf10d80829a0ec767785e2822e0433bd817116 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 20 May 2025 23:02:25 +0200 Subject: [PATCH 2/6] Update pyproject.toml --- pyproject.toml | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index b3cf2fd..2c4f928 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=61"] +requires = ["setuptools>=77.0.3"] build-backend = "setuptools.build_meta" [tool.isort] @@ -7,13 +7,14 @@ profile = "black" [project] name = "matadi" +description = "Material Definition with Automatic Differentiation" authors = [ - {email = "a.dutzler@gmail.com"}, - {name = "Andreas Dutzler"} + {name = "Andreas Dutzler", email = "a.dutzler@gmail.com"}, ] -description = "Material Definition with Automatic Differentiation" +requires-python = ">=3.9" readme = "README.md" -license = {file = "LICENSE"} +license = "GPL-3.0-or-later" +license-files = ["LICENSE"] keywords = [ "python", "constitution", @@ -43,7 +44,7 @@ classifiers = [ "Topic :: Utilities" ] dynamic = ["version"] -requires-python = ">=3.9" + dependencies = ["casadi", "numpy"] [project.optional-dependencies] @@ -53,5 +54,7 @@ all = ["matplotlib", "scipy"] version = {attr = "matadi.__about__.__version__"} [project.urls] -Code = "https://github.com/adtzlr/matadi" -Issues = "https://github.com/adtzlr/matadi/issues" \ No newline at end of file +Homepage = "https://github.com/adtzlr/matadi" +Documentation = "https://github.com/adtzlr/matadi" +Repository = "https://github.com/adtzlr/matadi" +Issues = "https://github.com/adtzlr/matadi/issues" From 3a5eae45ff0325dc3c1f98ef4c9fa065d5ff141f Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 20 May 2025 23:02:31 +0200 Subject: [PATCH 3/6] Update _lab_compressible.py --- src/matadi/_lab_compressible.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/matadi/_lab_compressible.py b/src/matadi/_lab_compressible.py index 7756615..373a0c9 100644 --- a/src/matadi/_lab_compressible.py +++ b/src/matadi/_lab_compressible.py @@ -41,7 +41,7 @@ def stability(stretch, stretches_23): # calculate linear solution of stretch 1 resulting from unit load dl = np.linalg.inv(B)[0, 0] - return np.all(dl > 0) + return dl > 0 def stress_free(stretches_23, stretch): s = stress(stretch, *stretches_23) @@ -81,7 +81,7 @@ def stability(stretch, stretch_3): # calculate linear solution of stretch 1 resulting from unit load dl = np.linalg.inv(B)[0, 0] - return np.all(dl > 0) + return dl > 0 def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] @@ -120,7 +120,7 @@ def stability(stretch, stretch_3): # calculate linear solution of stretch 1 resulting from unit load dl = np.linalg.inv(B)[0, 0] - return np.all(dl > 0) + return dl > 0 def stress_free(stretch_3, stretch): return [stress(stretch, *stretch_3)[2, 2]] From 20662136f983dfcf488ac35e0e69efc8b8e797d8 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 20 May 2025 23:02:33 +0200 Subject: [PATCH 4/6] Update _lab_incompressible.py --- src/matadi/_lab_incompressible.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/matadi/_lab_incompressible.py b/src/matadi/_lab_incompressible.py index 2111fbf..ef37a50 100644 --- a/src/matadi/_lab_incompressible.py +++ b/src/matadi/_lab_incompressible.py @@ -44,7 +44,6 @@ def stress_free(stretch): def stability(stretch): F = np.diag([*kinematics(stretch)]) - G = np.diag([*kinematics(stretch + 1e-6)]) P = self.material.gradient([F])[0] A = self.material.hessian([F])[0] @@ -68,17 +67,10 @@ def stability(stretch): ) # init unit force in direction 1 - df = np.zeros(2) - df[0] = 1 - # calculate linear solution of stretch 1 resulting from unit load - dl = (np.linalg.inv(B) @ df)[0] - - # check slope of force - Q = self.material.gradient([G])[0][0, 0] - P = self.material.gradient([F])[0][0, 0] + dl = np.linalg.inv(B)[0, 0] - return True if dl > 0 and (Q - P) > 0 else False + return dl > 0 return ( stress_free(stretch), From 63a70c88daca70b6e33190f649c1c6838b5cfe4f Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 20 May 2025 23:05:27 +0200 Subject: [PATCH 5/6] Update pyproject.toml --- pyproject.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 2c4f928..41fa596 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,6 @@ classifiers = [ "Development Status :: 5 - Production/Stable", "Programming Language :: Python", "Intended Audience :: Science/Research", - "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)", "Operating System :: OS Independent", "Programming Language :: Python", "Programming Language :: Python :: 3", From dfd4359ee3460b11b150d059ec471473bbf01df7 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 20 May 2025 23:06:33 +0200 Subject: [PATCH 6/6] Update __about__.py --- src/matadi/__about__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/matadi/__about__.py b/src/matadi/__about__.py index d31c31e..493f741 100644 --- a/src/matadi/__about__.py +++ b/src/matadi/__about__.py @@ -1 +1 @@ -__version__ = "0.2.3" +__version__ = "0.3.0"