Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/tests/test_tangential.py
24902e40d876maxwell_viscoelastic…

Last Change
Raw File

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()
c4science · Help