diff --git a/tests/conftest.py b/tests/conftest.py index d0c840b..f1fb9f7 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,289 +1,309 @@ # -*- coding: utf-8 -*- # @file # LICENSE # # Copyright (©) 2016-2021 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 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): pass # try: # parser.addoption( # '--with-mpi', action="store_true", default=False, # help="Run MPI tests, this should be paired with mpirun." # ) # except ValueError: # pass 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, 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() +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="package", autouse=True) 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="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, domain, modes, size): - self.E = 3. - self.nu = 0. - self.e_star = self.E / (1 - self.nu**2) - self.n = size - self.domain = domain + 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, self.model.nu = self.E, self.nu + self.model.E = 3. + self.model.nu = 0. - self.pressure = profile(np.cos, size, modes, 1) - self.solution = profile(np.cos, size, modes, - 1 / (np.pi * self.e_star * norm(modes))) + 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=[tm.model_type.basic_2d], +@pytest.fixture(scope="package", + params=set(type_list) - {tm.model_type.volume_1d}, ids=mtidfn) def patch_westergaard(request): - return PatchWestergaard(request.param, [1., 1.], [3, 1], [6, 6]) + 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) def patch_isotropic_plasticity(request): return UniformPlasticity(request.param, [1., 1., 1.], [4, 4, 4]) diff --git a/tests/test_dumper.py b/tests/test_dumper.py index 3914060..c813d57 100644 --- a/tests/test_dumper.py +++ b/tests/test_dumper.py @@ -1,177 +1,172 @@ # -*- coding: utf-8 -*- # @file # LICENSE # # Copyright (©) 2016-2021 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 import os from pathlib import Path import tamaas as tm import numpy as np import pytest -from conftest import mtidfn +from conftest import mtidfn, type_list from tamaas.dumpers import NumpyDumper, JSONDumper, FieldDumper from tamaas.dumpers._helper import directory_dump, step_dump @pytest.fixture(scope="module", - params=[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], + params=type_list, ids=mtidfn) def model_fixture(request): dim = tm.type_traits[request.param].dimension model = tm.ModelFactory.createModel(request.param, dim * [1.], dim * [4]) model.displacement[:] = 3 model.traction[:] = 1 model['additional'] = 2 * np.ones(model.shape + [tm.type_traits[model.type].components]) return model class Dumper(tm.ModelDumper): """Simple numpy dumper""" def __init__(self): tm.ModelDumper.__init__(self) def dump(self, model): np.savetxt('tractions.txt', np.ravel(model.traction)) np.savetxt('displacement.txt', np.ravel(model.displacement)) dumper_types = [NumpyDumper, JSONDumper] try: from tamaas.dumpers import H5Dumper dumper_types.append(H5Dumper) except ImportError: pass try: from tamaas.dumpers import UVWDumper dumper_types.append(UVWDumper) except ImportError: pass try: from tamaas.dumpers import NetCDFDumper dumper_types.append(NetCDFDumper) except ImportError: pass @pytest.mark.parametrize('dumper_class', dumper_types) def test_dumper(model_fixture, tmp_path, dumper_class): os.chdir(tmp_path) filename = 'test_{}'.format(dumper_class.__name__) if not issubclass(dumper_class, FieldDumper): dumper = dumper_class(filename) else: dumper = dumper_class(filename, all_fields=True) filename = dumper.file_path dumper << model_fixture # UVW does not read VTK files if dumper_class is UVWDumper: with pytest.raises(NotImplementedError): dumper.read(filename) assert Path(filename).is_file() return rmodel = dumper.read(filename) assert rmodel.type == model_fixture.type assert rmodel.shape == model_fixture.shape for field in model_fixture: # assert field in rmodel assert np.all(np.abs(rmodel[field] - model_fixture[field]) < 1e-15) def test_custom_dumper(tmp_path): os.chdir(tmp_path) model = tm.ModelFactory.createModel(tm.model_type.volume_2d, [1., 1., 1.], [16, 4, 8]) dumper = Dumper() dumper << model tractions = np.loadtxt('tractions.txt') displacements = np.loadtxt('displacement.txt') assert np.all(tractions.reshape(model.traction.shape) == model.traction) assert np.all(displacements.reshape(model.displacement.shape) == model.displacement) @pytest.fixture def simple_model(scope="module"): return tm.ModelFactory.createModel(tm.model_type.basic_2d, [1, 1], [1, 1]) def test_directory_dump(simple_model, tmp_path): os.chdir(tmp_path) @directory_dump("dummy") class Dummy(tm.ModelDumper): def __init__(self): super(Dummy, self).__init__() def dump(self, model): pass @property def file_path(self): return 'dummy' dumper = Dummy() dumper << simple_model assert Path('dummy').is_dir() def test_step_dump(simple_model, tmp_path): os.chdir(tmp_path) @step_dump class Dummy(FieldDumper): extension = 'dummy' def dump_to_file(self, file_descriptor, model): with open(file_descriptor, 'w'): pass files = [] dumper = Dummy('dummy') def dump(): files.append(dumper.file_path) dumper << simple_model dump() dump() for file in files: assert Path(file).is_file() diff --git a/tests/test_patch_westergaard.py b/tests/test_patch_westergaard.py index e3fcdf6..149c075 100644 --- a/tests/test_patch_westergaard.py +++ b/tests/test_patch_westergaard.py @@ -1,43 +1,66 @@ # -*- coding: utf-8 -*- # @file # LICENSE # # Copyright (©) 2016-2021 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 from numpy.linalg import norm +import tamaas as tm def test_patch_westergaard(patch_westergaard): model = patch_westergaard.model solution = patch_westergaard.solution pressure = patch_westergaard.pressure - model['traction'][:] = pressure[:] - model.solveNeumann() - output_displ = model['displacement'] + + shape = (model.boundary_shape + [tm.type_traits[model.type].components]) + + # Setting the correct pressure component + traction = model.traction.reshape(shape)[..., -1] + traction[:] = pressure + + # solveNeumann() doesn't do anything on volume models + if model.type == tm.model_type.volume_2d: + tm.ModelFactory.registerVolumeOperators(model) + model.operators['boussinesq'](model.traction, model.displacement) + else: + model.solveNeumann() + + if model.type in [tm.model_type.volume_1d, tm.model_type.volume_2d]: + output_displ = model.displacement[0, ..., -1] + else: + output_displ = model.displacement.reshape(shape)[..., -1] error = norm(solution - output_displ) / norm(solution) assert error < 1e-15 and norm(solution) > 1e-15, \ "Neumann error = {}".format(error) - output_displ[:, :] = solution[:, :] + # Removing model types where no dirichlet is defined yet + if model.type in {tm.model_type.surface_1d, + tm.model_type.surface_2d, + tm.model_type.volume_2d}: + return + + traction[:] = 0 + output_displ[:] = solution model.solveDirichlet() - error = norm(pressure - model['traction']) / norm(pressure) + error = norm(pressure - traction) / norm(pressure) assert error < 1e-15, "Dirichlet error = {}".format(error)