diff --git a/tests/SConscript b/tests/SConscript index 215f610..2c4c761 100644 --- a/tests/SConscript +++ b/tests/SConscript @@ -1,163 +1,164 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @file # @section LICENSE # # Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 print_function from os.path import join, abspath from SCons.Script import Split, Copy, Dir, Import from detect import FindGTest # ------------------------------------------------------------------------------ def copyComStr(env, main): if 'SHCXXCOMSTR' in main: env['CXXCOMSTR'] = main['SHCXXCOMSTR'] if 'SHLINKCOMSTR' in main: env['LINKCOMSTR'] = main['SHLINKCOMSTR'] # ------------------------------------------------------------------------------ def make_python_tests(env): """Copy python tests to build directory""" test_env = env.Clone( PRINT_CMD_LINE_FUNC=main_env['gen_print']("Copying", "red", main_env)) test_files = Split(""" run_tests.sh test_hertz.py test_westergaard.py test_patch_westergaard.py + test_patch_plasticity.py test_surface.py test_autocorrelation.py test_hertz_disp.py test_hertz_kato.py test_saturated_pressure.py test_bem_grid.py test_flood_fill.py test_integral_operators.py test_dumper.py test_gtest.py test_tangential.py test_boussinesq_surface.py test_voigt.py test_memory.py fftfreq.py conftest.py """) src_dir = "#/tests" for file in test_files: source = join(src_dir, file) test_env.Command(file, source, Copy("$TARGET", "$SOURCE")) test_env = env.Clone() # Helper module for integral operators test_env['SHLIBPREFIX'] = '' register = test_env.SharedLibrary(target="register_integral_operators", source=["register_integral_operators.cpp"], LIBS=['Tamaas'], RPATH=[abspath('../src')]) Import('libTamaas') test_env.Depends(register, libTamaas) # ------------------------------------------------------------------------------ def compile_google_test(env, gtest_path): gtest_obj = env.Object('gtest.o', [join(gtest_path, "src/gtest-all.cc")]) return env.StaticLibrary('gtest', gtest_obj) # ------------------------------------------------------------------------------ def make_google_tests(env): gtest_dir = Dir(env['GTEST_ROOT']) gtest_env = env.Clone(CPPPATH=[gtest_dir], CXXFLAGS=['-pthread', '-isystem', join(str(gtest_dir), "include")]) colors = env['COLOR_DICT'] if not env['verbose']: gtest_env['ARCOMSTR'] = u'{}[Ar]{} $TARGET'.format(colors['purple'], colors['end']) gtest_env['RANLIBCOMSTR'] = \ u'{}[Randlib]{} $TARGET'.format(colors['purple'], colors['end']) FindGTest(gtest_env) libgtest = None # Hugly hack to detect if we need to compile gtest submodule if env['GTEST_ROOT'] == '#third-party/googletest/googletest': gtest_path = str(gtest_dir) libgtest = compile_google_test(gtest_env, gtest_path) env.AppendUnique(LIBS=['Tamaas'], CXXFLAGS=gtest_env['CXXFLAGS']) google_test_files = Split(""" test_fft.cpp test_grid.cpp test_loop.cpp test_model.cpp test_static_types.cpp test_integration.cpp """) gtest_main = env.Object("tamaas_gtest_main.o", 'tamaas_gtest_main.cc') gtest_all = env.Program('test_gtest_all', google_test_files + [gtest_main], LIBS=(env['LIBS'] + ['gtest'])) Import('libTamaas') env.Depends(gtest_all, libTamaas) env.Depends(gtest_all, libgtest) # ------------------------------------------------------------------------------ def make_bare_tests(env): rough = env.Program("test_rough_surface.cpp") Import('libTamaas') env.Depends(rough, libTamaas) # ------------------------------------------------------------------------------ Import('main_env') # Setup of test environment test_env = main_env.Clone() test_env.AppendUnique(LIBS=['Tamaas'], LIBPATH=['.']) copyComStr(test_env, main_env) # Building tests that do not require any third party make_bare_tests(test_env) if test_env['build_python']: test_env.Tool(pybind11) test_env.ParseConfig("{}-config --ldflags".format(test_env['py_exec'])) test_env['CCFLAGS'] = [] make_python_tests(test_env) # Building google tests if test_env['use_googletest']: make_google_tests(test_env) diff --git a/tests/conftest.py b/tests/conftest.py index a684644..3f715a0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,172 +1,215 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 numpy as np import tamaas as tm import pytest from numpy.linalg import norm def profile(func, size, mode, amplitude): x = np.linspace(0, 1, size[0], endpoint=False, dtype=tm.dtype) y = np.linspace(0, 1, size[1], endpoint=False, dtype=tm.dtype) x, y = np.meshgrid(x, y, indexing='ij') surface = amplitude * func(2*np.pi*x*mode[0]) * func(2*np.pi*y*mode[1]) return surface.copy() @pytest.fixture(scope="session") def tamaas_fixture(): tm.initialize() yield None tm.finalize() 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="module") def hertz(tamaas_fixture): return HertzFixture(1024, 0.00001) 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="module") def westergaard(tamaas_fixture): return WestergaardFixture(19683, 0.1) class PatchWestergaard: def __init__(self, model_type, domain, modes, size): self.E = 3. self.nu = 0. self.e_star = self.E / (1 - self.nu**2) self.n = size self.domain = domain self.model = tm.ModelFactory.createModel(model_type, domain, size) self.model.setElasticity(self.E, self.nu) self.pressure = profile(np.cos, size, modes, 1) self.solution = profile(np.cos, size, modes, 1 / (np.pi * self.e_star * norm(modes))) @pytest.fixture(scope="module", params=[tm.model_type.basic_2d]) def patch_westergaard(tamaas_fixture, request): return PatchWestergaard(request.param, [1., 1.], [3, 1], [6, 6]) + + +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 + } + + +@pytest.fixture(scope="module", params=[tm.model_type.volume_2d]) +def patch_isotropic_plasticity(tamaas_fixture, request): + return UniformPlasticity(request.param, [1., 1., 1.], [4, 4, 4]) diff --git a/tests/test_patch_plasticity.py b/tests/test_patch_plasticity.py new file mode 100644 index 0000000..424b951 --- /dev/null +++ b/tests/test_patch_plasticity.py @@ -0,0 +1,42 @@ +# -*- coding: utf-8 -*- +# @file +# @section LICENSE +# +# Copyright (©) 2016-20 EPFL (École Polytechnique Fédérale de Lausanne), +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# 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 numpy.linalg import norm +from tamaas.nonlinear_solvers import DFSANESolver + + +def test_patch_plasticity(patch_isotropic_plasticity): + model = patch_isotropic_plasticity.model + residual = patch_isotropic_plasticity.residual + + applied_pressure = 0.1 + + solver = DFSANESolver(residual, model) + pressure = model['traction'][..., 2] + pressure[:] = applied_pressure + + solver.solve() + solver.updateState() + + solution = patch_isotropic_plasticity.solution(applied_pressure) + + for key in solution: + error = norm(model[key] - solution[key]) / applied_pressure + assert error < 3e-15