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