diff --git a/tests/conftest.py b/tests/conftest.py index c99abea..d6f60ca 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,306 +1,309 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2022 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 . from __future__ import division, print_function import subprocess import sys from functools import reduce from operator import mul as multiply import numpy as np import pytest from numpy.linalg import norm try: import tamaas as tm except ModuleNotFoundError as e: print(e, file=sys.stderr) print("Use 'scons test' to run tests", file=sys.stderr) sys.exit(1) type_list = [ tm.model_type.basic_1d, tm.model_type.basic_2d, tm.model_type.surface_1d, tm.model_type.surface_2d, tm.model_type.volume_1d, tm.model_type.volume_2d ] def get_libtamaas_file(): try: ldd_out = bytes(subprocess.check_output(['ldd', tm._tamaas.__file__])) for line in filter(lambda x: 'Tamaas' in x, ldd_out.decode().split('\n')): path = line.split(' => ')[1] return path.split(' ')[0] return "static link" except subprocess.CalledProcessError: return "N/A" def pytest_report_header(config): print('tamaas build type: {}\nmodule file: {}\nwrapper: {}\nlibTamaas: {}'. format(tm.TamaasInfo.build_type, tm.__file__, tm._tamaas.__file__, get_libtamaas_file())) def pytest_addoption(parser): parser.addoption('--log-debug', action='store_true', default=False) def mtidfn(val): return str(val).split('.')[-1] def pkridfn(val): return { tm.PolonskyKeerRey.pressure: "pkr_pressure", tm.PolonskyKeerRey.gap: "pkr_gap" }[val] def profile(func, domain, N, modes, amplitude): coords = [ np.linspace(0, d, n, endpoint=False, dtype=tm.dtype) for n, d in zip(N, domain) ] coords = np.meshgrid(*coords, indexing='ij') sines = [func(2 * np.pi * x * m) for x, m in zip(coords, modes)] return reduce(multiply, [amplitude] + sines) @pytest.fixture(scope='session', autouse=True) def log_level(pytestconfig): if pytestconfig.option.log_debug: tm.set_log_level(tm.LogLevel.debug) class HertzFixture: def __init__(self, n, load): self.domain_size = 1 self.n = n self.load = load self.curvature = 0.1 self.radius = 1. / self.curvature self.e_star = 1. self.a = (3 * load / (4 * self.curvature * self.e_star))**(1. / 3.) self.x = np.linspace(-self.domain_size / 2., self.domain_size / 2., self.n, dtype=tm.dtype) self.y = self.x.copy() self.x, self.y = np.meshgrid(self.x, self.y) self._computeSurface() self._computePressure() self._computeDisplacement() def _computeDisplacement(self): r = np.sqrt(self.x**2 + self.y**2) self.displacement = np.zeros_like(r) contact = r < self.a self.displacement[contact] = self.surface[contact] self.displacement[~contact] = \ (self.surface[~contact] + self.a / (np.pi * self.radius) * np.sqrt(r[~contact]**2 - self.a**2) + (r[~contact]**2 - 2 * self.a**2) / (np.pi * self.radius) * np.arccos(self.a / r[~contact])) def _computePressure(self): r = np.sqrt(self.x**2 + self.y**2) self.pressure = np.zeros_like(r) contact = np.where(r < self.a) self.pressure[contact] = \ 2 * self.e_star / (np.pi * self.radius) \ * np.sqrt(self.a**2 - r[contact]**2) def _computeSurface(self): self.surface = -1. / (2 * self.radius) * (self.x**2 + self.y**2) @pytest.fixture(scope="package") def hertz(): return HertzFixture(1024, 0.00001) @pytest.fixture(scope="package") def hertz_coarse(): return HertzFixture(512, 0.0001) @pytest.fixture(scope="package", params=[tm.PolonskyKeerRey.pressure], ids=pkridfn) def pkr(hertz, request): model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [hertz.domain_size, hertz.domain_size], [hertz.n, hertz.n]) model.E, model.nu = hertz.e_star, 0 solver = tm.PolonskyKeerRey(model, hertz.surface, 1e-12, request.param, request.param) solver.solve(hertz.load) return model, hertz class WestergaardFixture: def __init__(self, n, load): self.domain_size = 1. self.lamda = 1. self.delta = 0.1 self.e_star = 1. self.n = n self.p_star = np.pi * self.e_star * self.delta / self.lamda self.load = load * self.p_star self.a = self.lamda / np.pi \ * np.arcsin(np.sqrt(self.load / self.p_star)) self.x = np.linspace(-self.domain_size / 2., self.domain_size / 2., self.n, endpoint=False, dtype=tm.dtype) self._computeSurface() self._computePressure() self._computeDisplacement() def _computeSurface(self): self.surface = self.delta * np.cos(2 * np.pi * self.x / self.lamda) def _computePressure(self): self.pressure = np.zeros_like(self.surface) contact = np.where(np.abs(self.x) < self.a) self.pressure[contact] = 2 * self.load \ * (np.cos(np.pi * self.x[contact] / self.lamda) / np.sin(np.pi * self.a / self.lamda)**2) \ * np.sqrt(np.sin(np.pi * self.a / self.lamda)**2 - np.sin(np.pi * self.x[contact] / self.lamda)**2) def _computeDisplacement(self): psi = np.pi * np.abs(self.x) / self.lamda psi_a = np.pi * self.a / self.lamda with np.errstate(invalid='ignore'): # get some warnings out of the way self.displacement = ( np.cos(2 * psi) + 2 * np.sin(psi) * np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2) - 2 * np.sin(psi_a)**2 * np.log( (np.sin(psi) + np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2)) / np.sin(psi_a))) contact = np.where(np.abs(self.x) < self.a) self.displacement[contact] = np.cos(2 * psi[contact]) self.displacement *= self.load * self.lamda / (np.pi * self.e_star * np.sin(psi_a)**2) @pytest.fixture(scope="package") def westergaard(): return WestergaardFixture(19683, 0.1) class PatchWestergaard: def __init__(self, model_type): dim = tm.type_traits[model_type].dimension bdim = tm.type_traits[model_type].boundary_dimension domain = [1.] * dim size = [6] * dim self.modes = np.random.randint(1, 3, (bdim, )) self.model = tm.ModelFactory.createModel(model_type, domain, size) self.model.E = 3. self.model.nu = 0. self.pressure = profile(np.cos, self.model.boundary_system_size, self.model.boundary_shape, self.modes, 1) self.solution = profile( np.cos, self.model.boundary_system_size, self.model.boundary_shape, self.modes, 1 / (np.pi * self.model.E_star * norm(self.modes))) @pytest.fixture(scope="package", params=set(type_list) - {tm.model_type.volume_1d}, ids=mtidfn) def patch_westergaard(request): return PatchWestergaard(request.param) class UniformPlasticity: def __init__(self, model_type, domain, sizes): self.n = sizes self.domain = domain self.model = tm.ModelFactory.createModel(model_type, domain, sizes) self.E_h = 0.1 self.sigma_y = 0.01 self.residual = tm.ModelFactory.createResidual(self.model, sigma_y=self.sigma_y, hardening=self.E_h) self.model.E = 1. self.model.nu = 0. def solution(self, p): E, nu = self.model.E, self.model.nu E_h, sigma_y = self.E_h, self.sigma_y mu = E / (2 * (1 + nu)) strain = -1 / (mu + E_h) * (p * (3 * mu + E_h) / (2 * mu) - np.sign(p) * sigma_y) dep = (2 * mu * np.abs(strain) - sigma_y) / (3 * mu + E_h) plastic_strain = np.sign(strain) / 2 * dep * np.array([ -1, -1, 2, 0, 0, 0, ], dtype=tm.dtype) stress = 2 * mu * (np.array([0, 0, strain, 0, 0, 0], dtype=tm.dtype) - plastic_strain) return { "stress": stress, "plastic_strain": plastic_strain, "cumulated_plastic_strain": dep }, { "stress": mu, "plastic_strain": 1, "cumulated_plastic_strain": 1 } -@pytest.fixture(params=[tm.model_type.volume_2d], ids=mtidfn) +# Tests plane-strain situations with Ny = 1 +@pytest.fixture(params=[1, 4], ids=mtidfn) def patch_isotropic_plasticity(request): - return UniformPlasticity(request.param, [1., 1., 1.], [4, 4, 4]) + return UniformPlasticity(tm.model_type.volume_2d, + [1., 1., 1.], + [4, 4, request.param]) diff --git a/tests/test_patch_plasticity.py b/tests/test_patch_plasticity.py index 1d93c54..85fef2d 100644 --- a/tests/test_patch_plasticity.py +++ b/tests/test_patch_plasticity.py @@ -1,52 +1,51 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-20 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2022 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 . import pytest -from numpy.linalg import norm +from numpy.testing import assert_allclose import tamaas from tamaas.nonlinear_solvers import \ DFSANESolver as PySolver, \ DFSANECXXSolver as CXXSolver, \ NewtonKrylovSolver @pytest.mark.parametrize('solvers', [CXXSolver, NewtonKrylovSolver, PySolver]) def test_patch_plasticity(patch_isotropic_plasticity, solvers): "Test analyitical solution of 1d plasticity" tamaas.set_log_level(tamaas.LogLevel.info) Solver = solvers model = patch_isotropic_plasticity.model residual = patch_isotropic_plasticity.residual applied_pressure = 0.1 solver = Solver(residual) solver.tolerance = 1e-15 pressure = model['traction'][..., 2] pressure[:] = applied_pressure solver.solve() solver.updateState() solution, normal = patch_isotropic_plasticity.solution(applied_pressure) for key in solution: - error = norm(model[key] - solution[key]) / normal[key] - assert error < 3e-15 + assert_allclose((model[key] - solution[key]) / normal[key], 0, atol=3e-15)