From 7d1f59a047737053152b9eeb911ba52ebbef2012 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Thu, 27 Feb 2025 21:45:59 -0500 Subject: [PATCH] switching constitutive model interface to use grad_u in all kinematic variable interfaces. --- pancax/constitutive_models/base.py | 56 +++++---- pancax/constitutive_models/blatz_ko.py | 6 +- pancax/constitutive_models/gent.py | 7 +- pancax/constitutive_models/neohookean.py | 8 +- pancax/constitutive_models/swanson.py | 13 +- pancax/physics_kernels/solid_mechanics.py | 6 +- pyproject.toml | 1 + .../test_base_constitutive_model.py | 13 +- test/constitutive_models/test_gent.py | 19 +-- test/constitutive_models/test_neohookean.py | 19 +-- test/constitutive_models/test_swanson.py | 18 +-- test/constitutive_models/test_swanson4.py | 112 ----------------- test/constitutive_models/test_swanson6.py | 119 ------------------ 13 files changed, 88 insertions(+), 309 deletions(-) delete mode 100644 test/constitutive_models/test_swanson4.py delete mode 100644 test/constitutive_models/test_swanson6.py diff --git a/pancax/constitutive_models/base.py b/pancax/constitutive_models/base.py index 9dc9bf3..a7dd787 100644 --- a/pancax/constitutive_models/base.py +++ b/pancax/constitutive_models/base.py @@ -7,83 +7,97 @@ Scalar = float +State = Float[Array, "ns"] Tensor = Float[Array, "3 3"] class BaseConstitutiveModel(eqx.Module): - def cauchy_stress(self, F: Tensor) -> Tensor: - # F = grad_u + jnp.eye(3) - J = self.jacobian(F) - P = self.pk1_stress(F) + def cauchy_stress(self, grad_u: Tensor) -> Tensor: + F = grad_u + jnp.eye(3) + J = self.jacobian(grad_u) + P = self.pk1_stress(grad_u) return (1. / J) * P @ F.T + def deformation_gradient(self, grad_u: Tensor) -> Tensor: + F = grad_u + jnp.eye(3) + return F + @abstractmethod - def energy(self, F: Tensor) -> Scalar: + def energy(self, grad_u: Tensor) -> Scalar: """ This method returns the algorithmic strain energy density. """ pass - def invariants(self, F: Tensor) -> Tuple[Scalar, Scalar, Scalar]: - I1 = self.I1(F) - I2 = self.I2(F) - I3 = self.jacobian(F)**2 + def invariants(self, grad_u: Tensor) -> Tuple[Scalar, Scalar, Scalar]: + I1 = self.I1(grad_u) + I2 = self.I2(grad_u) + I3 = self.I3(grad_u) return jnp.array([I1, I2, I3]) - def I1(self, F: Tensor) -> Scalar: + def I1(self, grad_u: Tensor) -> Scalar: r""" Calculates the first invariant - - **F**: the deformation gradient + - **grad_u**: the displacement gradient $$ I_1 = tr\left(\mathbf{F}^T\mathbf{F}\right) $$ """ + F = self.deformation_gradient(grad_u) I1 = jnp.trace(F @ F.T) return I1 - def I1_bar(self, F: Tensor) -> Scalar: + def I1_bar(self, grad_u: Tensor) -> Scalar: r""" Calculates the first distortional invariant - - **F**: the deformation gradient + - **grad_u**: the displacement gradient $$ \bar{I}_1 = J^{-2/3}tr\left(\mathbf{F}^T\mathbf{F}\right) $$ """ + F = self.deformation_gradient(grad_u) I1 = jnp.trace(F @ F.T) - J = self.jacobian(F) + J = self.jacobian(grad_u) return jnp.power(J, -2. / 3.) * I1 - def I2(self, F: Tensor) -> Scalar: + def I2(self, grad_u: Tensor) -> Scalar: + F = self.deformation_gradient(grad_u) C = F.T @ F C2 = C @ C I1 = jnp.trace(C) I2 = 0.5 * (I1**2 - jnp.trace(C2)) return I2 - def I2_bar(self, F: Tensor) -> Scalar: + def I2_bar(self, grad_u: Tensor) -> Scalar: + F = self.deformation_gradient(grad_u) C = F.T @ F C2 = C @ C I1 = jnp.trace(C) I2 = 0.5 * (I1**2 - jnp.trace(C2)) - J = self.jacobian(F) + J = self.jacobian(grad_u) return jnp.power(J, -4. / 3.) * I2 - def jacobian(self, F: Tensor) -> Scalar: + def I3(self, grad_u: Tensor) -> Scalar: + J = self.jacobian(grad_u) + return J * J + + def jacobian(self, grad_u: Tensor) -> Scalar: r""" This simply calculate the jacobian but with guard rails to return nonsensical numbers if a non-positive jacobian is encountered during training. - - **F**: the deformation gradient + - **grad_u**: the displacement gradient $$ J = det(\mathbf{F}) $$ """ + F = self.deformation_gradient(grad_u) J = jnp.linalg.det(F) J = jax.lax.cond( J <= 0.0, @@ -93,8 +107,8 @@ def jacobian(self, F: Tensor) -> Scalar: ) return J - def pk1_stress(self, F: Tensor) -> Tensor: - return jax.grad(self.energy, argnums=0)(F) + def pk1_stress(self, grad_u: Tensor) -> Tensor: + return jax.grad(self.energy, argnums=0)(grad_u) def properties(self): return self.__dataclass_fields__ diff --git a/pancax/constitutive_models/blatz_ko.py b/pancax/constitutive_models/blatz_ko.py index bed8502..a965f5a 100644 --- a/pancax/constitutive_models/blatz_ko.py +++ b/pancax/constitutive_models/blatz_ko.py @@ -6,13 +6,13 @@ class BlatzKo(BaseConstitutiveModel): shear_modulus: Property - def energy(self, F): + def energy(self, grad_u): # unpack properties G = self.shear_modulus # kinematics - I2 = self.I2(F) - I3 = jnp.linalg.det(F.T @ F) + I2 = self.I2(grad_u) + I3 = self.I3(grad_u) # constitutive W = (G / 2.) * (I2 / I3 + 2 * jnp.sqrt(I3) - 5.) diff --git a/pancax/constitutive_models/gent.py b/pancax/constitutive_models/gent.py index 9fce7bb..6d0d0a7 100644 --- a/pancax/constitutive_models/gent.py +++ b/pancax/constitutive_models/gent.py @@ -16,14 +16,13 @@ class Gent(BaseConstitutiveModel): shear_modulus: Property Jm_parameter: Property - def energy(self, F): + def energy(self, grad_u): # unpack properties K, G, Jm = self.bulk_modulus, self.shear_modulus, self.Jm_parameter # kinematics - C = F.T @ F - J = self.jacobian(F) - I_1_bar = jnp.trace(jnp.power(J, -2. / 3.) * C) + J = self.jacobian(grad_u) + I_1_bar = self.I1_bar(grad_u) # guard rail check_value = I_1_bar > Jm + 3.0 - 0.001 diff --git a/pancax/constitutive_models/neohookean.py b/pancax/constitutive_models/neohookean.py index 077bca7..9a575c8 100644 --- a/pancax/constitutive_models/neohookean.py +++ b/pancax/constitutive_models/neohookean.py @@ -14,14 +14,12 @@ class NeoHookean(BaseConstitutiveModel): bulk_modulus: Property shear_modulus: Property - def energy(self, F: Tensor) -> Scalar: + def energy(self, grad_u: Tensor) -> Scalar: K, G = self.bulk_modulus, self.shear_modulus # kinematics - # F = jnp.eye(3) + grad_u - C = F.T @ F - J = self.jacobian(F) - I_1_bar = jnp.trace(1. / jnp.square(jnp.cbrt(J)) * C) + J = self.jacobian(grad_u) + I_1_bar = self.I1_bar(grad_u) # constitutive W_vol = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J)) diff --git a/pancax/constitutive_models/swanson.py b/pancax/constitutive_models/swanson.py index de54c54..878abc8 100644 --- a/pancax/constitutive_models/swanson.py +++ b/pancax/constitutive_models/swanson.py @@ -24,7 +24,7 @@ class Swanson(BaseConstitutiveModel): # hack because Swanson is a stupid model cutoff_strain: float = eqx.field(static=True) - def energy(self, F): + def energy(self, grad_u): K = self.bulk_modulus A1, P1 = self.A1, self.P1 B1, Q1 = self.B1, self.Q1 @@ -32,21 +32,24 @@ def energy(self, F): tau_cutoff = (1. / 3.) * (3. + self.cutoff_strain**2) - 1. # kinematics - J = self.jacobian(F) - C = F.T @ F - C_bar = jnp.power(J, -2. / 3.) * C - I_1_bar = jnp.trace(C_bar) + J = self.jacobian(grad_u) + I_1_bar = self.I1_bar(grad_u) + I_2_bar = self.I2_bar(grad_u) tau_1 = (1. / 3.) * I_1_bar - 1. + tau_2 = (1. / 3.) * I_2_bar - 1. tau_tilde_1 = tau_1 + tau_cutoff + tau_tilde_2 = tau_2 + tau_cutoff # constitutive W_vol = K * (J * jnp.log(J) - J + 1.) W_dev_tau = 3. / 2. * ( A1 / (P1 + 1.) * (tau_tilde_1**(P1 + 1.)) + + B1 / (Q1 + 1.) * (tau_tilde_2**(Q1 + 1.)) + C1 / (R1 + 1.) * (tau_tilde_1**(R1 + 1.)) ) W_dev_cutoff = 3. / 2. * ( A1 / (P1 + 1.) * (tau_cutoff**(P1 + 1.)) + + B1 / (Q1 + 1.) * (tau_cutoff**(Q1 + 1.)) + C1 / (R1 + 1.) * (tau_cutoff**(R1 + 1.)) ) W_dev = W_dev_tau - W_dev_cutoff diff --git a/pancax/physics_kernels/solid_mechanics.py b/pancax/physics_kernels/solid_mechanics.py index 6fbe205..64dae9b 100644 --- a/pancax/physics_kernels/solid_mechanics.py +++ b/pancax/physics_kernels/solid_mechanics.py @@ -68,9 +68,5 @@ def __init__(self, constitutive_model, formulation) -> None: self.formulation = formulation def energy(self, params, x, t, u, grad_u, *args): - # _, model = params grad_u = self.formulation.modify_field_gradient(grad_u) - F = grad_u + jnp.eye(3) - return self.constitutive_model.energy(F) - # return self.constitutive_model.energy(grad_u) - # return model.energy(grad_u) + return self.constitutive_model.energy(grad_u) diff --git a/pyproject.toml b/pyproject.toml index 474055c..aa3be20 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,6 +6,7 @@ authors = [ ] dependencies = [ 'matplotlib', + 'meshio', 'netCDF4', 'pandas', 'scipy', diff --git a/test/constitutive_models/test_base_constitutive_model.py b/test/constitutive_models/test_base_constitutive_model.py index f62f6f2..f83ce3f 100644 --- a/test/constitutive_models/test_base_constitutive_model.py +++ b/test/constitutive_models/test_base_constitutive_model.py @@ -1,6 +1,4 @@ -# from pancax import NeoHookean, NeoHookeanFixedBulkModulus from pancax import NeoHookean -# import jax import jax.numpy as jnp @@ -15,7 +13,8 @@ def test_jacobian(): [0., 2., 0.], [0., 0., 1.] ]) - J = model.jacobian(F) + grad_u = F - jnp.eye(3) + J = model.jacobian(grad_u) assert jnp.array_equal(J, jnp.linalg.det(F)) # TODO add better test. @@ -31,10 +30,6 @@ def test_jacobian_bad_value(): [0., 2., 0.], [0., 0., -1.] ]) - J = model.jacobian(F) + grad_u = F - jnp.eye(3) + J = model.jacobian(grad_u) assert jnp.array_equal(J, 1.e3) - - -# def test_bulk_modulus_init(): - # # model = NeoHookeanFixedBulkModulus(bulk_modulus=10.0) - # assert model.bulk_modulus == 10.0 diff --git a/test/constitutive_models/test_gent.py b/test/constitutive_models/test_gent.py index 1671f13..551f831 100644 --- a/test/constitutive_models/test_gent.py +++ b/test/constitutive_models/test_gent.py @@ -1,4 +1,3 @@ -# from pancax import FixedProperties, Gent, GentFixedBulkModulus from pancax import BoundedProperty, Gent from .utils import * import jax @@ -33,10 +32,11 @@ def gent_2(): def simple_shear_test(model): gammas = jnp.linspace(0.0, 1., 100) Fs = jax.vmap(simple_shear)(gammas) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, gamma, I1_bar, J) in zip(psis, sigmas, gammas, I1_bars, Js): psi_an = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J)) + \ @@ -61,10 +61,11 @@ def simple_shear_test(model): def uniaxial_strain_test(model): lambdas = jnp.linspace(1., 2., 100) Fs = jax.vmap(uniaxial_strain)(lambdas) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, lambda_, I1_bar, J) in zip(psis, sigmas, lambdas, I1_bars, Js): psi_an = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J)) + \ diff --git a/test/constitutive_models/test_neohookean.py b/test/constitutive_models/test_neohookean.py index 2e73139..4ac5305 100644 --- a/test/constitutive_models/test_neohookean.py +++ b/test/constitutive_models/test_neohookean.py @@ -29,10 +29,11 @@ def neohookean_2(): def simple_shear_test(model): gammas = jnp.linspace(0.0, 1., 100) Fs = jax.vmap(simple_shear)(gammas) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs)#, props()) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs)#, props()) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, gamma, I1_bar, J) in zip(psis, sigmas, gammas, I1_bars, Js): psi_an = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J)) + \ @@ -57,11 +58,11 @@ def simple_shear_test(model): def uniaxial_strain_test(model): lambdas = jnp.linspace(1., 4., 100) Fs = jax.vmap(uniaxial_strain)(lambdas) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs)#, props()) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs)#, props()) - # K, G = model.unpack_properties(props()) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, lambda_, I1_bar, J) in zip(psis, sigmas, lambdas, I1_bars, Js): psi_an = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J)) + \ diff --git a/test/constitutive_models/test_swanson.py b/test/constitutive_models/test_swanson.py index 0eba1e0..29c26aa 100644 --- a/test/constitutive_models/test_swanson.py +++ b/test/constitutive_models/test_swanson.py @@ -46,11 +46,12 @@ def swanson_2(): def simple_shear_test(model): gammas = jnp.linspace(0.05, 1., 100) Fs = jax.vmap(simple_shear)(gammas) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, I1_bar, J, B_bar) in zip(psis, sigmas, I1_bars, Js, B_bars): psi_an = K * (J * jnp.log(J) - J + 1.) + \ @@ -70,11 +71,12 @@ def simple_shear_test(model): def uniaxial_strain_test(model): lambdas = jnp.linspace(1.2, 4., 100) Fs = jax.vmap(uniaxial_strain)(lambdas) + grad_us = jax.vmap(lambda F: F - jnp.eye(3))(Fs) B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) - Js = jax.vmap(model.jacobian)(Fs) - I1_bars = jax.vmap(model.I1_bar)(Fs) - psis = jax.vmap(model.energy, in_axes=(0,))(Fs) - sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(Fs) + Js = jax.vmap(model.jacobian)(grad_us) + I1_bars = jax.vmap(model.I1_bar)(grad_us) + psis = jax.vmap(model.energy, in_axes=(0,))(grad_us) + sigmas = jax.vmap(model.cauchy_stress, in_axes=(0,))(grad_us) for (psi, sigma, I1_bar, J, B_bar) in zip(psis, sigmas, I1_bars, Js, B_bars): psi_an = K * (J * jnp.log(J) - J + 1.) + \ diff --git a/test/constitutive_models/test_swanson4.py b/test/constitutive_models/test_swanson4.py deleted file mode 100644 index 2076e71..0000000 --- a/test/constitutive_models/test_swanson4.py +++ /dev/null @@ -1,112 +0,0 @@ -# from pancax import FixedProperties -# from pancax import Swanson4, Swanson4FixedBulkModulus -# from .utils import * -# import jax -# import jax.numpy as jnp -# import pytest - - -# @pytest.fixture -# def swanson(): -# return Swanson4(cutoff_strain=0.01) - - -# @pytest.fixture -# def swanson_fbm(): -# return Swanson4FixedBulkModulus(bulk_modulus=10.0, cutoff_strain=0.01) - - -# @pytest.fixture -# def properties(): -# return FixedProperties([10.0, 0.93074321, -0.07673672, 0.10448312, 1.71691036]) - - -# @pytest.fixture -# def properties_fbm(): -# return FixedProperties([0.93074321, -0.07673672, 0.10448312, 1.71691036]) - - - -# def simple_shear_test(model, props): -# gammas = jnp.linspace(0.05, 1., 100) -# Fs = jax.vmap(simple_shear)(gammas) -# B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) -# Js = jax.vmap(model.jacobian)(Fs) -# I1_bars = jax.vmap(model.I1_bar)(Fs) -# psis = jax.vmap(model.energy, in_axes=(0, None))(Fs, props()) -# sigmas = jax.vmap(model.cauchy_stress, in_axes=(0, None))(Fs, props()) - -# K, A1, P1, C1, R1 = model.unpack_properties(props()) - -# for (psi, sigma, I1_bar, J, B_bar) in zip(psis, sigmas, I1_bars, Js, B_bars): -# psi_an = K * (J * jnp.log(J) - J + 1.) + \ -# 1.5 * A1 / (P1 + 1.) * (I1_bar / 3. - 1.)**(P1 + 1.) + \ -# 1.5 * C1 / (R1 + 1.) * (I1_bar / 3. - 1.)**(R1 + 1.) -# dUdI1_bar = 0.5 * A1 * (I1_bar / 3. - 1.)**P1 + \ -# 0.5 * C1 * (I1_bar / 3. - 1.)**R1 -# B_bar_dev = B_bar - (1. / 3.) * jnp.trace(B_bar) * jnp.eye(3) -# sigma_an = (2. / J) * dUdI1_bar * B_bar_dev + K * jnp.log(J) * jnp.eye(3) - -# assert jnp.allclose(psi, psi_an, atol=1e-3) -# for i in range(3): -# for j in range(3): -# assert jnp.allclose(sigma[i, j], sigma_an[i, j], atol=1e-3) - - -# def uniaxial_strain_test(model, props): -# lambdas = jnp.linspace(1.2, 4., 100) -# Fs = jax.vmap(uniaxial_strain)(lambdas) -# B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) -# Js = jax.vmap(model.jacobian)(Fs) -# I1_bars = jax.vmap(model.I1_bar)(Fs) -# psis = jax.vmap(model.energy, in_axes=(0, None))(Fs, props()) -# sigmas = jax.vmap(model.cauchy_stress, in_axes=(0, None))(Fs, props()) - -# K, A1, P1, C1, R1 = model.unpack_properties(props()) - -# for (psi, sigma, I1_bar, J, B_bar) in zip(psis, sigmas, I1_bars, Js, B_bars): -# psi_an = K * (J * jnp.log(J) - J + 1.) + \ -# 1.5 * A1 / (P1 + 1.) * (I1_bar / 3. - 1.)**(P1 + 1.) + \ -# 1.5 * C1 / (R1 + 1.) * (I1_bar / 3. - 1.)**(R1 + 1.) -# dUdI1_bar = 0.5 * A1 * (I1_bar / 3. - 1.)**P1 + \ -# 0.5 * C1 * (I1_bar / 3. - 1.)**R1 -# B_bar_dev = B_bar - (1. / 3.) * jnp.trace(B_bar) * jnp.eye(3) -# sigma_an = (2. / J) * dUdI1_bar * B_bar_dev + K * jnp.log(J) * jnp.eye(3) - -# assert jnp.allclose(psi, psi_an, atol=1e-3) -# for i in range(3): -# for j in range(3): -# assert jnp.allclose(sigma[i, j], sigma_an[i, j], atol=1e-3) - - -# def unpack_props_test(model, props): -# props = model.unpack_properties(props()) -# assert props[0] == 10.0 -# assert props[1] == 0.93074321 -# assert props[2] == -0.07673672 -# assert props[3] == 0.10448312 -# assert props[4] == 1.71691036 - - -# def test_swanson_unpack_props( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# unpack_props_test(swanson, properties) -# unpack_props_test(swanson_fbm, properties_fbm) - - -# def test_simple_shear( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# simple_shear_test(swanson, properties) -# simple_shear_test(swanson_fbm, properties_fbm) - - -# def test_uniaxial_strain( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# uniaxial_strain_test(swanson, properties) -# uniaxial_strain_test(swanson_fbm, properties_fbm) diff --git a/test/constitutive_models/test_swanson6.py b/test/constitutive_models/test_swanson6.py deleted file mode 100644 index df3e2a0..0000000 --- a/test/constitutive_models/test_swanson6.py +++ /dev/null @@ -1,119 +0,0 @@ -# from pancax import FixedProperties -# from pancax import Swanson6, Swanson6FixedBulkModulus -# from .utils import * -# import jax -# import jax.numpy as jnp -# import pytest - - -# @pytest.fixture -# def swanson(): -# return Swanson6(cutoff_strain=0.01) - - -# @pytest.fixture -# def swanson_fbm(): -# return Swanson6FixedBulkModulus(bulk_modulus=10.0, cutoff_strain=0.01) - - -# @pytest.fixture -# def properties(): -# return FixedProperties([10.0, 0.93074321, -0.07673672, 0.0, 0.5, 0.10448312, 1.71691036]) - - -# @pytest.fixture -# def properties_fbm(): -# return FixedProperties([0.93074321, -0.07673672, 0.0, 0.5, 0.10448312, 1.71691036]) - - -# def simple_shear_test(model, props): -# gammas = jnp.linspace(0.05, 1., 100) -# Fs = jax.vmap(simple_shear)(gammas) -# B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) -# Js = jax.vmap(model.jacobian)(Fs) -# I1_bars = jax.vmap(model.I1_bar)(Fs) -# I2_bars = jax.vmap(model.I2_bar)(Fs) -# psis = jax.vmap(model.energy, in_axes=(0, None))(Fs, props()) -# sigmas = jax.vmap(model.cauchy_stress, in_axes=(0, None))(Fs, props()) - -# K, A1, P1, B1, Q1, C1, R1 = model.unpack_properties(props()) - -# for (psi, sigma, I1_bar, I2_bar, J, B_bar) in zip(psis, sigmas, I1_bars, I2_bars, Js, B_bars): -# psi_an = K * (J * jnp.log(J) - J + 1.) + \ -# 1.5 * A1 / (P1 + 1.) * (I1_bar / 3. - 1.)**(P1 + 1.) + \ -# 1.5 * B1 / (Q1 + 1.) * (I2_bar / 3. - 1.)**(R1 + 1.) + \ -# 1.5 * C1 / (R1 + 1.) * (I1_bar / 3. - 1.)**(R1 + 1.) -# dUdI1_bar = 0.5 * A1 * (I1_bar / 3. - 1.)**P1 + \ -# 0.5 * C1 * (I1_bar / 3. - 1.)**R1 -# dUdI2_bar = 0.5 * B1 * (I2_bar / 3. - 1.)**Q1 -# B_bar_dev = B_bar - (1. / 3.) * jnp.trace(B_bar) * jnp.eye(3) -# sigma_an = (2. / J) * (dUdI1_bar + I1_bar * dUdI2_bar) * B_bar_dev + K * jnp.log(J) * jnp.eye(3) - -# assert jnp.allclose(psi, psi_an, atol=1e-3) -# for i in range(3): -# for j in range(3): -# assert jnp.allclose(sigma[i, j], sigma_an[i, j], atol=1e-3) - - -# def uniaxial_strain_test(model, props): -# lambdas = jnp.linspace(1.2, 4., 100) -# Fs = jax.vmap(uniaxial_strain)(lambdas) -# B_bars = jax.vmap(lambda F: jnp.power(jnp.linalg.det(F), -2. / 3.) * F @ F.T)(Fs) -# Js = jax.vmap(model.jacobian)(Fs) -# I1_bars = jax.vmap(model.I1_bar)(Fs) -# I2_bars = jax.vmap(model.I2_bar)(Fs) -# psis = jax.vmap(model.energy, in_axes=(0, None))(Fs, props()) -# sigmas = jax.vmap(model.cauchy_stress, in_axes=(0, None))(Fs, props()) - -# K, A1, P1, B1, Q1, C1, R1 = model.unpack_properties(props()) - -# for (psi, sigma, I1_bar, I2_bar, J, B_bar) in zip(psis, sigmas, I1_bars, I2_bars, Js, B_bars): -# psi_an = K * (J * jnp.log(J) - J + 1.) + \ -# 1.5 * A1 / (P1 + 1.) * (I1_bar / 3. - 1.)**(P1 + 1.) + \ -# 1.5 * B1 / (Q1 + 1.) * (I2_bar / 3. - 1.)**(R1 + 1.) + \ -# 1.5 * C1 / (R1 + 1.) * (I1_bar / 3. - 1.)**(R1 + 1.) -# dUdI1_bar = 0.5 * A1 * (I1_bar / 3. - 1.)**P1 + \ -# 0.5 * C1 * (I1_bar / 3. - 1.)**R1 -# dUdI2_bar = 0.5 * B1 * (I2_bar / 3. - 1.)**Q1 -# B_bar_dev = B_bar - (1. / 3.) * jnp.trace(B_bar) * jnp.eye(3) -# sigma_an = (2. / J) * (dUdI1_bar + I1_bar * dUdI2_bar) * B_bar_dev + K * jnp.log(J) * jnp.eye(3) - -# assert jnp.allclose(psi, psi_an, atol=1e-3) -# for i in range(3): -# for j in range(3): -# assert jnp.allclose(sigma[i, j], sigma_an[i, j], atol=1e-3) - - -# def unpack_props_test(model, props): -# props = model.unpack_properties(props()) -# assert props[0] == 10.0 -# assert props[1] == 0.93074321 -# assert props[2] == -0.07673672 -# assert props[3] == 0.0 -# assert props[4] == 0.5 -# assert props[5] == 0.10448312 -# assert props[6] == 1.71691036 - - -# def test_swanson_unpack_props( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# unpack_props_test(swanson, properties) -# unpack_props_test(swanson_fbm, properties_fbm) - - -# def test_simple_shear( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# simple_shear_test(swanson, properties) -# simple_shear_test(swanson_fbm, properties_fbm) - - -# def test_uniaxial_strain( -# swanson, properties, -# swanson_fbm, properties_fbm -# ): -# uniaxial_strain_test(swanson, properties) -# uniaxial_strain_test(swanson_fbm, properties_fbm)