diff --git a/.arclint b/.arclint new file mode 100644 index 0000000..73c2265 --- /dev/null +++ b/.arclint @@ -0,0 +1,18 @@ +{ + "linters": { + "python": { + "type": "flake8", + "include": [ "(\\.py$)", "(^hbm$)" ], + "exclude": ["(^third-party/)", "(legacy/)"], + "severity.rules": { + "(^E)": "error", + "(^W)": "warning", + "(^F)": "advice" + } + }, + + "merges": { + "type": "merge-conflict" + } + } +} diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..15fc7e3 --- /dev/null +++ b/.flake8 @@ -0,0 +1,2 @@ +[flake8] +max-line-length = 80 diff --git a/doc/Makefile b/doc/Makefile index 0c5906e..ed502aa 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -1,32 +1,32 @@ DOXYGEN = doxygen SPHINX = sphinx-build SPHINX_SOURCE_DIR = sphinx/source SPHINX_SOURCES := $(wildcard $(SPHINX_SOURCE_DIR)/*.rst) SPHINX_SOURCES += $(SPHINX_SOURCE_DIR)/conf.py BUILD_DIR = build DOXYGEN_BUILD_DIR = $(BUILD_DIR)/doxygen SPHINX_BUILD_DIR = $(BUILD_DIR)/sphinx/html # Default target -all:sphinx +all:doxygen sphinx # Phony targets for aliases .PHONY:doxygen sphinx clean # Target aliases doxygen:$(DOXYGEN_BUILD_DIR) sphinx:$(SPHINX_BUILD_DIR) # Build doxygen documentation $(DOXYGEN_BUILD_DIR):doxygen/Doxyfile @mkdir -p $@ $(DOXYGEN) $< $(SPHINX_BUILD_DIR):$(SPHINX_SOURCES) $(DOXYGEN_BUILD_DIR) @mkdir -p $@ $(SPHINX) -b html $(SPHINX_SOURCE_DIR) $@ clean: $(RM) -r $(BUILD_DIR) diff --git a/examples/SConstruct b/examples/SConstruct deleted file mode 100644 index af60e3b..0000000 --- a/examples/SConstruct +++ /dev/null @@ -1,31 +0,0 @@ -# -*- 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 . - -import os - -env = Environment(CC = 'g++', - CXXFLAGS = ['-std=c++11', '-g', '-fopenmp', '-O3'], - TOOLS = ['default', 'tamaas'], - TAMAAS_PATH = '../build-release/', - ENV = os.environ) - -env.Program('rough_surface.cc') diff --git a/examples/rough_contact.py b/examples/rough_contact.py index 10f145b..abfb80e 100644 --- a/examples/rough_contact.py +++ b/examples/rough_contact.py @@ -1,74 +1,74 @@ #!/usr/bin/env python3 # @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 . import tamaas as tm import matplotlib.pyplot as plt # Initialize threads and fftw tm.initialize() tm.set_log_level(tm.LogLevel.info) # Show progression of solver # Surface size n = 512 # Surface generator sg = tm.SurfaceGeneratorFilter2D() sg.setSizes([n, n]) sg.random_seed = 0 # Spectrum spectrum = tm.Isopowerlaw2D() sg.setFilter(spectrum) # Parameters spectrum.q0 = 16 spectrum.q1 = 16 spectrum.q2 = 64 spectrum.hurst = 0.8 # Generating surface surface = sg.buildSurface() -#surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface) +# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface) surface /= n -#print(spectrum.rmsSlopes()) -#print(tm.SurfaceStatistics.computeRMSSlope(surface)) +# print(spectrum.rmsSlopes()) +# print(tm.SurfaceStatistics.computeRMSSlope(surface)) plt.imshow(surface) # Creating model model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n]) # Solver solver = tm.PolonskyKeerRey(model, surface, 1e-12, tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.pressure) # Solve for target pressure p_target = 0.1 solver.solve(p_target) plt.figure() plt.imshow(model.getTraction()) print(model.getTraction().mean()) # Cleanup threads tm.finalize() plt.show() diff --git a/examples/site_scons/site_init.py b/examples/site_scons/site_init.py deleted file mode 100644 index a6745c5..0000000 --- a/examples/site_scons/site_init.py +++ /dev/null @@ -1,151 +0,0 @@ -# @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 SCons.Script import * -from os.path import * - -import os - - -def fftw(env): - """A Tool to search for fftw headers and libraries""" - - if env.GetOption('clean'): - return - env.SetDefault(FFTW_VERSION='3') - env.SetDefault(FFTW_LIBRARY_WISH=[]) - - print("Building for fftw version {}".format(env['FFTW_VERSION'])) - - if 'FFTW_INCLUDE_DIR' in env: - env.AppendUnique(CPPPATH=env['FFTW_INCLUDE_DIR']) - - if 'FFTW_LIBRARY_DIR' in env: - env.AppendUnique(LIBPATH=env['FFTW_LIBRARY_DIR']) - - version = env['FFTW_VERSION'] - if version == "2": - lib_names = {'main': 'fftw'} - inc_names = ['fftw.h'] - else: - lib_names = {'main': 'fftw3', - 'thread': 'fftw3_threads', - 'omp': 'fftw3_omp'} - inc_names = ['fftw3.h'] - - try: - libs = [lib_names[i] for i in env['FFTW_LIBRARY_WISH']] - except: - raise SCons.Errors.StopError( - 'Incompatible wishlist {0} from version {1}'.format( - env['FFTW_LIBRARY_WISH'], env['FFTW_VERSION'])) - - env.AppendUnique(LIBS=libs) - if version == "2": - env.Append(LIBS='m') - - conf = Configure(env) - if not conf.CheckLibWithHeader(libs, inc_names, 'c++'): - raise SCons.Errors.StopError( - 'Failed to find libraries {0} or ' - 'headers {1}.'.format(str(lib_names), str(inc_names))) - - env = conf.Finish() - - -def criterion(env): - """A tool to configure for Criterion (test suite)""" - if env.GetOption('clean'): - return - - env.AppendUnique(CPPPATH=[Dir("#third-party/Criterion/include")], - LIBPATH=[abspath(str(Dir("#third-party/Criterion/build")))]) - conf = Configure(env) - # import pdb - # pdb.set_trace() - if not conf.CheckLibWithHeader('criterion', - 'criterion/criterion.h', - 'c++'): - raise SCons.Errors.StopError( - 'Failed to find library Criterion or criterion.h header in third-party') - env = conf.Finish() - - -def boost(env): - """A tool to confugure for boost""" - if env.GetOption('clean'): - return - - if 'BOOST_INCLUDE_DIR' in env: - env.AppendUnique(CPPPATH=env['BOOST_INCLUDE_DIR']) - conf = Configure(env) - if not conf.CheckCXXHeader('boost/version.hpp'): - raise SCons.Errors.StopError( - 'Failed to find Boost library') - env = conf.Finish() - - -def thrust(env): - """A tool to confugure for thrust""" - if env.GetOption('clean'): - return - - if 'THRUST_INCLUDE_DIR' in env: - env.AppendUnique(CPPPATH=env['THRUST_INCLUDE_DIR']) - conf = Configure(env) - if not conf.CheckCXXHeader('thrust/version.h'): - raise SCons.Errors.StopError( - 'Failed to find Thrust library') - env = conf.Finish() - - -def tamaas(env): - """Sets correct environement variables for Tamaas. - - 'TAMAAS_PATH' is list of potential tamaas repositories - - 'TAMAAS_BUILD_TYPE' is the build type of tamaas (release, debug, profiling)""" - if env.GetOption("clean"): - return - - try: - build_type = env['TAMAAS_BUILD_TYPE'] - except: - build_type = "release" - - if 'TAMAAS_PATH' not in env: - print("Please set the 'TAMAAS_PATH' variable in your scons environment") - Exit(1) - - tm_path = env['TAMAAS_PATH'] - include_dir = tm_path + '/src' - subdirs = "core model solvers gpu bem surface".split(" ") - include_dirs = [abspath(join(include_dir, x)) for x in subdirs] - lib_dir = join(tm_path, 'build-{}/src'.format(build_type)) - - env.AppendUnique(CPPPATH = include_dirs) - env.AppendUnique(CPPPATH = ['/opt/cuda/include']) # TODO fix this dirty hack - env.AppendUnique(LIBPATH = [abspath(lib_dir)]) - env.AppendUnique(RPATH = [abspath(lib_dir)]) - env.AppendUnique(CXXFLAGS = ["-std=c++14"]) - - conf = Configure(env) - if not conf.CheckLibWithHeader("Tamaas", "tamaas.hh", language="CXX"): - raise SCons.Errors.StopError( - "Tamaas ({}) not found in {}".format(build_type, abspath(tm_path))) - env['TAMAAS_PATH'] = abspath(tm_path) - env = conf.Finish() diff --git a/examples/statistics.py b/examples/statistics.py index 1168807..528c00a 100644 --- a/examples/statistics.py +++ b/examples/statistics.py @@ -1,86 +1,87 @@ #!/usr/bin/env python3 # @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 . import os import numpy as np import matplotlib.pyplot as plt import tamaas as tm from matplotlib.colors import LogNorm # Initialize threads and fftw tm.initialize() tm.set_log_level(tm.LogLevel.info) # Show progression of solver # Surface size n = 512 # Surface generator sg = tm.SurfaceGeneratorFilter2D() sg.setSizes([n, n]) sg.random_seed = 0 # Spectrum spectrum = tm.Isopowerlaw2D() sg.setFilter(spectrum) # Parameters spectrum.q0 = 16 spectrum.q1 = 16 spectrum.q2 = 64 spectrum.hurst = 0.8 # Generating surface surface = sg.buildSurface() # Computing PSD and shifting for plot psd = tm.Statistics2D.computePowerSpectrum(surface) psd = np.fft.fftshift(psd, axes=0) plt.imshow(psd.real, norm=LogNorm()) plt.gca().set_title('Power Spectrum Density') plt.gcf().tight_layout() # Computing autocorrelation and shifting for plot acf = tm.Statistics2D.computeAutocorrelation(surface) acf = np.fft.fftshift(acf) plt.figure() plt.imshow(acf) plt.gca().set_title('Autocorrelation') plt.gcf().tight_layout() plt.show() # Write the rough surface in paraview's VTK format try: import uvw try: os.mkdir('paraview') - except FileExistsError as e: + except FileExistsError: pass x = np.linspace(0, 1, n, endpoint=True) y = x.copy() - with uvw.RectilinearGrid(os.path.join('paraview', 'surface.vtr'), (x, y)) as grid: + with uvw.RectilinearGrid(os.path.join('paraview', 'surface.vtr'), + (x, y)) as grid: grid.addPointData(uvw.DataArray(surface, range(2), 'surface')) -except ImportError as e: +except ImportError: print("uvw not installed, will not write VTK file") pass diff --git a/python/tamaas/__init__.py b/python/tamaas/__init__.py index 00e277c..7818ba5 100644 --- a/python/tamaas/__init__.py +++ b/python/tamaas/__init__.py @@ -1,71 +1,73 @@ # @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 . """ Tamaas is a library dedicated to the fast treatment of rough contact problems. See README.md and examples for guidance as how to use Tamaas. See __authors__, __license__, __copyright__ for extra information about Tamaas. """ __authors__ = [ "Guillaume Anciaux ", "Lucas Frérot ", "Son Pham-Ba ", "Valentine Rey ", ] __copyright__ = "Copyright (©) 2016-19 EPFL " \ + "(École Polytechnique Fédérale de Lausanne), " \ + "Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" __license__ = "AGPLv3" __maintainer__ = "Lucas Frérot" __email__ = "lucas.frerot@epfl.ch" -try: - from ._tamaas import * # noqa +try: + from . import _tamaas as _tm # Redefinition of model_type constants (for compatibility) - model_type_basic_1d = model_type.basic_1d - model_type_basic_2d = model_type.basic_2d - model_type_surface_1d = model_type.surface_1d - model_type_surface_2d = model_type.surface_2d - model_type_volume_1d = model_type.volume_1d - model_type_volume_2d = model_type.volume_2d + model_type_basic_1d = _tm.model_type.basic_1d + model_type_basic_2d = _tm.model_type.basic_2d + model_type_surface_1d = _tm.model_type.surface_1d + model_type_surface_2d = _tm.model_type.surface_2d + model_type_volume_1d = _tm.model_type.volume_1d + model_type_volume_2d = _tm.model_type.volume_2d + + from ._tamaas import * # noqa except ImportError as e: print("Error trying to import _tamaas:\n{}".format(e)) raise e class ParamHelper: """Legacy class to manage parameters/setters/getters""" def __init__(self, obj): self.obj = obj def set(self, params): for key in params: setter_name = 'set' + key try: accessor = getattr(self.obj, setter_name) accessor(params[key]) except Exception: print("Setter '{}({})' does not exist for object {}" .format(setter_name, type(params[key]), self.obj)) diff --git a/site_scons/site_tools/nvcc.py b/site_scons/site_tools/nvcc.py index a385268..7751922 100644 --- a/site_scons/site_tools/nvcc.py +++ b/site_scons/site_tools/nvcc.py @@ -1,205 +1,223 @@ # -*- 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 . """ SCons tool to setup compilation environment for NVidia's NVCC Based on version https://github.com/thrust/thrust/blob/master/site_scons/site_tools/nvcc.py """ from __future__ import print_function import SCons from SCons.Errors import StopError import platform from os.path import join known_components = { 'cufft': 'cufft.h', 'cufftw': 'cufftw.h' } def _find_cuda_paths(env): """Finding cuda paths""" if 'CUDA_TOOLKIT_PATH' not in env: raise StopError("CUDA_TOOLKIT_PATH variable not found in environment") cuda_path = env['CUDA_TOOLKIT_PATH'] inc_path = join(cuda_path, 'include') lib_path = join(cuda_path, 'lib') bin_path = join(cuda_path, 'bin') if platform.machine()[-2:] == '64': lib_path += '64' # add nvcc's location to PATH env.PrependENVPath('PATH', bin_path) # add the location of the cudart shared library to LD_LIBRARY_PATH as well # this allows us to execute CUDA programs during the build env.PrependENVPath('LD_LIBRARY_PATH', lib_path) # add to include paths env.AppendUnique(CPPPATH=[inc_path]) # add to library paths env.AppendUnique(LIBPATH=[lib_path]) def _configure_cuda_components(env): """Finding required components in cuda""" if env.GetOption('clean'): return if 'CUDA_COMPONENTS' in env: for component in env['CUDA_COMPONENTS']: if component not in known_components: raise StopError("Unknown cuda component '{}'".format(component)) conf = SCons.Script.Configure(env) - if not conf.CheckLibWithHeader(component, known_components[component], 'c++'): + if not conf.CheckLibWithHeader(component, + known_components[component], + 'c++'): raise StopError("Failed to find library {} or header {}".format( component, known_components[component])) env = conf.Finish() def _define_compile_cuda(env): """Define the compile string for files that need nvcc""" env['NVCC'] = env.Detect('nvcc') # Setting compilation command env['NVCCCOM'] = \ - '$NVCC $CUDA_ARCH_FLAG -o $TARGET -c $NVCC_CXXFLAGS $NVCC_CCFLAGS $_CCCOMCOM -x cu $SOURCES' + '$NVCC $CUDA_ARCH_FLAG -o $TARGET -c $NVCC_CXXFLAGS $NVCC_CCFLAGS ' \ + + '$_CCCOMCOM -x cu $SOURCES' env['SHNVCCCOM'] = \ - '$NVCC $CUDA_ARCH_FLAG -o $TARGET -dc -Xcompiler -fPIC $NVCC_CXXFLAGS $NVCC_CCFLAGS $_CCCOMCOM -x cu $SOURCES' + '$NVCC $CUDA_ARCH_FLAG -o $TARGET -dc -Xcompiler -fPIC ' \ + + '$NVCC_CXXFLAGS $NVCC_CCFLAGS $_CCCOMCOM -x cu $SOURCES' # Constructing proper cc flags env['_NVCC_BARE_CCFLAGS'] = \ '${_concat("", CCFLAGS, "", __env__, _NVCC_BARE_FILTER)}' env['_NVCC_PASSED_CCFLAGS'] = \ - '${_concat("-Xcompiler ", CCFLAGS, "", __env__, _NVCC_COMPILER_PASSED_FILTER)}' + '${_concat("-Xcompiler ", CCFLAGS, "", __env__, ' \ + + '_NVCC_COMPILER_PASSED_FILTER)}' # Constructing proper cxx flags env['_NVCC_BARE_CXXFLAGS'] = \ '${_concat("", CXXFLAGS, "", __env__, _NVCC_BARE_FILTER)}' env['_NVCC_PASSED_CXXFLAGS'] = \ - '${_concat("-Xcompiler ", CXXFLAGS, "", __env__, _NVCC_COMPILER_PASSED_FILTER)}' + '${_concat("-Xcompiler ", CXXFLAGS, "", __env__, ' \ + + '_NVCC_COMPILER_PASSED_FILTER)}' # Putting all together env['NVCC_CCFLAGS'] = '$_NVCC_BARE_CCFLAGS $_NVCC_PASSED_CCFLAGS' env['NVCC_CXXFLAGS'] = '$_NVCC_BARE_CXXFLAGS $_NVCC_PASSED_CXXFLAGS' def _define_link_cuda(env): """Define the link string""" # Fixing rpaths env['RPATHPREFIX'] = '-rpath=' env['__RPATH'] = '${_concat("-Xlinker ", _RPATH, "", __env__)}' env['LINK'] = env['NVCC'] env['SHLINK'] = env['NVCC'] # Replacing old link command strings env['LINKCOM'] = \ - '$LINK $CUDA_ARCH_FLAG -o $TARGET $NVCC_LINKFLAGS $__RPATH $SOURCES $_LIBDIRFLAGS $_LIBFLAGS' + '$LINK $CUDA_ARCH_FLAG -o $TARGET $NVCC_LINKFLAGS $__RPATH ' \ + + '$SOURCES $_LIBDIRFLAGS $_LIBFLAGS' env['SHLINKCOM'] = \ - '$SHLINK -shared $CUDA_ARCH_FLAG -o $TARGET $NVCC_SHLINKFLAGS $__SHLIBVERSIONFLAGS $__RPATH $SOURCES $_LIBDIRFLAGS $_LIBFLAGS' + '$SHLINK -shared $CUDA_ARCH_FLAG -o $TARGET $NVCC_SHLINKFLAGS ' \ + + '$__SHLIBVERSIONFLAGS $__RPATH $SOURCES $_LIBDIRFLAGS $_LIBFLAGS' # Constructing proper static linker flags env['_NVCC_BARE_LINKFLAGS'] = \ '${_concat("", LINKFLAGS, "", __env__, _NVCC_BARE_FILTER)}' env['_NVCC_COMPILER_LINKFLAGS'] = \ - '${_concat("-Xcompiler ", LINKFLAGS, "", __env__, _NVCC_COMPILER_FILTER)}' + '${_concat("-Xcompiler ", LINKFLAGS, "", __env__, ' \ + + '_NVCC_COMPILER_FILTER)}' env['_NVCC_PASSED_LINKFLAGS'] = \ - '${_concat("-Xlinker ", LINKFLAGS, "", __env__, _NVCC_LINK_PASSED_FILTER)}' + '${_concat("-Xlinker ", LINKFLAGS, "", __env__, ' \ + + '_NVCC_LINK_PASSED_FILTER)}' # Constructing proper shared linker flags env['_NVCC_BARE_SHLINKFLAGS'] = \ '${_concat("", SHLINKFLAGS, "", __env__, _NVCC_BARE_FILTER)}' env['_NVCC_COMPILER_SHLINKFLAGS'] = \ - '${_concat("-Xcompiler ", LINKFLAGS, "", __env__, _NVCC_COMPILER_FILTER)}' + '${_concat("-Xcompiler ", LINKFLAGS, "", __env__, ' \ + + '_NVCC_COMPILER_FILTER)}' env['_NVCC_PASSED_SHLINKFLAGS'] = \ - '${_concat("-Xlinker ", SHLINKFLAGS, "", __env__, _NVCC_LINK_PASSED_FILTER)}' + '${_concat("-Xlinker ", SHLINKFLAGS, "", __env__, ' \ + + '_NVCC_LINK_PASSED_FILTER)}' # Putting all together env['NVCC_LINKFLAGS'] = \ - '$_NVCC_BARE_LINKFLAGS $_NVCC_COMPILER_LINKFLAGS $_NVCC_PASSED_LINKFLAGS' + '$_NVCC_BARE_LINKFLAGS $_NVCC_COMPILER_LINKFLAGS ' \ + + '$_NVCC_PASSED_LINKFLAGS' env['NVCC_SHLINKFLAGS'] = env['NVCC_LINKFLAGS'] env['NVCC_SHLINKFLAGS'] += \ - ' $_NVCC_BARE_SHLINKFLAGS $_NVCC_COMPILER_SHLINKFLAGS $_NVCC_PASSED_SHLINKFLAGS' + ' $_NVCC_BARE_SHLINKFLAGS $_NVCC_COMPILER_SHLINKFLAGS ' \ + + '$_NVCC_PASSED_SHLINKFLAGS' def _define_commands(env): """Defining the command strings""" # Flags allowed by nvcc bare_flags = """-std=c++11 -O0 -O1 -O2 -O3 -g -pg -G -w""".split() # Experimental flags bare_flags += """-expt-extended-lambda -expt-relaxed-constexpr""".split() # Flags that must be passed to compiler at link time compiler_flags = """-fopenmp -qopenmp""".split() # Pass flags bare to nvcc - env['_NVCC_BARE_FILTER'] = lambda flags: list(filter(lambda flag: flag in bare_flags, flags)) + env['_NVCC_BARE_FILTER'] = lambda flags: list(filter( + lambda flag: flag in bare_flags, flags)) # Must prepend -Xcompiler, even at link time env['_NVCC_COMPILER_FILTER'] = lambda flags: list(filter( lambda flag: flag in compiler_flags, flags)) # Prepend -Xcompiler env['_NVCC_COMPILER_PASSED_FILTER'] = lambda flags: list(filter( lambda flag: flag not in set(bare_flags), flags)) # Prepend -Xlinker env['_NVCC_LINKED_PASSED_FILTER'] = lambda flags: list(filter( lambda flag: flag not in set(bare_flags) | set(compiler_flags), flags)) _define_compile_cuda(env) _define_link_cuda(env) def _add_actions_cuda(env): - """Adding actions to .cu files to compile with nvcc. Other files are not affected.""" + """ + Adding actions to .cu files to compile with nvcc. + Other files are not affected. + """ nvcc_action = SCons.Action.Action('$NVCCCOM', '$NVCCCOMSTR') shnvcc_action = SCons.Action.Action('$SHNVCCCOM', '$NVCCCOMSTR') static, shared = SCons.Tool.createObjBuilders(env) # Compiling with nvcc action added to detected .cu files static.add_action('.cu', nvcc_action) shared.add_action('.cu', shnvcc_action) # Emitter to qualify correctly object code type static.add_emitter('.cu', SCons.Defaults.StaticObjectEmitter) shared.add_emitter('.cu', SCons.Defaults.SharedObjectEmitter) env['CXXCOM'] = '$NVCCCOM' env['SHCXXCOM'] = '$SHNVCCCOM' # Scanner for dependency calculations SCons.Tool.SourceFileScanner.add_scanner('.cu', SCons.Scanner.C.CScanner()) def generate(env): """Setup environment for compiling .cu files with nvcc""" _find_cuda_paths(env) _add_actions_cuda(env) _define_commands(env) _configure_cuda_components(env) def exists(env): return env.Detect('nvcc') diff --git a/site_scons/site_tools/tamaas.py b/site_scons/site_tools/tamaas.py index 6c6065e..d018941 100644 --- a/site_scons/site_tools/tamaas.py +++ b/site_scons/site_tools/tamaas.py @@ -1,61 +1,63 @@ # -*- 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 . import os from SCons.Errors import StopError from SCons.Script import Configure def generate(env): """Sets correct environement variables for Tamaas. - 'TAMAAS_PATH' is the path to the build directory of Tamaas""" if env.GetOption("clean"): return if 'TAMAAS_PATH' not in env: raise StopError("Please set the 'TAMAAS_PATH' variable in your " + "scons environment") tm_path = env['TAMAAS_PATH'] include_dir = tm_path + '/src' subdirs = "core model solvers bem surface percolation".split(" ") - include_dirs = [os.path.abspath(os.path.join(include_dir, x)) for x in subdirs] + include_dirs = [os.path.abspath(os.path.join(include_dir, x)) + for x in subdirs] lib_dir = include_dir env.AppendUnique(CPPPATH=include_dirs) env.AppendUnique(LIBPATH=[os.path.abspath(lib_dir)]) env.AppendUnique(RPATH=[os.path.abspath(lib_dir)]) env.AppendUnique(CXXFLAGS=["-std=c++14"]) - env.AppendUnique(CPPDEFINES=["THRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP"]) + env.AppendUnique( + CPPDEFINES=["THRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP"]) env.AppendUnique(LIBS=['fftw3', 'fftw3_omp']) conf = Configure(env) if not conf.CheckLibWithHeader("Tamaas", "tamaas.hh", language="CXX"): raise StopError( "Tamaas not found in {}".format(os.path.abspath(tm_path))) env['TAMAAS_PATH'] = os.path.abspath(tm_path) env = conf.Finish() def exists(env): if 'TAMAAS_PATH' in env: return True return False diff --git a/site_scons/version.py b/site_scons/version.py index a5704db..d8e4a07 100644 --- a/site_scons/version.py +++ b/site_scons/version.py @@ -1,61 +1,62 @@ # -*- 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 . import subprocess import base64 from functools import reduce def git(*args): return subprocess.check_output(["git"] + list(args)) def write_info_file(file_name, build_type="release"): header = '#include "tamaas_info.hh"\n\n' file_content = '' branch = git("rev-parse", "--abbrev-ref", "HEAD")[:-1] commit = git("rev-parse", branch)[:-1] file_content += \ - 'std::string tamaas::TamaasInfo::build_type = "{}";\n'.format(build_type) + 'std::string tamaas::TamaasInfo::build_type = "{}";\n'.format( + build_type) file_content += \ 'std::string tamaas::TamaasInfo::branch = "{}";\n'.format(branch) file_content += \ 'std::string tamaas::TamaasInfo::commit = "{}";\n'.format(commit) diff = git("diff") diff += "|".encode('ascii') + git("diff", "--cached") file_content += \ 'std::string tamaas::TamaasInfo::diff = "{}";\n'.format( base64.b64encode(diff)) remotes = git('remote', '-v') remotes_strings = filter(lambda x: x != '', remotes.splitlines()) remotes_string = reduce(lambda x, y: x + '"{}\\n"\n'.format(y), remotes_strings, "") file_content += \ 'std::string tamaas::TamaasInfo::remotes = {};\n'.format(remotes_string) with open(file_name, 'w') as file: file.write(header + file_content) diff --git a/tests/fftfreq.py b/tests/fftfreq.py index bfa61c3..71aaebf 100644 --- a/tests/fftfreq.py +++ b/tests/fftfreq.py @@ -1,54 +1,53 @@ # -*- 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 __author__ = "Lucas Frérot" import numpy as np from numpy.fft import fftfreq -from tamaas import dtype def frequencies1D(frequencies): n = len(frequencies) frequencies[:] = n * fftfreq(n) def assign2DFrequencies(qx, qy, frequencies): qx, qy = np.meshgrid(qx, qy, indexing='ij') frequencies[:, :, 0] = qx frequencies[:, :, 1] = qy def frequencies2D(frequencies): n = frequencies.shape qx = fftfreq(n[0]) * n[0] qy = fftfreq(n[1]) * n[1] assign2DFrequencies(qx, qy, frequencies) def hfrequencies2D(frequencies): n = frequencies.shape qx = fftfreq(n[0]) * n[0] qy = np.arange(n[1]) assign2DFrequencies(qx, qy, frequencies) diff --git a/tests/test_autocorrelation.py b/tests/test_autocorrelation.py index 6998446..a520be3 100644 --- a/tests/test_autocorrelation.py +++ b/tests/test_autocorrelation.py @@ -1,101 +1,101 @@ # -*- 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 print_function import tamaas as tm import numpy as np from scipy.special import jv import scipy.integrate as integrate def spectrum(r, kl, kr, ks, C, H): if r > ks: return 0 elif r > kr: return C * (r/float(kr))**(-2.*(H+1)) elif r > kl: return C else: return 0 # n is normalization parameter def integrate_bessel(x, spectrum, n): result = 2 * np.pi * \ integrate.quad(lambda y: y*jv(0, y*x) * spectrum(n*y), 0, np.inf)[0] return result def test_autocorrelation(): tm.initialize() grid_size = 256 lambda_l = 8. lambda_r = 8. lambda_s = 1. L = 4 * lambda_l kl = int(L / lambda_l) kr = int(L / lambda_r) ks = int(L / lambda_s) # Computing ACF from integral with Bessel function - l = grid_size // 2 # points on curve - x = np.linspace(0, L/2, l) - bessel_acf = np.zeros(l) + n = grid_size // 2 # points on curve + x = np.linspace(0, L/2, n) + bessel_acf = np.zeros(n) integrate_bessel_v = np.vectorize(integrate_bessel, otypes=[tm.dtype]) bessel_acf = (L / (2*np.pi))**2 *\ integrate_bessel_v(x, lambda x: spectrum(x, kl, kr, ks, 1., 0.8), L / (2*np.pi)) nsamples = 1000 etages = np.zeros((grid_size, grid_size, nsamples)) # Computing ACF with Hu & Tonder algo reimplemented for i in range(nsamples): SG = tm.SurfaceGeneratorFilterFFT() # generate surface SG.getGridSize().assign(grid_size) SG.getHurst().assign(0.8) SG.getQ0().assign(kl) SG.getQ1().assign(kr) SG.getQ2().assign(ks) SG.getRandomSeed().assign(i) SG.Init() a = SG.buildSurface() etages[:, :, i] = a[:, :] - cov = np.zeros(l) + cov = np.zeros(n) - for i in np.arange(l, grid_size, 1): - cov[i-l] = np.sum((etages[l, i, :]-np.mean(etages[l, i, :])) * - (etages[l, l, :]-np.mean(etages[l, l, :])))/nsamples + for i in np.arange(n, grid_size, 1): + cov[i-n] = np.sum((etages[n, i, :]-np.mean(etages[n, i, :])) * + (etages[n, n, :]-np.mean(etages[n, n, :])))/nsamples error = np.sum((cov-bessel_acf)*(cov-bessel_acf)) \ / np.sum(bessel_acf*bessel_acf) print(error) assert error < 1e-1 if __name__ == "__main__": test_autocorrelation() diff --git a/tests/test_bem_grid.py b/tests/test_bem_grid.py index 04bee0f..95d8379 100644 --- a/tests/test_bem_grid.py +++ b/tests/test_bem_grid.py @@ -1,113 +1,112 @@ # -*- 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 print_function import tamaas as tm import numpy as np from numpy.linalg import norm -from numpy.fft import fft2, ifft2 def make_profile(f, profile, k): x = np.linspace(0, 1, profile.shape[0], endpoint=False, dtype=tm.dtype) y = np.linspace(0, 1, profile.shape[1], endpoint=False, dtype=tm.dtype) x, y = np.meshgrid(x, y) profile[:, :] = f(2*k*np.pi*y) def sin_profile(profile, k): make_profile(np.sin, profile, k) def cos_profile(profile, k): make_profile(np.cos, profile, k) def test_bem_grid(): """Checking that all the pipes are clean""" tm.initialize(1) size = 128 k = 1 E = 0.91 nu = 0.3 mu = E / (2*(1+nu)) pressure = np.zeros((size, size), dtype=tm.dtype) surface = np.zeros((size, size), dtype=tm.dtype) westergaard_disp = np.zeros((size, size), dtype=tm.dtype) cos_profile(pressure, k) cos_profile(surface, k) bem = tm.BemGridPolonski(surface) bem.setElasticity(E, nu) bem.computeInfluence() bem.getTractions()[:, :, 2] = pressure[:, :] bem.computeDisplacementsFromTractions() displacements = bem.getDisplacements() # Constructing vertical displacement solution cos_profile(westergaard_disp, k) westergaard_disp *= (1-nu**2)/(E*np.pi*k) westergaard_norm = (1-nu**2)/(E*np.pi*k)*np.sqrt(0.5) * size # Computing error error = norm(displacements[:, :, 2] - westergaard_disp)/westergaard_norm assert error < 1e-10, "Error in normal displacement" # Constructing horizontal displacement solution sin_profile(westergaard_disp, k) westergaard_disp *= (2*nu-1)/(4*np.pi*k*mu) westergaard_norm = np.abs(2*nu-1)/(4*np.pi*k*mu)*np.sqrt(0.5)*size error = norm(displacements[:, :, 1] - westergaard_disp)/westergaard_norm assert error < 1e-10, "Error in horizontal displacement" # Looking at normal vectors normals = np.zeros((size, size, 3), dtype=tm.dtype) normals_norm = np.zeros((size, size, 3), dtype=tm.dtype) sin_profile(normals[:, :, 1], k) normals[:, :, 1] *= 2*np.pi*k normals[:, :, 2] = 1 x = np.linspace(0, 1, size, endpoint=False, dtype=tm.dtype) y = np.linspace(0, 1, size, endpoint=False, dtype=tm.dtype) x, y = np.meshgrid(x, y) normals_norm[:, :, 0] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2) normals_norm[:, :, 1] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2) normals_norm[:, :, 2] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2) normals /= normals_norm bem.computeSurfaceNormals() bem_normals = bem.getSurfaceNormals() error = norm(bem_normals - normals) assert error < 1e-9, "Error in surface normals" tm.finalize() if __name__ == "__main__": test_bem_grid() diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py index a414f63..8a685e7 100644 --- a/tests/test_flood_fill.py +++ b/tests/test_flood_fill.py @@ -1,71 +1,71 @@ # -*- 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 print_function, division import tamaas as tm import numpy as np def test_flood_fill2d(tamaas_fixture): field = np.zeros((18, 18), np.bool) # Cluster of permeter 18 area 13 field[2:4, 3:6] = True - field[4 , 4:6] = True - field[5 , 5 ] = True - field[3:5, 6:8] = True + field[4 , 4:6] = True # noqa + field[5 , 5 ] = True # noqa + field[3:5, 6:8] = True # noqa # Cluster of perimeter 8 area 4 field[14:16, 3:5] = True # Cluster of perimeter 20 area 11 - field[9:11, 8:11 ] = True - field[7:9 , 9 ] = True - field[10 , 11:14] = True + field[9:11, 8:11 ] = True # noqa + field[7:9 , 9 ] = True # noqa + field[10 , 11:14] = True # noqa # Cluster of perimeter 18 area 9 field[3:5, 14:16] = True field[5:10, 15] = True clusters = tm.FloodFill.getClusters(field, False) # Dummy class for clusters class dummy: def __init__(self, area, perimeter): self.area = area self.perimeter = perimeter # Solution ref = [dummy(13, 18), dummy(9, 18), dummy(11, 20), dummy(4, 8)] for x, y in zip(clusters, ref): assert x.getArea() == y.area assert x.getPerimeter() == y.perimeter def test_flood_fill3d(tamaas_fixture): field = np.zeros((5, 5, 5), np.bool) field[2:4, 3, 1:3] = True clusters = tm.FloodFill.getVolumes(field, False) assert len(clusters) == 1 assert clusters[0].getArea() == 4 diff --git a/tests/test_gtest.py b/tests/test_gtest.py index ab1d5ec..5fdcf15 100644 --- a/tests/test_gtest.py +++ b/tests/test_gtest.py @@ -1,34 +1,35 @@ # -*- 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 . import subprocess import os def test_gtest(): me_path = os.path.realpath(__file__) current_dir = os.path.dirname(me_path) my_env = os.environ.copy() my_env['PYTHONPATH'] = current_dir + ":" + os.getenv('PYTHONPATH', "") subprocess.check_call( [os.path.join(current_dir, 'test_gtest_all'), - '--gtest_output=xml:' + os.path.join(current_dir, 'gtest_results.xml')], + '--gtest_output=xml:' + + os.path.join(current_dir, 'gtest_results.xml')], env=my_env ) diff --git a/tests/test_hertz_disp.py b/tests/test_hertz_disp.py index d984c6d..aac1d77 100644 --- a/tests/test_hertz_disp.py +++ b/tests/test_hertz_disp.py @@ -1,130 +1,136 @@ # -*- 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 print_function, division import tamaas as tm import numpy as np def constructHertzProfile(size, curvature): radius = 1. / curvature x = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) y = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) x, y = np.meshgrid(x, y) surface = np.sqrt(radius**2 - x**2 - y**2) surface -= surface.min() return surface.copy() 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)) + disp = np.pi * max_pressure / (4 * contact_size * e_star) \ + * (2 * contact_size**2 - (x**2 + y**2)) return disp.copy() def test_hertz_disp(): tm.initialize() grid_size = 512 curvature = 0.5 effective_modulus = 1. load = 0.12 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.gap, tm.PolonskyKeerRey.gap) solver.solve(load - surface.mean()) tractions = model.getTraction() true_displacements = model.getDisplacement() true_gap = true_displacements - surface pressure = np.mean(tractions) contact_points = np.where(1.0e-10 >= true_gap) contact_area = (np.size(contact_points)/2)/float(grid_size*grid_size) print('The contact area computed with the gap map is '+str(contact_area)) - hertz_contact_size = (3 * pressure / (4 * curvature * effective_modulus))**(1. / 3.) + hertz_contact_size = (3 * pressure + / (4 * curvature * effective_modulus))**(1. / 3.) print('Exact Hertz contact radius is '+str(hertz_contact_size)) hertz_area = np.pi * hertz_contact_size**2 print('Exact Hertz contact area is '+str(hertz_area)) 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 * pressure * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi - pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure + hertz_max_pressure = (6 * pressure * effective_modulus**2 + * curvature ** 2)**(1. / 3.) / np.pi + pressure_error = np.abs(hertz_max_pressure - max_pressure) \ + / hertz_max_pressure print('Exact Hertz max pressure is '+str(hertz_max_pressure)) print('max pressure is '+str(np.max(tractions))) 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 true_gap[i, j] < 1e-12] # 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 += true_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 # Computation of error of displacement in contact zone error = 0. hertz_norm = 0. for index in contact_indexes: - error += (hertz_displacements[index] - true_displacements[index] - correction)**2 + error += (hertz_displacements[index] - + true_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)) - assert area_error < 1e-2 and pressure_error < 1e-2 and displacement_error < 2e-3 + assert area_error < 1e-2 \ + and pressure_error < 1e-2 \ + and displacement_error < 2e-3 tm.finalize() if __name__ == "__main__": test_hertz_disp() diff --git a/tests/test_hertz_kato.py b/tests/test_hertz_kato.py index b224ad6..cb70274 100644 --- a/tests/test_hertz_kato.py +++ b/tests/test_hertz_kato.py @@ -1,117 +1,123 @@ # -*- 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 def constructHertzProfile(size, curvature): radius = 1. / curvature x = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) y = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) x, y = np.meshgrid(x, y) surface = np.sqrt(radius**2 - x**2 - y**2) surface -= surface.mean() return surface.copy() def computeHertzDisplacement(e_star, contact_size, max_pressure, size): x = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) y = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) x, y = np.meshgrid(x, y) - disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - - (x**2 + y**2)) + disp = np.pi * max_pressure / (4 * contact_size * e_star) \ + * (2 * contact_size**2 - (x**2 + y**2)) return disp.copy() def test_hertz_kato(): tm.initialize() grid_size = 512 curvature = 0.1 effective_modulus = 1. load = 0.0001 surface = constructHertzProfile(grid_size, curvature) bem = tm.BemKato(surface) bem.setEffectiveModulus(effective_modulus) bem.setMaxIterations(100000) bem.computeEquilibrium(1e-7, load) tractions = bem.getTractions() displacements = bem.getDisplacements() # 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.) + hertz_contact_size = (3 * load / (4 * curvature + * effective_modulus))**(1. / 3.) 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 + 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 # 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 + 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)) - assert area_error < 6e-3 and pressure_error < 4e-3 and displacement_error < 3e-4 + assert area_error < 6e-3 \ + and pressure_error < 4e-3 \ + and displacement_error < 3e-4 tm.finalize() if __name__ == "__main__": test_hertz_kato() diff --git a/tests/test_tangential.py b/tests/test_tangential.py index 791bdfa..760b3c2 100644 --- a/tests/test_tangential.py +++ b/tests/test_tangential.py @@ -1,102 +1,106 @@ # -*- 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 from numpy.linalg import norm import tamaas as tm def test_pressure(): tm.initialize() E = 1.0 nu = 0.5 mu = 0.5 n = 81 L = 1.0 p_target = np.array([0.5e-4, 0, 2e-4], dtype=tm.dtype) - g_target = np.array([-0.000223, 0, 9.99911], dtype=tm.dtype) + # g_target = np.array([-0.000223, 0, 9.99911], dtype=tm.dtype) x = np.linspace(-L/2, L/2, n, dtype=tm.dtype) y = np.copy(x) x, y = np.meshgrid(x, y, indexing="ij") r_sqrd = (x**2 + y**2) # Spherical contact R = 10.0 surface = R * np.sqrt((r_sqrd < R**2) * (1 - r_sqrd / R**2)) # Creating model - model = tm.ModelFactory.createModel(tm.model_type_surface_2d, [L, L], [n, n]) + model = tm.ModelFactory.createModel(tm.model_type_surface_2d, + [L, L], [n, n]) model.setElasticity(E, nu) p = model.getTraction() # Theoretical solution Estar = E / (1.0 - nu**2) Fn = p_target[2] * L**2 Fx = p_target[0] * L**2 d_theory = np.power(3 / 4 * Fn / Estar / np.sqrt(R), 2/3) a_theory = np.sqrt(R * d_theory) c_theory = a_theory * (1 - Fx / (mu * Fn)) ** (1/3) p0_theory = np.power(6 * Fn * Estar**2 / np.pi**3 / R**2, 1/3) t1_theory = mu * p0_theory t2_theory = t1_theory * c_theory / a_theory - # p_theory = p0_theory * np.sqrt((r_sqrd < a_theory**2) * (1 - r_sqrd / a_theory**2)) - t_theory = t1_theory * np.sqrt((r_sqrd < a_theory**2) * (1 - r_sqrd / a_theory**2)) - t_theory = t_theory - t2_theory * np.sqrt((r_sqrd < c_theory**2) * (1 - r_sqrd / c_theory**2)) + # p_theory = p0_theory * np.sqrt((r_sqrd < a_theory**2) + # * (1 - r_sqrd / a_theory**2)) + t_theory = t1_theory * np.sqrt((r_sqrd < a_theory**2) + * (1 - r_sqrd / a_theory**2)) + t_theory = t_theory - t2_theory * np.sqrt((r_sqrd < c_theory**2) + * (1 - r_sqrd / c_theory**2)) def assert_error(err): error = norm(p[..., 0] - t_theory) / norm(t_theory) print(error) assert error < err # Test Kato solver solver = tm.Kato(model, surface, 1e-12, mu) solver.setMaxIterations(1000) solver.solve(p_target, 100) assert_error(4e-2) # solver.solveRelaxed(g_target) # assert_error(1e-2) # # Test BeckTeboulle solver # solver = tm.BeckTeboulle(model, surface, 1e-12, mu) # solver.solve(g_target) # assert_error(1e-2) # Test Condat solver solver = tm.Condat(model, surface, 1e-12, mu) - solver.setMaxIterations(5000) # or 10000 + solver.setMaxIterations(5000) # or 10000 solver.solve(p_target) - assert_error(7e-2) # 4e-2 for 10000 iterations + assert_error(7e-2) # 4e-2 for 10000 iterations # Test tangential Polonsky Keer solver solver = tm.PolonskyKeerTan(model, surface, 1e-12, mu) solver.setMaxIterations(1000) solver.solve(p_target) assert_error(4e-2) tm.finalize() if __name__ == "__main__": test_pressure() diff --git a/tests/test_voigt.py b/tests/test_voigt.py index e1cdfd5..8abfb5d 100644 --- a/tests/test_voigt.py +++ b/tests/test_voigt.py @@ -1,34 +1,34 @@ # -*- 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 . import tamaas as tm import numpy as np def test_voigt(tamaas_fixture): shape = [3, 3, 3] stress = np.ones(shape + [9], dtype=tm.dtype) - #voigt = np.zeros(shape + [6]) + # voigt = np.zeros(shape + [6]) voigt = tm.to_voigt(stress) sol = np.ones(shape + [6], dtype=tm.dtype) sol[:, :, :, 3:] = np.sqrt(2) print(voigt) print(sol) assert np.all(np.abs(sol - voigt) < 1e-15) diff --git a/tests/test_westergaard.py b/tests/test_westergaard.py index 900ad4c..e66ba23 100644 --- a/tests/test_westergaard.py +++ b/tests/test_westergaard.py @@ -1,80 +1,78 @@ # -*- 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 sys -import numpy as np from numpy.linalg import norm import tamaas as tm import pytest @pytest.fixture(scope="module", params=[tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.gap]) def pkr(westergaard, request): loads = { tm.PolonskyKeerRey.pressure: westergaard.load, tm.PolonskyKeerRey.gap: 0.06697415181446396 } model = tm.ModelFactory.createModel(tm.model_type.basic_1d, [1.], [westergaard.n]) model.setElasticity(westergaard.e_star, 0) solver = tm.PolonskyKeerRey(model, westergaard.surface, 1e-12, request.param, request.param) solver.setMaxIterations(5000) solver.solve(loads[request.param]) return model, westergaard def test_energy(pkr): model, sol = pkr traction, displacement = model.getTraction(), model.getDisplacement() error = norm((traction - sol.pressure)*(displacement - sol.displacement)) error /= norm(sol.pressure * sol.displacement) assert error < 4e-6 def test_pkr_multi(westergaard): model_basic = tm.ModelFactory.createModel(tm.model_type.basic_1d, [1.], [westergaard.n]) model_volume = tm.ModelFactory.createModel(tm.model_type.volume_1d, [1., 1.], [20, westergaard.n]) pressures = {} for model in [model_basic, model_volume]: print(model) solver = tm.PolonskyKeerRey(model, westergaard.surface, 1e-12, tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.pressure) solver.solve(westergaard.load) pressures[model] = model.getTraction() error = norm(pressures[model_basic] - pressures[model_volume][:, 1])\ / westergaard.load assert error < 1e-16 if __name__ == "__main__": print("Test is meant to run with pytest")