diff --git a/README.md b/README.md index 0ced2a9..f4b62e3 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ Material Definition with Automatic Differentiation (AD) [![PyPI version shields.io](https://img.shields.io/pypi/v/matadi.svg)](https://pypi.python.org/pypi/matadi/) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) ![Made with love in Graz (Austria)](https://img.shields.io/badge/Made%20with%20%E2%9D%A4%EF%B8%8F%20in-Graz%20(Austria)-0c674a) [![codecov](https://codecov.io/gh/adtzlr/matadi/branch/main/graph/badge.svg?token=2EY2U4ZL35)](https://codecov.io/gh/adtzlr/matadi) [![DOI](https://zenodo.org/badge/408564756.svg)](https://zenodo.org/badge/latestdoi/408564756) ![Codestyle black](https://img.shields.io/badge/code%20style-black-black) ![GitHub Repo stars](https://img.shields.io/github/stars/adtzlr/matadi?logo=github) ![PyPI - Downloads](https://img.shields.io/pypi/dm/matadi) -matADi is a simple Python module which acts as a wrapper on top of [casADi](https://web.casadi.org/) for easy definitions of hyperelastic strain energy functions. Gradients (stresses) and hessians (elasticity tensors) are carried out by casADi's powerful and fast **Automatic Differentiation (AD)** capabilities. It is designed to handle inputs with trailing axes which is especially useful for the application in Python-based finite element modules like [scikit-fem](https://scikit-fem.readthedocs.io/en/latest/) or [FElupe](https://adtzlr.github.io/felupe/). Mixed-field formulations are supported as well as single-field formulations. +matADi is a simple Python module which acts as a wrapper on top of [casADi](https://web.casadi.org/) [1] for easy definitions of hyperelastic strain energy functions. Gradients (stresses) and hessians (elasticity tensors) are carried out by casADi's powerful and fast **Automatic Differentiation (AD)** capabilities. It is designed to handle inputs with trailing axes which is especially useful for the application in Python-based finite element modules like [scikit-fem](https://scikit-fem.readthedocs.io/en/latest/) or [FElupe](https://adtzlr.github.io/felupe/). Mixed-field formulations are supported as well as single-field formulations. ## Installation Install `matADi` from PyPI via pip. @@ -114,10 +114,18 @@ from matadi import Lab lab = Lab(NH) data = lab.run(ux=True, bx=True, ps=True) -fix, ax = lab.plot(data) +fig, ax = lab.plot(data) ``` ![Lab experiments(Neo-Hooke)](https://raw.githubusercontent.com/adtzlr/matadi/main/docs/images/plot_lab-nh.svg) +Unstable states of deformation can be indicated as dashed lines with the stability argument `lab.plot(data, stability=True)`. This checks if +a) the volume ratio is greater zero, +b) the slope of stress vs. stretch and +c) the sign of the resulting stretch from a small superposed force in one direction. + ## Hints -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). \ No newline at end of file +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). + +## References +[1] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, *CasADi - A software framework for nonlinear optimization and optimal control*, Math. Prog. Comp., vol. 11, no. 1, pp. 1–36, 2019, [![DOI:10.1007/s12532-018-0139-4](https://zenodo.org/badge/10.1007/s12532-018-0139-4.svg)](https://doi.org/10.1007/s12532-018-0139-4) \ No newline at end of file diff --git a/matadi/__about__.py b/matadi/__about__.py index a73339b..00ec2dc 100644 --- a/matadi/__about__.py +++ b/matadi/__about__.py @@ -1 +1 @@ -__version__ = "0.0.8" +__version__ = "0.0.9" diff --git a/matadi/_lab.py b/matadi/_lab.py index a76e660..b9f451a 100644 --- a/matadi/_lab.py +++ b/matadi/_lab.py @@ -16,79 +16,206 @@ def _get_title(self): ) def _uniaxial(self, stretch): - def stress(stretch_2, stretch_3): + def stress(stretch, stretch_2, stretch_3): F = np.diag([stretch, stretch_2, stretch_3]) return self.material.gradient([F])[0] - def stress_free(stretches_23): - return [stress(*stretches_23)[1, 1], stress(*stretches_23)[2, 2]] + def stability(stretch, stretches_23, stretches_23_eps): + F = np.diag([stretch, *stretches_23]) + G = np.diag([stretch + 1e-6, *stretches_23_eps]) + A = self.material.hessian([F])[0] - res = root(stress_free, np.ones(2)) - stretch_2, stretch_3 = res.x + # convert hessian to (6,6) matrix + B = np.zeros((6, 6)) + c = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (2, 0)] - return stress(stretch_2, stretch_3)[0, 0], stretch_2, stretch_3 + for i, a in enumerate(c): + for j, b in enumerate(c): + B[i, j] = A[(*a, *b)] + + # init unit force in direction 1 + df = np.zeros(6) + df[0] = 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] + + return True if dl > 0 and J > 0 and (Q - P) > 0 else False + + def stress_free(stretches_23, stretch): + s = stress(stretch, *stretches_23) + return [s[1, 1], s[2, 2]] + + res = root(stress_free, np.ones(2), args=(stretch,)) + res_eps = root(stress_free, np.ones(2), args=(stretch + 1e-6,)) + + return ( + stress(stretch, *res.x)[0, 0], + *res.x, + stability(stretch, res.x, res_eps.x), + ) def _biaxial(self, stretch): - def stress(stretch_3): + def stress(stretch, stretch_3): F = np.diag([stretch, stretch, stretch_3]) return self.material.gradient([F])[0] - def stress_free(stretch_3): - return [stress(*stretch_3)[2, 2]] + def stability(stretch, stretch_3, stretch_3_eps): + 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 (6,6) matrix + B = np.zeros((6, 6)) + c = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (2, 0)] - res = root(stress_free, np.ones(1)) + for i, a in enumerate(c): + for j, b in enumerate(c): + B[i, j] = A[(*a, *b)] + + # init unit force in direction 1 + df = np.zeros(6) + df[0] = 1 + + # calculate linear solution of stretch 1 resulting from unit load + dl = (np.linalg.inv(B) @ df)[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 + + def stress_free(stretch_3, stretch): + return [stress(stretch, *stretch_3)[2, 2]] + + res = root(stress_free, np.ones(1), args=(stretch,)) stretch_3 = res.x[0] - return stress(stretch_3)[0, 0], stretch, stretch_3 + 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), + ) def _planar(self, stretch): - def stress(stretch_3): + def stress(stretch, stretch_3): F = np.diag([stretch, 1, stretch_3]) return self.material.gradient([F])[0] - def stress_free(stretch_3): - return [stress(*stretch_3)[2, 2]] + def stability(stretch, stretch_3, stretch_3_eps): + F = np.diag([stretch, 1, stretch_3]) + G = np.diag([stretch + 1e-6, 1, stretch_3_eps]) + A = self.material.hessian([F])[0] - res = root(stress_free, np.ones(1)) + # convert hessian to (6,6) matrix + B = np.zeros((6, 6)) + c = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (2, 0)] + + for i, a in enumerate(c): + for j, b in enumerate(c): + B[i, j] = A[(*a, *b)] + + # init unit force in direction 1 + df = np.zeros(6) + df[0] = 1 + + # calculate linear solution of stretch 1 resulting from unit load + dl = (np.linalg.inv(B) @ df)[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 + + def stress_free(stretch_3, stretch): + return [stress(stretch, *stretch_3)[2, 2]] + + res = root(stress_free, np.ones(1), args=(stretch,)) stretch_3 = res.x[0] - return stress(stretch_3)[0, 0], 1, stretch_3 + 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), + ) def run(self, ux=True, bx=True, ps=True, stretch_max=2.5, num=50): out = [] - Data = namedtuple("Data", "label stretch stretch_2 stretch_3 stress") + Data = namedtuple("Data", "label stretch stretch_2 stretch_3 stress stability") if ux: stretch_ux = np.linspace(1 - (stretch_max - 1) / 5, stretch_max, num) - stress_ux, stretch_2_ux, stretch_3_ux = np.array( + stress_ux, stretch_2_ux, stretch_3_ux, stability_ux = np.array( [self._uniaxial(s11) for s11 in stretch_ux] ).T ux_data = Data( - "uniaxial", stretch_ux, stretch_2_ux, stretch_3_ux, stress_ux + "uniaxial", + stretch_ux, + stretch_2_ux, + stretch_3_ux, + stress_ux, + stability_ux, ) out.append(ux_data) if bx: stretch_bx = np.linspace(1, (stretch_max - 1) / 2 + 1, num) - stress_bx, stretch_2_bx, stretch_3_bx = np.array( + stress_bx, stretch_2_bx, stretch_3_bx, stability_bx = np.array( [self._biaxial(s11) for s11 in stretch_bx] ).T - bx_data = Data("biaxial", stretch_bx, stretch_2_bx, stretch_3_bx, stress_bx) + bx_data = Data( + "biaxial", + stretch_bx, + stretch_2_bx, + stretch_3_bx, + stress_bx, + stability_bx, + ) out.append(bx_data) if ps: stretch_ps = np.linspace(1, stretch_max, num) - stress_ps, stretch_2_ps, stretch_3_ps = np.array( + stress_ps, stretch_2_ps, stretch_3_ps, stability_ps = np.array( [self._planar(s11) for s11 in stretch_ps] ).T - ps_data = Data("planar", stretch_ps, stretch_2_ps, stretch_3_ps, stress_ps) + ps_data = Data( + "planar", + stretch_ps, + stretch_2_ps, + stretch_3_ps, + stress_ps, + stability_ps, + ) out.append(ps_data) return out - def plot(self, data): + def plot(self, data, stability=False): fig, ax = plt.subplots() @@ -99,7 +226,28 @@ def plot(self, data): } for d in data: - ax.plot(d.stretch, d.stress, **lineargs[d.label], label=d.label) + + if stability: + + stable = np.array(d.stability, dtype=bool) + + unstable_idx = np.arange(len(d.stretch))[~stable] + stable_idx = np.arange(len(d.stretch))[stable] + + stress_stable = d.stress.copy() + stress_stable[unstable_idx] = np.nan + + stress_unstable = d.stress.copy() + stress_unstable[stable_idx] = np.nan + + ax.plot(d.stretch, stress_stable, **lineargs[d.label], label=d.label) + ax.plot( + d.stretch, stress_unstable, **lineargs[d.label], linestyle="--", + ) + + else: + + ax.plot(d.stretch, d.stress, **lineargs[d.label], label=d.label) ax.grid() ax.set_title(self.title) diff --git a/matadi/models/_isotropic_hyperelasticity.py b/matadi/models/_isotropic_hyperelasticity.py index e1d9fcd..e00f3de 100644 --- a/matadi/models/_isotropic_hyperelasticity.py +++ b/matadi/models/_isotropic_hyperelasticity.py @@ -12,7 +12,7 @@ def mooney_rivlin(F, C10, C01, bulk): J = det(F) C = transpose(F) @ F I1 = J ** (-2 / 3) * trace(C) - I2 = J ** (-4 / 3) * trace(C @ C) + I2 = J ** (-4 / 3) * (trace(C) ** 2 - trace(C @ C)) / 2 return C10 * (I1 - 3) + C01 * (I2 - 3) + bulk * (J - 1) ** 2 / 2 @@ -32,7 +32,7 @@ def third_order_deformation(F, C10, C01, C11, C20, C30, bulk): J = det(F) C = transpose(F) @ F I1 = J ** (-2 / 3) * trace(C) - I2 = J ** (-4 / 3) * trace(C @ C) + I2 = J ** (-4 / 3) * (trace(C) ** 2 - trace(C @ C)) / 2 return ( C10 * (I1 - 3) + C01 * (I2 - 3) diff --git a/pyproject.toml b/pyproject.toml index b7fb381..9435f98 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "matadi" -version = "0.0.8" +version = "0.0.9" description = "Material Definition with Automatic Differentiation" readme = "README.md" requires-python = ">=3.6" diff --git a/setup.cfg b/setup.cfg index 8ec329c..25357da 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = matadi -version = 0.0.8 +version = 0.0.9 author = Andreas Dutzler author_email = a.dutzler@gmail.com description = Material Definition with Automatic Differentiation diff --git a/tests/test_lab.py b/tests/test_lab.py index 750a23b..72d0b6f 100644 --- a/tests/test_lab.py +++ b/tests/test_lab.py @@ -6,11 +6,7 @@ def test_lab(): - mat = MaterialHyperelastic( - neo_hooke, - C10=0.5, - bulk=20.0, - ) + mat = MaterialHyperelastic(neo_hooke, C10=0.5, bulk=5000.0,) lab = Lab(mat) data = lab.run(ux=False, bx=False, ps=False, num=20) @@ -19,7 +15,7 @@ def test_lab(): data = lab.run(ux=True, bx=True, ps=True, num=20) fig, ax = lab.plot(data) - #plt.close(fig) + plt.close(fig) # dW and DW are always lists... assert len(data[0].stress) == 20 diff --git a/tests/test_lab_stability.py b/tests/test_lab_stability.py new file mode 100644 index 0000000..27f61ae --- /dev/null +++ b/tests/test_lab_stability.py @@ -0,0 +1,23 @@ +from matadi import MaterialHyperelastic, Lab +from matadi.models import mooney_rivlin + +import matplotlib.pyplot as plt + + +def test_lab(): + + mat = MaterialHyperelastic(mooney_rivlin, C10=0.5, C01=0.5, bulk=5000.0,) + + lab = Lab(mat) + data = lab.run(ux=True, bx=True, ps=True, num=50) + + fig, ax = lab.plot(data, stability=True) + plt.close(fig) + + # dW and DW are always lists... + assert len(data[0].stress) == 50 + assert len(data[0].stretch) == 50 + + +if __name__ == "__main__": + test_lab()