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)