Skip to content

Commit e5f27fa

Browse files
adtzlradtzlr
authored
Merge pull request #146 from adtzlr/simplify-stability-check
Simplify the stability checks in `Lab` and `LabIncompressible`
2 parents 633f36b + dfd4359 commit e5f27fa

File tree

5 files changed

+29
-85
lines changed

5 files changed

+29
-85
lines changed

‎README.md‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -186,10 +186,7 @@ fig2, ax2 = lab.plot_shear(data)
186186

187187
![Lab experiments shear(Microsphere)](https://raw.githubusercontent.com/adtzlr/matadi/main/docs/images/plot_shear_lab-microsphere.svg)
188188

189-
Unstable states of deformation can be indicated as dashed lines with the stability argument `lab.plot(data, stability=True)`. This checks whether if
190-
a) the volume ratio is greater zero,
191-
b) the monotonic increasing slope of force per undeformed area vs. stretch and
192-
c) the sign of the resulting stretch from a small superposed force in one direction.
189+
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.
193190

194191
## Hints and usage in FEM modules
195192
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.

‎pyproject.toml‎

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,20 @@
11
[build-system]
2-
requires = ["setuptools>=61"]
2+
requires = ["setuptools>=77.0.3"]
33
build-backend = "setuptools.build_meta"
44

55
[tool.isort]
66
profile = "black"
77

88
[project]
99
name = "matadi"
10+
description = "Material Definition with Automatic Differentiation"
1011
authors = [
11-
{email = "a.dutzler@gmail.com"},
12-
{name = "Andreas Dutzler"}
12+
{name = "Andreas Dutzler", email = "a.dutzler@gmail.com"},
1313
]
14-
description = "Material Definition with Automatic Differentiation"
14+
requires-python = ">=3.9"
1515
readme = "README.md"
16-
license = {file = "LICENSE"}
16+
license = "GPL-3.0-or-later"
17+
license-files = ["LICENSE"]
1718
keywords = [
1819
"python",
1920
"constitution",
@@ -29,7 +30,6 @@ classifiers = [
2930
"Development Status :: 5 - Production/Stable",
3031
"Programming Language :: Python",
3132
"Intended Audience :: Science/Research",
32-
"License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
3333
"Operating System :: OS Independent",
3434
"Programming Language :: Python",
3535
"Programming Language :: Python :: 3",
@@ -43,7 +43,7 @@ classifiers = [
4343
"Topic :: Utilities"
4444
]
4545
dynamic = ["version"]
46-
requires-python = ">=3.9"
46+
4747
dependencies = ["casadi", "numpy"]
4848

4949
[project.optional-dependencies]
@@ -53,5 +53,7 @@ all = ["matplotlib", "scipy"]
5353
version = {attr = "matadi.__about__.__version__"}
5454

5555
[project.urls]
56-
Code = "https://github.com/adtzlr/matadi"
57-
Issues = "https://github.com/adtzlr/matadi/issues"
56+
Homepage = "https://github.com/adtzlr/matadi"
57+
Documentation = "https://github.com/adtzlr/matadi"
58+
Repository = "https://github.com/adtzlr/matadi"
59+
Issues = "https://github.com/adtzlr/matadi/issues"

‎src/matadi/__about__.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
__version__ = "0.2.3"
1+
__version__ = "0.3.0"

‎src/matadi/_lab_compressible.py‎

Lines changed: 14 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,8 @@ def stress(stretch, stretch_2, stretch_3):
2525
F = np.diag([stretch, stretch_2, stretch_3])
2626
return self.material.gradient([F])[0]
2727

28-
def stability(stretch, stretches_23, stretches_23_eps):
28+
def stability(stretch, stretches_23):
2929
F = np.diag([stretch, *stretches_23])
30-
G = np.diag([stretch + 1e-6, *stretches_23_eps])
3130
A = self.material.hessian([F])[0]
3231

3332
# convert hessian to (3, 3) matrix
@@ -38,21 +37,11 @@ def stability(stretch, stretches_23, stretches_23_eps):
3837
for j, b in enumerate(c):
3938
B[i, j] = A[(*a, *b)]
4039

41-
# init unit force in direction 1
42-
df = np.zeros(3)
43-
df[0] = 1
44-
40+
# unit force in direction 1
4541
# calculate linear solution of stretch 1 resulting from unit load
46-
dl = (np.linalg.inv(B) @ df)[0]
47-
48-
# check volume ratio
49-
J = stretch * stretches_23[0] * stretches_23[1]
50-
51-
# check slope of force
52-
Q = self.material.gradient([G])[0][0, 0]
53-
P = self.material.gradient([F])[0][0, 0]
42+
dl = np.linalg.inv(B)[0, 0]
5443

55-
return True if dl > 0 and J > 0 and (Q - P) > 0 else False
44+
return dl > 0
5645

5746
def stress_free(stretches_23, stretch):
5847
s = stress(stretch, *stretches_23)
@@ -62,18 +51,10 @@ def stress_free(stretches_23, stretch):
6251
if not res.success:
6352
res = root(stress_free, np.ones(2) * 1 / np.sqrt(stretch), args=(stretch,))
6453

65-
res_eps = root(stress_free, np.ones(2), args=(stretch + 1e-6,))
66-
if not res_eps.success:
67-
res_eps = root(
68-
stress_free,
69-
np.ones(2) * 1 / np.sqrt(stretch + 1e-6),
70-
args=(stretch + 1e-6,),
71-
)
72-
7354
return (
7455
stress(stretch, *res.x)[0, 0],
7556
*res.x,
76-
stability(stretch, res.x, res_eps.x),
57+
stability(stretch, res.x),
7758
)
7859

7960
def _biaxial(self, stretch):
@@ -84,9 +65,8 @@ def stress(stretch, stretch_3):
8465
F = np.diag([stretch, stretch, stretch_3])
8566
return self.material.gradient([F])[0]
8667

87-
def stability(stretch, stretch_3, stretch_3_eps):
68+
def stability(stretch, stretch_3):
8869
F = np.diag([stretch, stretch, stretch_3])
89-
G = np.diag([stretch + 1e-6, stretch + 1e-6, stretch_3_eps])
9070
A = self.material.hessian([F])[0]
9171

9272
# convert hessian to (3, 3) matrix
@@ -97,36 +77,23 @@ def stability(stretch, stretch_3, stretch_3_eps):
9777
for j, b in enumerate(c):
9878
B[i, j] = A[(*a, *b)]
9979

100-
# init unit force in direction 1
101-
df = np.zeros(3)
102-
df[0] = 1
103-
80+
# unit force in direction 1
10481
# calculate linear solution of stretch 1 resulting from unit load
105-
dl = (np.linalg.inv(B) @ df)[0]
82+
dl = np.linalg.inv(B)[0, 0]
10683

107-
# check volume ratio
108-
J = stretch**2 * stretch_3
109-
110-
# check slope of force
111-
Q = self.material.gradient([G])[0][0, 0]
112-
P = self.material.gradient([F])[0][0, 0]
113-
114-
return True if dl > 0 and J > 0 and (Q - P) > 0 else False
84+
return dl > 0
11585

11686
def stress_free(stretch_3, stretch):
11787
return [stress(stretch, *stretch_3)[2, 2]]
11888

11989
res = root(stress_free, np.ones(1), args=(stretch,))
12090
stretch_3 = res.x[0]
12191

122-
res_eps = root(stress_free, np.ones(1), args=(stretch + 1e-6,))
123-
stretch_3_eps = res_eps.x[0]
124-
12592
return (
12693
stress(stretch, stretch_3)[0, 0],
12794
stretch,
12895
stretch_3,
129-
stability(stretch, stretch_3, stretch_3_eps),
96+
stability(stretch, stretch_3),
13097
)
13198

13299
def _planar(self, stretch):
@@ -137,9 +104,8 @@ def stress(stretch, stretch_3):
137104
F = np.diag([stretch, 1, stretch_3])
138105
return self.material.gradient([F])[0]
139106

140-
def stability(stretch, stretch_3, stretch_3_eps):
107+
def stability(stretch, stretch_3):
141108
F = np.diag([stretch, 1, stretch_3])
142-
G = np.diag([stretch + 1e-6, 1, stretch_3_eps])
143109
A = self.material.hessian([F])[0]
144110

145111
# convert hessian to (3, 3) matrix
@@ -151,35 +117,22 @@ def stability(stretch, stretch_3, stretch_3_eps):
151117
B[i, j] = A[(*a, *b)]
152118

153119
# init unit force in direction 1
154-
df = np.zeros(3)
155-
df[0] = 1
156-
157120
# calculate linear solution of stretch 1 resulting from unit load
158-
dl = (np.linalg.inv(B) @ df)[0]
121+
dl = np.linalg.inv(B)[0, 0]
159122

160-
# check volume ratio
161-
J = stretch * stretch_3
162-
163-
# check slope of force
164-
Q = self.material.gradient([G])[0][0, 0]
165-
P = self.material.gradient([F])[0][0, 0]
166-
167-
return True if dl > 0 and J > 0 and (Q - P) > 0 else False
123+
return dl > 0
168124

169125
def stress_free(stretch_3, stretch):
170126
return [stress(stretch, *stretch_3)[2, 2]]
171127

172128
res = root(stress_free, np.ones(1), args=(stretch,))
173129
stretch_3 = res.x[0]
174130

175-
res_eps = root(stress_free, np.ones(1), args=(stretch + 1e-6,))
176-
stretch_3_eps = res_eps.x[0]
177-
178131
return (
179132
stress(stretch, stretch_3)[0, 0],
180133
1,
181134
stretch_3,
182-
stability(stretch, stretch_3, stretch_3_eps),
135+
stability(stretch, stretch_3),
183136
)
184137

185138
def _shear(self, shear):

‎src/matadi/_lab_incompressible.py‎

Lines changed: 2 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,6 @@ def stress_free(stretch):
4444

4545
def stability(stretch):
4646
F = np.diag([*kinematics(stretch)])
47-
G = np.diag([*kinematics(stretch + 1e-6)])
4847

4948
P = self.material.gradient([F])[0]
5049
A = self.material.hessian([F])[0]
@@ -68,17 +67,10 @@ def stability(stretch):
6867
)
6968

7069
# init unit force in direction 1
71-
df = np.zeros(2)
72-
df[0] = 1
73-
7470
# calculate linear solution of stretch 1 resulting from unit load
75-
dl = (np.linalg.inv(B) @ df)[0]
76-
77-
# check slope of force
78-
Q = self.material.gradient([G])[0][0, 0]
79-
P = self.material.gradient([F])[0][0, 0]
71+
dl = np.linalg.inv(B)[0, 0]
8072

81-
return True if dl > 0 and (Q - P) > 0 else False
73+
return dl > 0
8274

8375
return (
8476
stress_free(stretch),

0 commit comments

Comments
 (0)