diff --git a/tests/SConscript b/tests/SConscript index 42e6ca7..b835d98 100644 --- a/tests/SConscript +++ b/tests/SConscript @@ -1,112 +1,113 @@ from __future__ import print_function import os from os.path import join, abspath, basename # ------------------------------------------------------------------------------ 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_pressure.py test_westergaard.py test_patch_westergaard.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 fftfreq.py + conftest.py """) src_dir = "#/tests" build_dir = 'build-' + main_env['build_type'] + '/tests' for file in test_files: source = join(src_dir, file) test_env.Command(file, source, Copy("$TARGET", "$SOURCE")) # ------------------------------------------------------------------------------ def compile_google_test(env, gtest_path): gtest_obj = env.Object('gtest.o', [join(gtest_path, "src/gtest-all.cc")]) env.StaticLibrary('gtest', gtest_obj) # ------------------------------------------------------------------------------ def make_google_tests(env): gtest_dir = Dir('#/third-party/googletest/googletest') 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']) gtest_env.Tool(gtest) gtest_path = str(gtest_dir) 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 """) google_test_files += ['libgtest.a'] gtest_main = env.Object("gtest_main.o", 'tamaas_gtest_main.cc') env.Program('test_gtest_all', google_test_files + [gtest_main]) # ------------------------------------------------------------------------------ def make_bare_tests(env): env.Program("test_rough_surface.cpp") # ------------------------------------------------------------------------------ Import('main_env') # Setup of test environment test_env = main_env.Clone() test_env.AppendUnique(LIBS=['Tamaas'], LIBPATH=[abspath('build-' + main_env['build_type'] + '/src')]) 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.AppendUnique(LIBS=['python3.6m']) 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 new file mode 100644 index 0000000..17c0da0 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python +# coding: utf-8 +# ----------------------------------------------------------------------------- +# @author Lucas Frérot +# +# @section LICENSE +# +# Copyright (©) 2016 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 . +# ----------------------------------------------------------------------------- +from __future__ import division, print_function + +import numpy as np +import tamaas as tm +import pytest + + +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) + 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 = 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(): + tm.initialize(1) + yield HertzFixture(1024, 0.0001) + tm.finalize() diff --git a/tests/test_hertz_pressure.py b/tests/test_hertz_pressure.py index 08b48ef..66b552d 100644 --- a/tests/test_hertz_pressure.py +++ b/tests/test_hertz_pressure.py @@ -1,123 +1,70 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 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 . # ----------------------------------------------------------------------------- from __future__ import division, print_function import numpy as np import tamaas as tm +import pytest -def constructHertzProfile(size, curvature): - radius = 1. / curvature - x = np.linspace(-0.5, 0.5, size) - y = np.linspace(-0.5, 0.5, size) - x, y = np.meshgrid(x, y) - surface = np.sqrt(radius**2 - x**2 - y**2) - surface -= surface.mean() - return surface.copy() +@pytest.fixture(scope="module", + params=[tm.PolonskyKeerRey.pressure]) +def pkr(hertz, request): + model = tm.ModelFactory.createModel(tm.model_type.basic_2d, + [hertz.domain_size, hertz.domain_size], + [hertz.n, hertz.n]) + model.setElasticity(hertz.e_star, 0) + solver = tm.PolonskyKeerRey(model, hertz.surface, 1e-12, + request.param, + request.param) + solver.solve(hertz.load) + return model, hertz -def computeHertzDisplacement(e_star, contact_size, max_pressure, size): - x = np.linspace(-0.5, 0.5, size) - y = np.linspace(-0.5, 0.5, size) - x, y = np.meshgrid(x, y) - disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - - (x**2 + y**2)) - return disp.copy() - - -def test_hert_pressure(): - tm.initialize(1) - grid_size = 1024 - curvature = 0.1 - effective_modulus = 1. - load = 0.0001 - - surface = constructHertzProfile(grid_size, curvature) - model = tm.ModelFactory.createModel(tm.model_type_basic_2d, - [1., 1.], - [grid_size, grid_size]) - model.setElasticity(1, 0) - solver = tm.PolonskyKeerRey(model, surface, 1e-12, - tm.PolonskyKeerRey.pressure, - tm.PolonskyKeerRey.pressure) - solver.solve(load) - - tractions = model.getTraction() - displacements = model.getDisplacement() +def test_contact_area(pkr): # Testing contact area against Hertz solution for solids of revolution - contact_area = tm.SurfaceStatistics.computeContactArea(tractions) - hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.) + model, hertz = pkr + contact_area = np.sum(model.getTraction() > 0) / float(hertz.n**2) + hertz_contact_size = hertz.a hertz_area = np.pi * hertz_contact_size**2 area_error = np.abs(hertz_area - contact_area) / hertz_area - print("Area error: {}".format(area_error)) - - # Testing maximum pressure - max_pressure = tractions.max() - hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi - pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure - print("Max pressure error: {}".format(pressure_error)) - - # Testing displacements - hertz_displacements = computeHertzDisplacement(effective_modulus, - hertz_contact_size, - hertz_max_pressure, - grid_size) - - # Selecing only the points that are in contact - contact_indexes = [(i, j) for i in range(grid_size) for j in range(grid_size) - if tractions[i, j] > 0] - - # Displacements of bem are centered around the mean of the whole surface - # and Hertz displacements are not centered, so we need to compute mean - # on the contact zone for both arrays - bem_mean = 0. - hertz_mean = 0. - for index in contact_indexes: - bem_mean += displacements[index] - hertz_mean += hertz_displacements[index] - - bem_mean /= len(contact_indexes) - hertz_mean /= len(contact_indexes) - # Correction applied when computing error - correction = hertz_mean - bem_mean + assert area_error < 1e-3 - # Computation of error of displacement in contact zone - error = 0. - hertz_norm = 0. - for index in contact_indexes: - error += (hertz_displacements[index] - displacements[index] - correction)**2 - hertz_norm += (hertz_displacements[index] - hertz_mean)**2 - displacement_error = np.sqrt(error / hertz_norm) - print("Displacement error (in contact zone): {}".format(displacement_error)) +def test_pressure(pkr): + model, hertz = pkr + error = np.sum(model.getTraction() - hertz.pressure) / hertz.n**2 + assert error < 1e-6 - assert area_error < 1e-3 and pressure_error < 3e-3 and displacement_error < 1e-4 - tm.finalize() +def test_displacement(pkr): + model, hertz = pkr + error = np.sum(model.getDisplacement() - hertz.displacement) / hertz.n**2 + assert error < 2e-5 if __name__ == "__main__": - test_hert_pressure() + print("Test meant to be run with pytest")