Skip to content

Commit 7a4e2ac

Browse files
finsbergfinsberg
committed
Fix tests
1 parent 1ec5d85 commit 7a4e2ac

8 files changed

Lines changed: 121 additions & 82 deletions

File tree

‎src/pulse/cardiac_model.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,8 @@ def strain_energy(
6161
if C_dot is not None:
6262
psi += self.viscoelasticity.strain_energy(C_dot)
6363

64+
return psi
65+
6466
def S(
6567
self,
6668
C: ufl.core.expr.Expr,

‎src/pulse/material_models/guccione.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@ class Guccione(HyperElasticMaterial):
6363
bf: Variable = field(default_factory=lambda: Variable(8.0, "dimensionless"))
6464
bt: Variable = field(default_factory=lambda: Variable(2.0, "dimensionless"))
6565
bfs: Variable = field(default_factory=lambda: Variable(4.0, "dimensionless"))
66+
deviatoric: bool = True
6667

6768
def __post_init__(self):
6869
# Check that all values are positive
@@ -109,10 +110,11 @@ def is_isotropic(self) -> bool:
109110
return False
110111

111112
def _Q(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
112-
J = ufl.sqrt(ufl.det(C))
113113
dim = C.ufl_shape[0]
114-
C_dev = pow(J, -2.0 / dim) * C
115-
E = 0.5 * (C_dev - ufl.Identity(dim))
114+
if self.deviatoric:
115+
Jm23 = pow(ufl.det(C), -1.0 / dim)
116+
C *= Jm23 # Make C deviatoric
117+
E = 0.5 * (C - ufl.Identity(dim))
116118

117119
bt = self.bt.to_base_units()
118120
bf = self.bf.to_base_units()

‎src/pulse/material_models/holzapfelogden.py‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -107,6 +107,7 @@ class HolzapfelOgden(HyperElasticMaterial):
107107
b_fs: Variable = field(default_factory=lambda: Variable(0.0, "dimensionless"))
108108
use_subplus: bool = field(default=True, repr=False)
109109
use_heaviside: bool = field(default=True, repr=False)
110+
deviatoric: bool = field(default=True, repr=False)
110111

111112
_W1_func: Invariant = field(
112113
init=False,
@@ -305,7 +306,15 @@ def _W8fs(self, I8fs: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
305306

306307
def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
307308
dim = C.ufl_shape[0]
308-
I1 = pow(invariants.I3(C), -1 / dim) * invariants.I1(C)
309+
310+
I1 = invariants.I1(C)
311+
312+
if self.deviatoric:
313+
Jm23 = pow(ufl.det(C), -1.0 / dim)
314+
C *= Jm23 # Make C deviatoric
315+
Jm23 = pow(invariants.I3(C), -1 / dim)
316+
I1 *= Jm23 # Convert to deviatoric I1 - see https://arxiv.org/pdf/2009.08754
317+
309318
I4f = self._I4f(C)
310319
I4s = self._I4s(C)
311320
I8fs = self._I8fs(C)

‎tests/test_active_model.py‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,9 @@ def test_transversely_active_stress(eta, Ta, mesh, u) -> None:
2727

2828
u.interpolate(lambda x: x)
2929
F = kinematics.DeformationGradient(u)
30-
W = active_model.strain_energy(F)
30+
C = F.T * F
31+
32+
W = active_model.strain_energy(C)
3133

3234
assert active_model.Fe(F) is F
3335
active_model.activation.assign(Ta)
@@ -53,7 +55,8 @@ def test_Passive(u) -> None:
5355

5456
u.interpolate(lambda x: x)
5557
F = kinematics.DeformationGradient(u)
58+
C = F.T * F
5659

5760
assert active_model.Fe(F) is F
58-
W = active_model.strain_energy(F)
61+
W = active_model.strain_energy(C)
5962
assert np.isclose(dolfinx.fem.assemble_scalar(dolfinx.fem.form(W * ufl.dx)), 0.0)

‎tests/test_cardiac_model.py‎

Lines changed: 54 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -11,51 +11,84 @@
1111

1212
@pytest.mark.parametrize("isotropy", (pulse.active_stress.ActiveStressModels.transversely,))
1313
@pytest.mark.parametrize(
14-
"comp_model",
15-
(pulse.compressibility.Incompressible(), pulse.compressibility.Compressible()),
14+
"comp_model_cls",
15+
(pulse.compressibility.Incompressible, pulse.compressibility.Compressible),
1616
)
17-
def test_CardiacModel_HolzapfelOgden(comp_model, isotropy, mesh, u):
18-
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
17+
def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u):
18+
# material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
1919
f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
2020
s0 = dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0))
21-
material = pulse.HolzapfelOgden(f0=f0, s0=s0, **material_params)
22-
21+
material = pulse.HolzapfelOgden(
22+
f0=f0,
23+
s0=s0,
24+
a=1.0,
25+
b=0.0,
26+
a_f=1.0,
27+
b_f=0.0,
28+
a_fs=0.0,
29+
b_fs=0.0,
30+
deviatoric=False,
31+
)
32+
comp_model = comp_model_cls()
2333
active_model = pulse.ActiveStress(f0, isotropy=isotropy)
2434
model = pulse.CardiacModel(
2535
material=material,
2636
active=active_model,
2737
compressibility=comp_model,
2838
)
39+
comp_model.register(p=1000.0)
2940
u.interpolate(lambda x: x / 10.0)
3041
F = pulse.kinematics.DeformationGradient(u)
31-
psi = model.strain_energy(F, p=1.0)
42+
C = F.T * F
43+
psi = model.strain_energy(C)
3244
value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
3345

46+
# value_mat = dolfinx.fem.assemble_scalar(dolfinx.fem.form(material.strain_energy(C) * ufl.dx))
47+
# value_comp = dolfinx.fem.assemble_scalar(
48+
# dolfinx.fem.form(comp_model.strain_energy(C) * ufl.dx),
49+
# )
50+
# value_active = dolfinx.fem.assemble_scalar(
51+
# dolfinx.fem.form(active_model.strain_energy(C) * ufl.dx)
52+
# )
53+
54+
# F = I + 0.1 I, C = 1.21 I, I4f = 1.21
55+
# J = det(F) = 1.1 ** 3
56+
3457
if isinstance(comp_model, pulse.compressibility.Incompressible):
35-
assert math.isclose(value, 53647.14346074856)
58+
# psi = 0.5 * a * (I1 - 3) + 0.5 * a_f * (I4f - 1)**2 + p (J - 1)
59+
# psi = 0.5 * 1000*1 * (3 * 1.21 - 3) + 0.5 * 1000.0 * 1.0 * (1.21 - 1)**2 +
60+
# 1000.0 * (1.1 ** 3 - 1) = 668.05
61+
assert math.isclose(value, 668.05)
3662
else:
37-
assert math.isclose(value, 103220.36041941364)
63+
# psi = 0.5 * a * (I1 - 3) + 0.5 * a_f * (I4f - 1)**2 + kappa * (J * ln(J) - J + 1)
64+
# psi = 0.5 * 1000*1 * (3 * 1.21 - 3) + 0.5 * 1000.0 * 1.0 * (1.21 - 1)**2 +
65+
# 1e6 * (1.1 ** 3 * math.log(1.1 ** 3) - 1.1 ** 3 + 1) = 49910.5979586692
66+
assert math.isclose(value, 49910.5979586692)
3867

3968

4069
@pytest.mark.parametrize("isotropy", (pulse.active_stress.ActiveStressModels.transversely,))
4170
@pytest.mark.parametrize(
42-
"comp_model",
43-
(pulse.compressibility.Incompressible(), pulse.compressibility.Compressible()),
71+
"comp_model_cls",
72+
(pulse.compressibility.Incompressible, pulse.compressibility.Compressible),
4473
)
45-
def test_CardiacModel_NeoHookean(comp_model, isotropy, mesh, u):
74+
def test_CardiacModel_NeoHookean(comp_model_cls, isotropy, mesh, u):
4675
material = pulse.NeoHookean(
4776
mu=dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(15.0)),
4877
)
4978
f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
79+
comp_model = comp_model_cls()
5080
active_model = pulse.ActiveStress(f0, isotropy=isotropy)
5181
model = pulse.CardiacModel(
5282
material=material,
5383
active=active_model,
5484
compressibility=comp_model,
5585
)
86+
comp_model.register(p=dolfinx.fem.Constant(mesh, 1.0))
5687
u.interpolate(lambda x: x / 10.0)
5788
F = pulse.kinematics.DeformationGradient(u)
58-
psi = model.strain_energy(F, p=1.0)
89+
C = F.T * F
90+
psi = model.strain_energy(C)
91+
5992
value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
6093

6194
if isinstance(comp_model, pulse.compressibility.Incompressible):
@@ -66,24 +99,27 @@ def test_CardiacModel_NeoHookean(comp_model, isotropy, mesh, u):
6699

67100
@pytest.mark.parametrize("isotropy", (pulse.active_stress.ActiveStressModels.transversely,))
68101
@pytest.mark.parametrize(
69-
"comp_model",
70-
(pulse.compressibility.Incompressible(), pulse.compressibility.Compressible()),
102+
"comp_model_cls",
103+
(pulse.compressibility.Incompressible, pulse.compressibility.Compressible),
71104
)
72-
def test_CardiacModel_Guccione(comp_model, isotropy, mesh, u):
105+
def test_CardiacModel_Guccione(comp_model_cls, isotropy, mesh, u):
73106
f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
74107
s0 = dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0))
75108
n0 = dolfinx.fem.Constant(mesh, (0.0, 0.0, 1.0))
76109
material_params = pulse.Guccione.default_parameters()
77-
material = pulse.Guccione(f0=f0, s0=s0, n0=n0, **material_params)
110+
material = pulse.Guccione(f0=f0, s0=s0, n0=n0, **material_params, deviatoric=False)
78111
active_model = pulse.ActiveStress(f0, isotropy=isotropy)
112+
comp_model = comp_model_cls()
79113
model = pulse.CardiacModel(
80114
material=material,
81115
active=active_model,
82116
compressibility=comp_model,
83117
)
118+
comp_model.register(p=dolfinx.fem.Constant(mesh, 1.0))
84119
u.interpolate(lambda x: x / 10.0)
85120
F = pulse.kinematics.DeformationGradient(u)
86-
psi = model.strain_energy(F, p=1.0)
121+
C = F.T * F
122+
psi = model.strain_energy(C)
87123
value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
88124

89125
if isinstance(comp_model, pulse.compressibility.Incompressible):

‎tests/test_compressibility.py‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -12,29 +12,29 @@ def test_Incompressible(u, P1) -> None:
1212
p.x.array[:] = 3.14
1313
u.interpolate(lambda x: x)
1414
F = pulse.kinematics.DeformationGradient(u)
15-
J = pulse.kinematics.Jacobian(F)
15+
C = F.T * F
1616
comp = pulse.compressibility.Incompressible()
1717
comp.register(p)
18-
psi = comp.strain_energy(J)
18+
psi = comp.strain_energy(C)
1919
value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
2020
assert math.isclose(value, 3.14 * (8 - 1))
2121

2222

2323
def test_Incompressible_with_missing_pressure_raises_MissingModelAttribute(u) -> None:
2424
u.interpolate(lambda x: x)
2525
F = pulse.kinematics.DeformationGradient(u)
26-
J = pulse.kinematics.Jacobian(F)
26+
C = F.T * F
2727
comp = pulse.compressibility.Incompressible()
2828
with pytest.raises(pulse.exceptions.MissingModelAttribute):
29-
comp.strain_energy(J)
29+
comp.strain_energy(C)
3030

3131

3232
def test_Compressible(u) -> None:
3333
u.interpolate(lambda x: x)
3434
F = pulse.kinematics.DeformationGradient(u)
35-
J = pulse.kinematics.Jacobian(F)
35+
C = F.T * F
3636
kappa = 1234
3737
comp = pulse.compressibility.Compressible(kappa=pulse.Variable(kappa, "Pa"))
38-
psi = comp.strain_energy(J)
38+
psi = comp.strain_energy(C)
3939
value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
4040
assert math.isclose(value, kappa * (8 * math.log(8) - 8 + 1))

‎tests/test_invariants.py‎

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,8 @@
1616
def test_I1(cls, expected, u) -> None:
1717
u.interpolate(lambda x: x)
1818
F = cls(u)
19-
I1 = invariants.I1(F)
19+
C = F.T * F
20+
I1 = invariants.I1(C)
2021

2122
assert np.isclose(
2223
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I1 * ufl.dx)),
@@ -37,7 +38,8 @@ def test_I1(cls, expected, u) -> None:
3738
def test_I2(cls, expected, u) -> None:
3839
u.interpolate(lambda x: x)
3940
F = cls(u)
40-
I2 = invariants.I2(F)
41+
C = F.T * F
42+
I2 = invariants.I2(C)
4143
assert np.isclose(
4244
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I2 * ufl.dx)),
4345
expected,
@@ -54,7 +56,8 @@ def test_I2(cls, expected, u) -> None:
5456
def test_I3(cls, expected, u) -> None:
5557
u.interpolate(lambda x: x)
5658
F = cls(u)
57-
I3 = invariants.I3(F)
59+
C = F.T * F
60+
I3 = invariants.I3(C)
5861

5962
assert np.isclose(
6063
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I3 * ufl.dx)),
@@ -72,8 +75,9 @@ def test_I3(cls, expected, u) -> None:
7275
def test_I4(cls, expected, u, mesh) -> None:
7376
u.interpolate(lambda x: x)
7477
F = cls(u)
78+
C = F.T * F
7579
a0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
76-
I4 = invariants.I4(F, a0)
80+
I4 = invariants.I4(C, a0)
7781

7882
assert np.isclose(
7983
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I4 * ufl.dx)),
@@ -91,8 +95,9 @@ def test_I4(cls, expected, u, mesh) -> None:
9195
def test_I5(cls, expected, u, mesh) -> None:
9296
u.interpolate(lambda x: x)
9397
F = cls(u)
98+
C = F.T * F
9499
a0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
95-
I5 = invariants.I5(F, a0)
100+
I5 = invariants.I5(C, a0)
96101

97102
assert np.isclose(
98103
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I5 * ufl.dx)),
@@ -107,9 +112,10 @@ def test_I5(cls, expected, u, mesh) -> None:
107112
def test_I8(cls, expected, u, mesh) -> None:
108113
u.interpolate(lambda x: x)
109114
F = cls(u)
115+
C = F.T * F
110116
a0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
111117
b0 = dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0))
112-
I8 = invariants.I8(F, a0, b0)
118+
I8 = invariants.I8(C, a0, b0)
113119
assert np.isclose(
114120
dolfinx.fem.assemble_scalar(dolfinx.fem.form(I8 * ufl.dx)),
115121
expected,

0 commit comments

Comments
 (0)