From 1a4052f32400b4ff7b59447bbc6ba607b5c08339 Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 20:06:43 +0200 Subject: [PATCH 01/10] add hgo model --- matadi/models/_anisotropic_hyperelasticity.py | 46 +++++++++++++++++-- 1 file changed, 43 insertions(+), 3 deletions(-) diff --git a/matadi/models/_anisotropic_hyperelasticity.py b/matadi/models/_anisotropic_hyperelasticity.py index 96bc746..53850f1 100644 --- a/matadi/models/_anisotropic_hyperelasticity.py +++ b/matadi/models/_anisotropic_hyperelasticity.py @@ -1,7 +1,7 @@ -from ..math import transpose, det, DM, cos, sin, pi, sqrt, if_else +from ..math import transpose, det, DM, exp, cos, sin, pi, sqrt, if_else -def fiber(F, E, k=1, angle=30, axis=2, compression=False): +def fiber(F, E, angle, k=1, axis=2, compression=False): a = angle * pi / 180 plane = [(1, 2), (2, 0), (0, 1)][axis] @@ -21,9 +21,49 @@ def fiber(F, E, k=1, angle=30, axis=2, compression=False): return E * strain ** 2 / 2 -def fiber_family(F, E, k=1, angle=30, axis=2, compression=False): +def fiber_family(F, E, angle, k=1, axis=2, compression=False): f1 = fiber(F, E, k, angle, axis, compression) f2 = fiber(F, E, k, -angle, axis, compression) return f1 + f2 + + +def holzapfel_gasser_ogden(F, c, k1, k2, angle, axis=2): + + C = transpose(F) @ F + J1, J2, J3 = invariants(C) + + Cu = J3 ** (-1 / 3) * C + + I1 = J3 ** (-1 / 3) * J1 + # I2 = J3 ** (-2 / 3) * J2 + + alpha = angle * pi / 180 + plane = [(1, 2), (2, 0), (0, 1)][axis] + + N1 = DM.zeros(3) + N1[plane, :] = DM([cos(alpha), sin(alpha)]) + + N2 = DM.zeros(3) + N2[plane, :] = DM([cos(alpha), -sin(alpha)]) + + A1 = N1 @ transpose(N1) + A2 = N2 @ transpose(N2) + + I4 = trace(Cu @ A1) + I6 = trace(Cu @ A2) + + # I5 = trace(C @ C @ A1) + # I7 = trace(C @ C @ A2) + + # I8 = (transpose(N1) @ N2) * transpose(N1) @ C @ N2 + + W_iso = c / 2 * (I1 - 3) + + w1 = exp(k2 * (I4 - 1) ** 2) - 1 + w2 = exp(k2 * (I6 - 1) ** 2) - 1 + + W_aniso = k1 / (2 * k2) * (w1 + w2) + + return W_iso + W_aniso From 88ef1968c3782275b4955713a98562f9fb39989b Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 20:08:36 +0200 Subject: [PATCH 02/10] Update _anisotropic_hyperelasticity.py --- matadi/models/_anisotropic_hyperelasticity.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/matadi/models/_anisotropic_hyperelasticity.py b/matadi/models/_anisotropic_hyperelasticity.py index 53850f1..f22e4db 100644 --- a/matadi/models/_anisotropic_hyperelasticity.py +++ b/matadi/models/_anisotropic_hyperelasticity.py @@ -33,11 +33,7 @@ def holzapfel_gasser_ogden(F, c, k1, k2, angle, axis=2): C = transpose(F) @ F J1, J2, J3 = invariants(C) - - Cu = J3 ** (-1 / 3) * C - I1 = J3 ** (-1 / 3) * J1 - # I2 = J3 ** (-2 / 3) * J2 alpha = angle * pi / 180 plane = [(1, 2), (2, 0), (0, 1)][axis] @@ -51,13 +47,8 @@ def holzapfel_gasser_ogden(F, c, k1, k2, angle, axis=2): A1 = N1 @ transpose(N1) A2 = N2 @ transpose(N2) - I4 = trace(Cu @ A1) - I6 = trace(Cu @ A2) - - # I5 = trace(C @ C @ A1) - # I7 = trace(C @ C @ A2) - - # I8 = (transpose(N1) @ N2) * transpose(N1) @ C @ N2 + I4 = trace(J3 ** (-1 / 3) * C @ A1) + I6 = trace(J3 ** (-1 / 3) * C @ A2) W_iso = c / 2 * (I1 - 3) From 16d2e155da000d056b9b6879c9c1a1e5c4e7ebcd Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 20:09:15 +0200 Subject: [PATCH 03/10] Update math.py --- matadi/math.py | 1 + 1 file changed, 1 insertion(+) diff --git a/matadi/math.py b/matadi/math.py index a887cac..a7359eb 100644 --- a/matadi/math.py +++ b/matadi/math.py @@ -25,6 +25,7 @@ asinh, acosh, # math + exp, sqrt, sum1, sum2, From b41d77007ee5eca03b985919dd05f85ceeb8c3cb Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 20:10:10 +0200 Subject: [PATCH 04/10] Update __init__.py --- matadi/models/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/matadi/models/__init__.py b/matadi/models/__init__.py index 4d2eac4..f553ab2 100644 --- a/matadi/models/__init__.py +++ b/matadi/models/__init__.py @@ -10,4 +10,5 @@ from ._anisotropic_hyperelasticity import ( fiber, fiber_family, + holzapfel_gasser_ogden, ) From 553ed5a8f6f42805d2499497eb0da62ac7ce17f4 Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 20:11:09 +0200 Subject: [PATCH 05/10] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 6c0b467..f88a7d3 100644 --- a/README.md +++ b/README.md @@ -96,6 +96,7 @@ Available isotropic hyperelastic material models: Available anisotropic hyperelastic material models: - Fiber - Fiber-family (+/- combination of single Fiber) +- [Holzapfel Gasser Ogden](https://link.springer.com/article/10.1023/A:1010835316564) Any user-defined isotropic hyperelastic strain energy density function may be passed as the `fun` argument of `MaterialHyperelastic` by using the following template: From 6da3a55006a6b4a0dc8ea519c7e41627cc151202 Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 22:54:33 +0200 Subject: [PATCH 06/10] add missing imports --- matadi/models/_anisotropic_hyperelasticity.py | 25 +++++++++++++++---- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/matadi/models/_anisotropic_hyperelasticity.py b/matadi/models/_anisotropic_hyperelasticity.py index f22e4db..ba23b45 100644 --- a/matadi/models/_anisotropic_hyperelasticity.py +++ b/matadi/models/_anisotropic_hyperelasticity.py @@ -1,7 +1,20 @@ -from ..math import transpose, det, DM, exp, cos, sin, pi, sqrt, if_else +from ..math import ( + transpose, + det, + DM, + exp, + cos, + sin, + pi, + sqrt, + if_else, + invariants, + trace, +) def fiber(F, E, angle, k=1, axis=2, compression=False): + "Fiber" a = angle * pi / 180 plane = [(1, 2), (2, 0), (0, 1)][axis] @@ -22,6 +35,7 @@ def fiber(F, E, angle, k=1, axis=2, compression=False): def fiber_family(F, E, angle, k=1, axis=2, compression=False): + "Fiber Family" f1 = fiber(F, E, k, angle, axis, compression) f2 = fiber(F, E, k, -angle, axis, compression) @@ -29,7 +43,8 @@ def fiber_family(F, E, angle, k=1, axis=2, compression=False): return f1 + f2 -def holzapfel_gasser_ogden(F, c, k1, k2, angle, axis=2): +def holzapfel_gasser_ogden(F, c, k1, k2, kappa, angle, axis=2): + "Holzapfel-Gasser-Ogden" C = transpose(F) @ F J1, J2, J3 = invariants(C) @@ -52,9 +67,9 @@ def holzapfel_gasser_ogden(F, c, k1, k2, angle, axis=2): W_iso = c / 2 * (I1 - 3) - w1 = exp(k2 * (I4 - 1) ** 2) - 1 - w2 = exp(k2 * (I6 - 1) ** 2) - 1 + w4 = exp(k2 * (kappa * I1 + (1 - 3 * kappa) * I4 - 1) ** 2) - 1 + w6 = exp(k2 * (kappa * I1 + (1 - 3 * kappa) * I6 - 1) ** 2) - 1 - W_aniso = k1 / (2 * k2) * (w1 + w2) + W_aniso = k1 / (2 * k2) * (w4 + w6) return W_iso + W_aniso From f0cb49647e0d1bd89d3c1d545807058956ffdeac Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 22:54:49 +0200 Subject: [PATCH 07/10] update tests --- tests/test_fiber.py | 50 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 3 deletions(-) diff --git a/tests/test_fiber.py b/tests/test_fiber.py index ce981a6..12c6b3a 100644 --- a/tests/test_fiber.py +++ b/tests/test_fiber.py @@ -1,7 +1,7 @@ import numpy as np from matadi import Variable, MaterialHyperelastic -from matadi.models import fiber, fiber_family +from matadi.models import fiber, fiber_family, holzapfel_gasser_ogden from matadi.math import det, transpose, trace, invariants, sqrt @@ -22,9 +22,53 @@ def test_fiber(): DW = M.hessian([FF]) assert W0[0].shape == (2,) - assert dW[0].shape == (3, 3, 2,) - assert DW[0].shape == (3, 3, 3, 3, 2,) + assert dW[0].shape == ( + 3, + 3, + 2, + ) + assert DW[0].shape == ( + 3, + 3, + 3, + 3, + 2, + ) + + +def test_hgo(): + + # data + FF = np.zeros((3, 3, 2)) + for a in range(3): + FF[a, a] += 1 + + for model in [holzapfel_gasser_ogden]: + + # init Material + M = MaterialHyperelastic( + model, c=0.0764, k1=996.6, k2=524.6, kappa=0.2, angle=49.98, axis=2 + ) + + W0 = M.function([FF]) + dW = M.gradient([FF]) + DW = M.hessian([FF]) + + assert W0[0].shape == (2,) + assert dW[0].shape == ( + 3, + 3, + 2, + ) + assert DW[0].shape == ( + 3, + 3, + 3, + 3, + 2, + ) if __name__ == "__main__": test_fiber() + test_hgo() From 975a793a1258e816461fee44f11f43f7bfb79593 Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 22:54:57 +0200 Subject: [PATCH 08/10] lint black --- tests/test_simple.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/test_simple.py b/tests/test_simple.py index 7789605..0467dce 100644 --- a/tests/test_simple.py +++ b/tests/test_simple.py @@ -6,7 +6,7 @@ def neohooke(x, mu=1.0, bulk=200.0): """Strain energy density function of nearly-incompressible - Neo-Hookean isotropic hyperelastic material formulation.""" + Neo-Hookean isotropic hyperelastic material formulation.""" F = x[0] C = transpose(F) @ F @@ -30,7 +30,11 @@ def test_simple(): FF[a, a] += 1 # init Material - W = Material(x=[F], fun=neohooke, kwargs={"mu": 1.0, "bulk": 10.0},) + W = Material( + x=[F], + fun=neohooke, + kwargs={"mu": 1.0, "bulk": 10.0}, + ) W0 = W.function([FF]) dW = W.gradient([FF]) From d64eb5c9e850e74bfb39e8c08be70b00dc7af4d7 Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 22:58:16 +0200 Subject: [PATCH 09/10] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f88a7d3..ddb2f8b 100644 --- a/README.md +++ b/README.md @@ -96,7 +96,7 @@ Available isotropic hyperelastic material models: Available anisotropic hyperelastic material models: - Fiber - Fiber-family (+/- combination of single Fiber) -- [Holzapfel Gasser Ogden](https://link.springer.com/article/10.1023/A:1010835316564) +- [Holzapfel Gasser Ogden](https://royalsocietypublishing.org/doi/full/10.1098/rsif.2005.0073) Any user-defined isotropic hyperelastic strain energy density function may be passed as the `fun` argument of `MaterialHyperelastic` by using the following template: From 14a41ee2eb87f51f5cbb0880b646b27bf1cc076d Mon Sep 17 00:00:00 2001 From: Andreas D Date: Thu, 23 Sep 2021 23:07:18 +0200 Subject: [PATCH 10/10] bump version --- matadi/__about__.py | 2 +- setup.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/matadi/__about__.py b/matadi/__about__.py index 034f46c..6526deb 100644 --- a/matadi/__about__.py +++ b/matadi/__about__.py @@ -1 +1 @@ -__version__ = "0.0.6" +__version__ = "0.0.7" diff --git a/setup.cfg b/setup.cfg index f56fa06..a9c149f 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = matadi -version = 0.0.6 +version = 0.0.7 author = Andreas Dutzler author_email = a.dutzler@gmail.com description = Material Definition with Automatic Differentiation