Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88629037
test_tangential.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Oct 19, 20:39
Size
3 KB
Mime Type
text/x-python
Expires
Mon, Oct 21, 20:39 (2 d)
Engine
blob
Format
Raw Data
Handle
21482009
Attached To
rTAMAAS tamaas
test_tangential.py
View Options
# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
# Copyright (©) 2020-2023 Lucas Frérot
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from
__future__
import
division
,
print_function
import
numpy
as
np
from
numpy.linalg
import
norm
import
tamaas
as
tm
def
test_pressure
():
E
=
1.0
nu
=
0.5
mu
=
0.5
n
=
81
L
=
1.0
p_target
=
np
.
array
([
0.5e-4
,
0
,
2e-4
],
dtype
=
tm
.
dtype
)
# g_target = np.array([-0.000223, 0, 9.99911], dtype=tm.dtype)
x
=
np
.
linspace
(
-
L
/
2
,
L
/
2
,
n
,
dtype
=
tm
.
dtype
)
y
=
np
.
copy
(
x
)
x
,
y
=
np
.
meshgrid
(
x
,
y
,
indexing
=
"ij"
)
r_sqrd
=
(
x
**
2
+
y
**
2
)
# Spherical contact
R
=
10.0
surface
=
R
*
np
.
sqrt
((
r_sqrd
<
R
**
2
)
*
(
1
-
r_sqrd
/
R
**
2
))
# Creating model
model
=
tm
.
ModelFactory
.
createModel
(
tm
.
model_type
.
surface_2d
,
[
L
,
L
],
[
n
,
n
])
model
.
E
,
model
.
nu
=
E
,
nu
p
=
model
[
'traction'
]
# Theoretical solution
Estar
=
E
/
(
1.0
-
nu
**
2
)
Fn
=
p_target
[
2
]
*
L
**
2
Fx
=
p_target
[
0
]
*
L
**
2
d_theory
=
np
.
power
(
3
/
4
*
Fn
/
Estar
/
np
.
sqrt
(
R
),
2
/
3
)
a_theory
=
np
.
sqrt
(
R
*
d_theory
)
c_theory
=
a_theory
*
(
1
-
Fx
/
(
mu
*
Fn
))
**
(
1
/
3
)
p0_theory
=
np
.
power
(
6
*
Fn
*
Estar
**
2
/
np
.
pi
**
3
/
R
**
2
,
1
/
3
)
t1_theory
=
mu
*
p0_theory
t2_theory
=
t1_theory
*
c_theory
/
a_theory
# p_theory = p0_theory * np.sqrt((r_sqrd < a_theory**2)
# * (1 - r_sqrd / a_theory**2))
t_theory
=
t1_theory
*
np
.
sqrt
((
r_sqrd
<
a_theory
**
2
)
*
(
1
-
r_sqrd
/
a_theory
**
2
))
t_theory
=
t_theory
-
t2_theory
*
np
.
sqrt
((
r_sqrd
<
c_theory
**
2
)
*
(
1
-
r_sqrd
/
c_theory
**
2
))
def
assert_error
(
err
):
error
=
norm
(
p
[
...
,
0
]
-
t_theory
)
/
norm
(
t_theory
)
print
(
error
)
assert
error
<
err
# Test Kato solver
solver
=
tm
.
Kato
(
model
,
surface
,
1e-12
,
mu
)
solver
.
max_iter
=
1000
solver
.
solve
(
p_target
,
100
)
assert_error
(
4e-2
)
# solver.solveRelaxed(g_target)
# assert_error(1e-2)
# # Test BeckTeboulle solver
# solver = tm.BeckTeboulle(model, surface, 1e-12, mu)
# solver.solve(g_target)
# assert_error(1e-2)
# Test Condat solver
solver
=
tm
.
Condat
(
model
,
surface
,
1e-12
,
mu
)
solver
.
max_iter
=
5000
# or 10000
solver
.
solve
(
p_target
)
assert_error
(
7e-2
)
# 4e-2 for 10000 iterations
# Test tangential Polonsky Keer solver
solver
=
tm
.
PolonskyKeerTan
(
model
,
surface
,
1e-12
,
mu
)
solver
.
max_iter
=
1000
solver
.
solve
(
p_target
)
assert_error
(
4e-2
)
if
__name__
==
"__main__"
:
test_pressure
()
Event Timeline
Log In to Comment