Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_tangential.py
c6be98825283features/netcdf-mpi

Last Change
Raw File

test_tangential.py
View Options

#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @author Son Pham-Ba <son.phamba@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Tamaas 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Tamaas. If not, see <http://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():
tm.initialize()
E = 1.0
nu = 0.5
mu = 0.5
n = 81
L = 1.0
p_target = np.array([0.5e-4, 0, 2e-4]);
g_target = np.array([-0.000223, 0, 9.99911]);
x = np.linspace(-L/2, L/2, n)
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.setElasticity(E, nu)
p = model.getTraction();
# 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.setMaxIterations(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.setMaxIterations(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.setMaxIterations(1000)
solver.solve(p_target)
assert_error(4e-2)
tm.finalize()
if __name__ == "__main__":
test_pressure()
c4science · Help