diff --git a/INFOS.py b/INFOS.py
index 13128c0..3d6497f 100755
--- a/INFOS.py
+++ b/INFOS.py
@@ -1,55 +1,55 @@
#!/usr/bin/env python3
"""Defines the information to be used throughout the builds."""
# -*- mode:python; coding: utf-8 -*-
from collections import namedtuple
import versioneer
TamaasInfo = namedtuple('TamaasInfo',
['version',
'authors',
'maintainer',
'email',
'copyright',
'description',
'license'])
TAMAAS_INFOS = TamaasInfo(
version=versioneer.get_version(),
authors=([
u'Lucas Frérot',
'Guillaume Anciaux',
'Valentine Rey',
'Son Pham-Ba',
u'Jean-François Molinari'
]),
maintainer=u'Lucas Frérot',
email='lucas.frerot@imtek.uni-freiburg.de',
copyright=(
- "Copyright (©) 2016-2022 EPFL "
+ "Copyright (©) 2016-2023 EPFL "
"(École Polytechnique Fédérale de Lausanne), "
"Laboratory (LSMS - Laboratoire de Simulation en "
"Mécanique des Solides)\\n"
- "Copyright (©) 2020-2022 Lucas Frérot"
+ "Copyright (©) 2020-2023 Lucas Frérot"
),
description='A high-performance library for periodic rough surface contact',
license="SPDX-License-Identifier: AGPL-3.0-or-later",
)
def main():
import argparse
parser = argparse.ArgumentParser(description="Print Tamaas info")
parser.add_argument("--version", action="store_true")
args = parser.parse_args()
if args.version:
print(TAMAAS_INFOS.version)
if __name__ == "__main__":
main()
diff --git a/SConstruct b/SConstruct
index e060a42..f4e8961 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,476 +1,476 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
from __future__ import print_function
import sys
import os
from subprocess import check_output
from warnings import warn
# Import below not strictly necessary, but good for pep8
from SCons.Script import (
EnsurePythonVersion,
EnsureSConsVersion,
Help,
Environment,
Variables,
EnumVariable,
PathVariable,
BoolVariable,
ListVariable,
Split,
Export,
)
from SCons.Errors import StopError
from SCons import __version__ as scons_version
from version import get_git_subst
from detect import (FindFFTW, FindBoost, FindThrust, FindCuda, FindExpolit,
FindPybind11)
from INFOS import TAMAAS_INFOS
# ------------------------------------------------------------------------------
EnsurePythonVersion(3, 6)
EnsureSConsVersion(3, 0)
# ------------------------------------------------------------------------------
def detect_dependencies(env):
"Detect all dependencies"
fftw_comp = {
'omp': ['omp'],
'threads': ['threads'],
'none': [],
}
fftw_components = fftw_comp[env['fftw_threads']]
if main_env['use_mpi']:
fftw_components.append('mpi')
if main_env['use_fftw']:
FindFFTW(env, fftw_components, precision=env['real_type'])
if main_env['use_cuda']:
FindCuda(env)
FindBoost(env, [
'boost/version.hpp',
'boost/preprocessor/seq.hpp',
'boost/variant.hpp',
])
# Use thrust shipped with cuda if cuda is requested
thrust_var = 'CUDA_ROOT' if env['use_cuda'] else 'THRUST_ROOT'
FindThrust(env, env['backend'], thrust_var)
if env['build_python']:
FindPybind11(env)
FindExpolit(env)
def subdir(env, dir):
"Building a sub-directory"
return env.SConscript(env.File('SConscript', dir),
variant_dir=env.Dir(dir, env['build_dir']),
duplicate=True)
def print_build_info(env):
info = ("-- Tamaas ${version}\n" + "-- SCons {} (Python {}.{})\n".format(
scons_version, sys.version_info.major, sys.version_info.minor) +
"-- Build type: ${build_type}\n" +
"-- Thrust backend: ${backend}\n" +
("-- FFTW threads: ${fftw_threads}\n" if env['use_fftw'] else '') +
"-- MPI: ${use_mpi}\n" + "-- Build directory: ${build_dir}\n" +
("-- Python version (bindings): $py_version"
if env['build_python'] else ''))
print(env.subst(info))
# ------------------------------------------------------------------------------
# Main compilation
# ------------------------------------------------------------------------------
# Compilation colors
colors = {
'cyan': '\033[96m',
'purple': '\033[95m',
'blue': '\033[94m',
'green': '\033[92m',
'yellow': '\033[93m',
'gray': '\033[38;5;8m',
'orange': '\033[38;5;208m',
'red': '\033[91m',
'end': '\033[0m'
}
# Inherit all environment variables (for CXX detection, etc.)
main_env = Environment(ENV=os.environ, )
# Set tamaas information
for k, v in TAMAAS_INFOS._asdict().items():
main_env[k] = v
main_env['COLOR_DICT'] = colors
main_env.AddMethod(subdir, 'SubDirectory')
# Build variables
vars = Variables('build-setup.conf')
vars.AddVariables(
EnumVariable('build_type',
'Build type',
'release',
allowed_values=('release', 'profiling', 'debug'),
ignorecase=2),
EnumVariable('backend',
'Thrust backend',
'cpp',
allowed_values=('cpp', 'omp', 'tbb', 'cuda'),
ignorecase=2),
EnumVariable('fftw_threads',
'Threads FFTW library preference',
'none',
allowed_values=('omp', 'threads', 'none'),
ignorecase=2),
EnumVariable('sanitizer',
'Sanitizer type',
'none',
allowed_values=('none', 'memory', 'leaks', 'address'),
ignorecase=2),
PathVariable('prefix', 'Prefix where to install', '/usr/local',
PathVariable.PathAccept),
# Dependencies paths
PathVariable('FFTW_ROOT', 'FFTW custom path', os.getenv('FFTW_ROOT', ''),
PathVariable.PathAccept),
PathVariable('THRUST_ROOT', 'Thrust custom path',
os.getenv('THRUST_ROOT', ''), PathVariable.PathAccept),
PathVariable('BOOST_ROOT', 'Boost custom path',
os.getenv('BOOST_ROOT', ''), PathVariable.PathAccept),
PathVariable('CUDA_ROOT', 'Cuda custom path', os.getenv('CUDA_ROOT', ''),
PathVariable.PathAccept),
PathVariable('GTEST_ROOT', 'Googletest custom path',
os.getenv('GTEST_ROOT', ''), PathVariable.PathAccept),
PathVariable('PYBIND11_ROOT', 'Pybind11 custom path',
os.getenv('PYBIND11_ROOT', ''), PathVariable.PathAccept),
# Executables
('CXX', 'Compiler', os.getenv('CXX', 'g++')),
('MPICXX', 'MPI Compiler wrapper', os.getenv('MPICXX', 'mpicxx')),
('py_exec', 'Python executable', 'python3'),
# Compiler flags
('CXXFLAGS', 'C++ compiler flags', os.getenv('CXXFLAGS', "")),
# Cosmetic
BoolVariable('verbose', 'Activate verbosity',
os.getenv('VERBOSE') in {'1', 'True', 'true'}),
BoolVariable('color', 'Color the non-verbose compilation output', False),
# Tamaas components
BoolVariable('build_tests', 'Build test suite', False),
BoolVariable('build_python', 'Build python wrapper', True),
# Documentation
ListVariable('doc_builders',
'Generated documentation formats',
default='html',
names=Split("html man")), # TODO include latex
# Dependencies
BoolVariable('use_mpi', 'Builds multi-process parallelism', False),
# Distribution options
BoolVariable('strip_info', 'Strip binary of added information', True),
BoolVariable('build_static_lib', "Build a static libTamaas", False),
# Type variables
EnumVariable('real_type',
'Type for real precision variables',
'double',
allowed_values=('double', 'long double')),
EnumVariable('integer_type',
'Type for integer variables',
'int',
allowed_values=('int', 'long')),
)
# Set variables of environment
vars.Update(main_env)
help_text = vars.GenerateHelpText(main_env)
help_text += """
Commands:
scons [build] [options]... Compile Tamaas (and additional modules/tests)
scons install [prefix=/your/prefix] [options]... Install Tamaas to prefix
scons dev Install symlink to Tamaas python module (useful to development purposes)
scons test Run tests with pytest
scons doc Compile documentation with Doxygen and Sphinx+Breathe
scons archive Create a gzipped archive from source
""" # noqa
Help(help_text)
# Save all options, not just those that differ from default
with open('build-setup.conf', 'w') as setup:
for option in vars.options:
setup.write("# " + option.help.replace('\n', '\n# ') + "\n")
setup.write("{} = '{}'\n".format(option.key, main_env[option.key]))
# Printing unknown variables
if vars.UnknownVariables():
warn(f'unknown variables provided:\n\t{vars.UnknownVariables()}')
main_env['should_configure'] = \
not main_env.GetOption('clean') and not main_env.GetOption('help')
build_type = main_env['build_type']
build_dir = 'build-${build_type}'
main_env['build_dir'] = main_env.Dir(build_dir)
# Setting up the python name with version
if main_env['build_python']:
args = (main_env.subst("${py_exec} -c").split() + [
"from sysconfig import get_python_version;"
"print(get_python_version())"
])
main_env['py_version'] = bytes(check_output(args)).decode()
verbose = main_env['verbose']
# Remove colors if not set
if not main_env['color']:
for key in colors:
colors[key] = ''
if not verbose:
main_env['CXXCOMSTR'] = main_env['SHCXXCOMSTR'] = \
u'{0}[Compiling ($SHCXX)] {1}$SOURCE'.format(colors['green'],
colors['end'])
main_env['LINKCOMSTR'] = main_env['SHLINKCOMSTR'] = \
u'{0}[Linking] {1}$TARGET'.format(colors['purple'],
colors['end'])
main_env['ARCOMSTR'] = u'{}[Ar]{} $TARGET'.format(colors['purple'],
colors['end'])
main_env['RANLIBCOMSTR'] = \
u'{}[Randlib]{} $TARGET'.format(colors['purple'],
colors['end'])
main_env['PRINT_CMD_LINE_FUNC'] = pretty_cmd_print
main_env['INSTALLSTR'] = \
u'{}[Installing] {}$SOURCE to $TARGET'.format(colors['blue'],
colors['end'])
# Include paths
main_env.AppendUnique(CPPPATH=[
'#/src', '#/src/core', '#/src/surface', '#/src/percolation', '#/src/model',
'#/src/model/elasto_plastic', '#/src/solvers', '#/src/gpu', '#/python',
'#/third-party/expolit/include'
])
# Changing the shared object extension
main_env['SHOBJSUFFIX'] = '.o'
# Variables for clarity
main_env['use_cuda'] = main_env['backend'] == "cuda"
main_env['use_fftw'] = not main_env['use_cuda']
main_env['use_mpi'] = main_env['use_mpi'] and not main_env['use_cuda']
if not main_env['use_fftw']:
main_env['fftw_threads'] = 'none'
# Back to gcc if cuda is activated
if main_env['use_cuda'] and "g++" not in main_env['CXX']:
raise StopError('GCC should be used when compiling with CUDA')
# Printing some build infos
if main_env['should_configure']:
print_build_info(main_env)
# OpenMP flags - compiler dependent
omp_flags = {
"g++": ["-fopenmp"],
"clang++": ["-fopenmp"],
"icpc": ["-qopenmp"]
}
def cxx_alias(cxx):
for k in omp_flags.keys():
if k in cxx:
return k
raise StopError('Unsupported compiler: ' + cxx)
cxx = cxx_alias(main_env['CXX'])
# Setting main compilation flags
main_env['CXXFLAGS'] = Split(main_env['CXXFLAGS'])
main_env['LINKFLAGS'] = main_env['CXXFLAGS']
main_env.AppendUnique(
CXXFLAGS=Split('-std=c++14 -Wall -Wextra'),
CPPDEFINES={
'TAMAAS_LOOP_BACKEND': 'TAMAAS_LOOP_BACKEND_${backend.upper()}',
'TAMAAS_FFTW_BACKEND': 'TAMAAS_FFTW_BACKEND_${fftw_threads.upper()}'
},
)
if main_env['backend'] != 'cuda':
main_env.AppendUnique(CXXFLAGS=['-pedantic'])
# Adding OpenMP flags
if main_env['backend'] == 'omp':
main_env.AppendUnique(CXXFLAGS=omp_flags[cxx])
main_env.AppendUnique(LINKFLAGS=omp_flags[cxx])
else:
main_env.AppendUnique(CXXFLAGS=['-Wno-unknown-pragmas'])
# Correct bug in clang?
if main_env['backend'] == 'omp' and cxx == "clang++":
main_env.AppendUnique(LIBS=["atomic"])
elif main_env['backend'] == 'tbb':
main_env.AppendUnique(LIBS=['tbb'])
# Manage MPI compiler
if main_env['use_mpi']:
main_env['CXX'] = '$MPICXX'
main_env.AppendUnique(CPPDEFINES=['TAMAAS_USE_MPI'])
main_env.AppendUnique(CXXFLAGS=['-Wno-cast-function-type'])
# Flags and options
if main_env['build_type'] == 'debug':
main_env.AppendUnique(CPPDEFINES=['TAMAAS_DEBUG'])
# Define the scalar types
main_env.AppendUnique(CPPDEFINES={
'TAMAAS_REAL_TYPE': '${real_type}',
'TAMAAS_INT_TYPE': '${integer_type}'
})
# Compilation flags
cxxflags_dict = {
"debug": Split("-g -O0"),
"profiling": Split("-g -O3 -fno-omit-frame-pointer"),
"release": Split("-O3")
}
if main_env['sanitizer'] != 'none':
if main_env['backend'] == 'cuda':
raise StopError("Sanitizers with cuda are not yet supported!")
cxxflags_dict[build_type].append('-fsanitize=${sanitizer}')
main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type])
main_env.AppendUnique(SHLINKFLAGS=cxxflags_dict[build_type])
main_env.AppendUnique(LINKFLAGS=cxxflags_dict[build_type])
if main_env['should_configure']:
basic_checks(main_env)
detect_dependencies(main_env)
# Writing information file
main_env.Tool('textfile')
main_env['SUBST_DICT'] = get_git_subst()
# Empty values if requested
if main_env['strip_info']:
for k in main_env['SUBST_DICT']:
main_env['SUBST_DICT'][k] = ""
# Substitution of environment file
main_env['SUBST_DICT'].update({
'@build_type@': '$build_type',
'@build_dir@': '${build_dir.abspath}',
'@build_version@': '$version',
'@backend@': '$backend',
})
# Environment file content
env_content = """export PYTHONPATH=@build_dir@/python:$$PYTHONPATH
export LD_LIBRARY_PATH=@build_dir@/src:$$LD_LIBRARY_PATH
"""
# Writing environment file
env_file = main_env.Textfile(
main_env.File('tamaas_environment.sh', main_env['build_dir']), env_content)
# Default targets
build_targets = ['build-cpp', env_file]
install_targets = ['install-lib']
if main_env._get_major_minor_revision(scons_version)[0] >= 4:
main_env.Tool('compilation_db')
build_targets.append(
main_env.CompilationDatabase(PRINT_CMD_LINE_FUNC=pretty_cmd_print
if not main_env['verbose'] else None))
# Building Tamaas library
Export('main_env')
main_env.SubDirectory('src')
# Building Tamaas extra components
for dir in ['python', 'tests']:
if main_env['build_{}'.format(dir)] and not main_env.GetOption('help'):
main_env.SubDirectory(dir)
build_targets.append('build-{}'.format(dir))
# Building API + Sphinx documentation if requested
main_env.SubDirectory('doc')
main_env.Alias('doc', 'build-doc')
install_targets.append('install-doc')
# Define dummy dev command when python is deactivated
if not main_env['build_python']:
dummy_command(
main_env, 'dev', 'Command "dev" does not do anything' +
' without python activated ("build_python=True")')
else:
install_targets.append('install-python')
# Define dummy test command when tests are deactivated
if not main_env['build_tests']:
dummy_command(
main_env, 'test', 'Command "test" does not do anything' +
' without tests activated ("build_tests=True")')
# Definition of target aliases, a.k.a. sub-commands
main_env.Alias('build', build_targets)
# Define proper install targets
main_env.Alias('install', install_targets)
# Default target is to build stuff
main_env.Default('build')
# Building a tar archive
archive = main_env.Command(
'tamaas-${version}.tar.gz',
'',
('git archive '
'--format=tar.gz '
'--prefix=tamaas/ '
'-o $TARGET HEAD'),
)
main_env.Alias('archive', archive)
diff --git a/doc/SConscript b/doc/SConscript
index b180071..122c3d0 100644
--- a/doc/SConscript
+++ b/doc/SConscript
@@ -1,107 +1,107 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 SCons.Script import Import, Glob, Dir, File
def add_sources(source_list, dirname, fnames):
source_list += Dir(dirname).glob('*.hh')
source_list += Dir(dirname).glob('*.cpp')
Import('main_env')
Import('libTamaas')
doc_env = main_env.Clone()
doc_env['man_section'] = 7
doc_targets = []
# Generate Doxygen API documentation
doc_env.Tool('doxygen')
if doc_env['has_doxygen']:
# Generating Doxyfile
doxygen_verbose = {True: "NO", False: "YES"}
doxyfile_target = doc_env.Substfile('doxygen/Doxyfile',
'doxygen/Doxyfile.in',
SUBST_DICT={
'@version@': '$version',
'@build_dir@': Dir('doxygen').path,
'@logo@': File('icon.svg').path,
'@src_dir@': Dir('#src').abspath,
'@verbose@':
doxygen_verbose[doc_env["verbose"]],
})
doxygen_target = doc_env.Doxygen('doxygen/xml/index.xml',
[doxyfile_target, 'icon.svg', 'icon.ico'])
# Adding all source files as dependencies
sources = []
Dir('#src').walk(add_sources, sources)
doc_env.Depends(doxygen_target, sources)
doc_targets.append(doxygen_target)
# Generate Sphinx User documentation
sphinx_targets_map = {
"html": "sphinx/html/index.html",
"man": "sphinx/tamaas.${man_section}",
"latex": "sphinx/latex/Tamaas.tex"
}
doc_env.Tool('sphinx')
sphinx_targets = None
if doc_env['has_sphinx'] and len(doc_targets) != 0:
doc_env['IMPLICIT_COMMAND_DEPENDENCIES'] = 0
sphinx_sources = [Glob('sphinx/source/*'), Glob('sphinx/source/figures/*')]
sphinx_targets = {
builder: doc_env.Sphinx(sphinx_targets_map[builder], sphinx_sources)
for builder in doc_env['doc_builders']
}
for target in sphinx_targets.values():
doc_env.Depends(target, doxygen_target)
if main_env['build_python']:
doc_env.Depends(target, 'build-python')
doc_targets += list(sphinx_targets.values())
# Alias for both docs
main_env.Alias('build-doc', doc_targets)
# Install target for documentation
share = Dir('share', doc_env['prefix'])
doc_install = []
sphinx_prefix_map = {}
if sphinx_targets is not None:
if "html" in doc_env['doc_builders']:
doc_install.append(doc_env.Install(
target=share.Dir('doc').Dir('tamaas'),
source=sphinx_targets['html'][0].dir))
if "man" in doc_env['doc_builders']:
doc_install.append(doc_env.Install(
target=share.Dir('man') .Dir(doc_env.subst('man${man_section}')),
source=sphinx_targets['man']))
main_env.Alias('install-doc', doc_install)
diff --git a/examples/adhesion.py b/examples/adhesion.py
index 808f153..4ca2a8a 100644
--- a/examples/adhesion.py
+++ b/examples/adhesion.py
@@ -1,123 +1,123 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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
import numpy as np
import argparse
from tamaas.utils import publications
from matplotlib.colors import ListedColormap
class AdhesionPython(tm.Functional):
"""
Functional class that extends a C++ class and implements the virtual
methods
"""
def __init__(self, rho, gamma):
tm.Functional.__init__(self)
self.rho = rho
self.gamma = gamma
def computeF(self, gap, pressure):
return -self.gamma * np.sum(np.exp(-gap / self.rho))
def computeGradF(self, gap, gradient):
gradient += self.gamma * np.exp(-gap / self.rho) / self.rho
parser = argparse.ArgumentParser()
parser.add_argument('--local-functional', dest="py_adh", action="store_true",
help="use the adhesion functional written in python")
args = parser.parse_args()
# Surface size
n = 1024
# Surface generator
sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 1
# Spectrum
sg.spectrum = tm.Isopowerlaw2D()
# Parameters
sg.spectrum.q0 = 16
sg.spectrum.q1 = 16
sg.spectrum.q2 = 64
sg.spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
surface /= n
plt.imshow(surface)
plt.title('Rough 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.gap,
tm.PolonskyKeerRey.gap)
adhesion_params = {
"rho": 2e-3,
"surface_energy": 2e-5
}
# Use the python derived from C++ functional class
if args.py_adh:
adhesion = AdhesionPython(adhesion_params["rho"],
adhesion_params["surface_energy"])
# Use the C++ class
else:
adhesion = tm.ExponentialAdhesionFunctional(surface)
adhesion.parameters = adhesion_params
solver.addFunctionalTerm(adhesion)
# Solve for target pressure
g_target = 5e-2
solver.solve(g_target)
tractions = model.traction
plt.figure()
plt.imshow(tractions)
plt.colorbar()
plt.title('Contact tractions')
plt.figure()
zones = np.zeros_like(tractions)
tol = 1e-6
zones[tractions > tol] = 1
zones[tractions < -tol] = -1
plt.imshow(zones, cmap=ListedColormap(['white', 'gray', 'black']))
plt.colorbar(ticks=[-2/3, 0, 2/3]).set_ticklabels([
'Adhesion', 'No Contact', 'Contact'
])
plt.title('Contact and Adhesion Zones')
plt.show()
publications()
diff --git a/examples/pipe_tools/contact b/examples/pipe_tools/contact
index 634c3a1..59d339c 100755
--- a/examples/pipe_tools/contact
+++ b/examples/pipe_tools/contact
@@ -1,58 +1,58 @@
#!/usr/bin/env python3
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
"""
Read from stdin a surface profile and solve an elastic contact problem with
load given as an argument.
"""
import argparse
import sys
import tamaas as tm
import numpy as np
from tamaas.dumpers import NumpyDumper
__author__ = "Lucas Frérot"
__copyright__ = (
- "Copyright (©) 2019-2022, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2023, EPFL (École Polytechnique Fédérale de Lausanne),"
"\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "SPDX-License-Identifier: AGPL-3.0-or-later"
tm.set_log_level(tm.LogLevel.error)
parser = argparse.ArgumentParser(
description="Compute the elastic contact solution with a given surface")
parser.add_argument("--input", "-i",
help="Rough surface file (default read from stdin)")
parser.add_argument("--tol",
type=float,
default=1e-12,
help="Solver tolerance")
parser.add_argument("load",
type=float,
help="Applied average pressure")
args = parser.parse_args()
if not args.input:
input = sys.stdin
else:
input = args.input
surface = np.loadtxt(input)
discretization = surface.shape
system_size = [1., 1.]
model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
system_size, discretization)
solver = tm.PolonskyKeerRey(model, surface, args.tol)
solver.solve(args.load)
dumper = NumpyDumper('numpy', 'traction', 'displacement')
dumper.dump_to_file(sys.stdout.buffer, model)
diff --git a/examples/pipe_tools/plot b/examples/pipe_tools/plot
index b6c3f6e..6e2d6fc 100755
--- a/examples/pipe_tools/plot
+++ b/examples/pipe_tools/plot
@@ -1,49 +1,49 @@
#!/usr/bin/env python3
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
"""
Read Numpy data from standard input, plot contact tractions and displacements.
"""
import sys
import io
import matplotlib.pyplot as plt
import numpy as np
__author__ = "Lucas Frérot"
__copyright__ = (
- "Copyright (©) 2019-2022, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2023, EPFL (École Polytechnique Fédérale de Lausanne),"
"\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "SPDX-License-Identifier: AGPL-3.0-or-later"
def load_stream(stream):
"""
Load numpy from binary stream (allows piping)
Code from
https://gist.github.com/CMCDragonkai/3c99fd4aabc8278b9e17f50494fcc30a
"""
np_magic = stream.read(6)
# use the sys.stdin.buffer to read binary data
np_data = stream.read()
# read it all into an io.BytesIO object
return io.BytesIO(np_magic + np_data)
fig, (ax_traction, ax_displacement) = plt.subplots(1, 2)
ax_traction.set_title('Traction')
ax_displacement.set_title('Displacement')
with load_stream(sys.stdin.buffer) as f_np:
data = np.load(f_np)
ax_traction.imshow(data['traction'])
ax_displacement.imshow(data['displacement'])
fig.set_size_inches(10, 6)
fig.tight_layout()
plt.show()
diff --git a/examples/pipe_tools/surface b/examples/pipe_tools/surface
index 0997600..87db867 100755
--- a/examples/pipe_tools/surface
+++ b/examples/pipe_tools/surface
@@ -1,77 +1,77 @@
#!/usr/bin/env python3
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
"""
Create a random self-affine rough surface from command line parameters.
"""
import argparse
import sys
import time
import tamaas as tm
import numpy as np
__author__ = "Lucas Frérot"
__copyright__ = (
- "Copyright (©) 2019-2022, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2023, EPFL (École Polytechnique Fédérale de Lausanne),"
"\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "SPDX-License-Identifier: AGPL-3.0-or-later"
parser = argparse.ArgumentParser(
description="Generate a self-affine rough surface"
)
parser.add_argument("--cutoffs", "-K",
nargs=3,
type=int,
help="Long, rolloff and short wavelength cutoffs",
metavar=('k_l', 'k_r', 'k_s'),
required=True)
parser.add_argument("--sizes",
nargs=2,
type=int,
help="Number of points",
metavar=('nx', 'ny'),
required=True)
parser.add_argument("--hurst", "-H",
type=float,
help="Hurst exponent",
required=True)
parser.add_argument("--rms",
type=float,
help="Root-mean-square of slopes",
default=1.)
parser.add_argument("--seed",
type=int,
help="Random seed",
default=int(time.time()))
parser.add_argument("--generator",
help="Generation method",
choices=('random_phase', 'filter'),
default='random_phase')
parser.add_argument("--output", "-o",
help="Output file name (compressed if .gz)")
args = parser.parse_args()
if args.generator == 'random_phase':
generator = tm.SurfaceGeneratorRandomPhase2D(args.sizes)
elif args.generator == 'filter':
generator = tm.SurfaceGeneratorFilter2D(args.sizes)
else:
raise ValueError('Unknown generator method {}'.format(args.generator))
generator.spectrum = tm.Isopowerlaw2D()
generator.spectrum.q0 = args.cutoffs[0]
generator.spectrum.q1 = args.cutoffs[1]
generator.spectrum.q2 = args.cutoffs[2]
generator.spectrum.hurst = args.hurst
generator.random_seed = args.seed
surface = generator.buildSurface() / generator.spectrum.rmsSlopes() * args.rms
output = args.output if args.output is not None else sys.stdout
np.savetxt(output, surface)
diff --git a/examples/plasticity.py b/examples/plasticity.py
index 8c48f80..4e4ba95 100644
--- a/examples/plasticity.py
+++ b/examples/plasticity.py
@@ -1,85 +1,85 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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 numpy as np
import tamaas as tm
from tamaas.dumpers import H5Dumper as Dumper
from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver, ToleranceManager
from tamaas.utils import publications
# Definition of modeled domain
model_type = tm.model_type.volume_2d
discretization = [32, 51, 51]
flat_domain = [1, 1]
system_size = [0.5] + flat_domain
# Creation of model
model = tm.ModelFactory.createModel(model_type,
system_size,
discretization)
model.E = 1.
model.nu = 0.3
# Setup for plasticity
residual = tm.ModelFactory.createResidual(model,
sigma_y=0.1 * model.E,
hardening=0.01 * model.E)
# Possibly change integration method
# residual.setIntegrationMethod(tm.integration_method.cutoff, 1e-12)
# Setup non-linear solver with variable tolerance
epsolver = ToleranceManager(1e-5, 1e-9, 1/4)(Solver)(residual)
# Setup for contact
x = np.linspace(0, system_size[1], discretization[1],
endpoint=False, dtype=tm.dtype)
y = np.linspace(0, system_size[2], discretization[2],
endpoint=False, dtype=tm.dtype)
xx, yy = np.meshgrid(x, y, indexing='ij')
R = 0.2
surface = -((xx - flat_domain[0] / 2)**2 +
(yy - flat_domain[1] / 2)**2) / (2 * R)
# Scatter surface across MPI processes
local_surface = tm.mpi.scatter(surface)
csolver = tm.PolonskyKeerRey(model, local_surface, 1e-12,
tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.pressure)
# EPIC setup
epic = tm.EPICSolver(csolver, epsolver, 1e-7)
# Dumper
dumper_helper = Dumper('hertz', 'displacement', 'stress', 'plastic_strain')
model.addDumper(dumper_helper)
loads = np.linspace(0.001, 0.005, 3)
for i, load in enumerate(loads):
epic.acceleratedSolve(load)
model.dump()
tm.Logger().get(tm.LogLevel.info) \
<< "---> Solved load step {}/{}".format(i+1, len(loads))
# Print list of relevant publications
publications()
diff --git a/examples/rough_contact.py b/examples/rough_contact.py
index 878c191..28ff8f6 100644
--- a/examples/rough_contact.py
+++ b/examples/rough_contact.py
@@ -1,66 +1,66 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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
from tamaas.utils import publications
# Initialize threads and fftw
tm.set_log_level(tm.LogLevel.info) # Show progression of solver
# Surface size
n = 512
# Surface generator
sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 1
# Spectrum
sg.spectrum = tm.Isopowerlaw2D()
# Parameters
sg.spectrum.q0 = 16
sg.spectrum.q1 = 16
sg.spectrum.q2 = 64
sg.spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
surface /= tm.Statistics2D.computeSpectralRMSSlope(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)
# Solve for target pressure
p_target = 0.1
solver.solve(p_target)
plt.figure()
plt.imshow(model.traction)
plt.title('Contact tractions')
print(model.traction.mean())
plt.show()
publications()
diff --git a/examples/saturated.py b/examples/saturated.py
index c48ab5e..48f0f04 100644
--- a/examples/saturated.py
+++ b/examples/saturated.py
@@ -1,68 +1,68 @@
# -*- coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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
import matplotlib.pyplot as plt
from tamaas.utils import publications
grid_size = 256
load = 1.6
p_sat = 100
iso = tm.Isopowerlaw2D()
iso.q0 = 4
iso.q1 = 4
iso.q2 = 16
iso.hurst = 0.8
sg = tm.SurfaceGeneratorFilter2D([grid_size, grid_size])
sg.random_seed = 2
sg.spectrum = iso
surface = sg.buildSurface()
surface -= np.max(surface)
model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
[1., 1.], [grid_size, grid_size])
model.E = 1.
model.nu = 0
esolver = tm.PolonskyKeerRey(model, surface, 1e-12,
tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.pressure)
esolver.solve(load)
elastic_tractions = model.traction.copy()
plt.imshow(elastic_tractions)
plt.title('Elastic contact tractions')
plt.colorbar()
model.traction[:] = 0.
solver = tm.KatoSaturated(model, surface, 1e-12, p_sat)
solver.dump_freq = 1
solver.max_iter = 200
solver.solve(load)
plt.figure()
plt.imshow(model.traction)
plt.title('Saturated contact tractions')
plt.colorbar()
plt.show()
publications()
diff --git a/examples/scipy_penalty.py b/examples/scipy_penalty.py
index efded2e..31edbbb 100644
--- a/examples/scipy_penalty.py
+++ b/examples/scipy_penalty.py
@@ -1,119 +1,119 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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
import matplotlib.pyplot as plt
import scipy.optimize
from tamaas.utils import hertz_surface
class GapPenaltyFunctional(tm.Functional):
"""Quadratic penalty on negative gap functional."""
def __init__(self, penalty):
super().__init__()
self.penalty = penalty
@staticmethod
def clip(gap):
"""Clip negative gap."""
return np.clip(gap, -np.inf, 0)
def computeF(self, gap, dual):
"""Compute penalty."""
g = self.clip(gap)
return 0.5 / self.penalty * np.vdot(g, g) / g.size
def computeGradF(self, gap, gradient):
"""Compute gradient."""
gradient[:] += (1 / self.penalty) * self.clip(gap)
class ScipyContactSolver(tm.ContactSolver):
"""Gap-based penalty solver using https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html"""
def __init__(self, model, surface, tolerance, penalty):
super().__init__(model, surface, tolerance)
if model.type != tm.model_type.basic_2d:
raise TypeError("model type is not basic_2d")
self.shape = model.boundary_shape
self.N = np.prod(self.shape)
model['gap'] = np.zeros(list(self.shape) + [1])
model.be_engine.registerDirichlet()
self.elastic = tm.ElasticFunctionalGap(
self.model.operators['Westergaard::dirichlet'], self.surface)
self.penalty = GapPenaltyFunctional(penalty)
self.addFunctionalTerm(self.elastic)
self.addFunctionalTerm(self.penalty)
def solve(self, mean_gap):
"""Solve with mean gap constraint."""
def cost(gap):
gap = gap.reshape(self.shape)
grad = np.zeros_like(gap).reshape(self.shape)
self.functional.computeGradF(gap, grad)
return self.functional.computeF(gap, grad)
def jac(gap):
gap = gap.reshape(self.shape)
grad = np.zeros_like(gap).reshape(self.shape)
self.functional.computeGradF(gap, grad)
return np.ravel(grad)
opt = scipy.optimize.minimize(
cost,
np.full(self.N, mean_gap),
jac=jac,
constraints=scipy.optimize.LinearConstraint(
np.full((1, self.N), 1 / self.N), mean_gap, mean_gap),
method='trust-constr',
tol=self.tolerance * mean_gap)
print(opt)
# Setup solved model
self.model['gap'][:] = opt.x.reshape(self.shape)
self.model.displacement[:] = self.model['gap'] + self.surface.reshape(
self.shape)
self.model.solveDirichlet()
self.model.traction[:] -= self.model.traction.min()
def main():
"""Run a simple Hertzian contact with penalty."""
N = 64
model = tm.Model(tm.model_type.basic_2d, [1, 1], [N] * 2)
surface = hertz_surface(model.system_size, model.shape, 1)
solver = ScipyContactSolver(model, surface, 1e-10, 1)
solver.solve(0.04)
plt.imshow(model.traction)
plt.colorbar()
plt.show()
if __name__ == "__main__":
main()
diff --git a/examples/statistics.py b/examples/statistics.py
index ca0ad5a..32c77ba 100644
--- a/examples/statistics.py
+++ b/examples/statistics.py
@@ -1,84 +1,84 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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 tamaas.utils import publications
from matplotlib.colors import LogNorm
tm.set_log_level(tm.LogLevel.info) # Show progression of solver
# Surface size
n = 512
# Surface generator
sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 1
# Spectrum
sg.spectrum = tm.Isopowerlaw2D()
# Parameters
sg.spectrum.q0 = 16
sg.spectrum.q1 = 16
sg.spectrum.q2 = 64
sg.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(np.clip(psd.real, 1e-10, np.inf), 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:
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:
grid.addPointData(uvw.DataArray(surface, range(2), 'surface'))
except ImportError:
print("uvw not installed, will not write VTK file")
pass
publications()
diff --git a/examples/stresses.py b/examples/stresses.py
index a557794..1114ee2 100644
--- a/examples/stresses.py
+++ b/examples/stresses.py
@@ -1,123 +1,123 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 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 argparse
import os
import numpy as np
import tamaas as tm
from tamaas.dumpers import H5Dumper as Dumper
from tamaas.dumpers._helper import hdf5toVTK
from tamaas.utils import publications
parser = argparse.ArgumentParser(
description="Hertzian tractios applied on elastic half-space")
parser.add_argument("radius", type=float, help="Radius of sphere")
parser.add_argument("load", type=float, help="Applied normal force")
parser.add_argument("name", help="Output file name")
parser.add_argument("--plots", help='Show surface normal stress',
action="store_true")
args = parser.parse_args()
# Definition of modeled domain
model_type = tm.model_type.volume_2d
discretization = [32, 127, 127]
system_size = [0.25, 1., 1.]
# Material contants
E = 1. # Young's modulus
nu = 0.3 # Poisson's ratio
E_star = E/(1-nu**2) # Hertz modulus
# Creation of model
model = tm.ModelFactory.createModel(model_type,
system_size,
discretization)
model.E = E
model.nu = nu
# Setup for integral operators
residual = tm.ModelFactory.createResidual(model, 0, 0)
# Coordinates
x = np.linspace(0, system_size[1], discretization[1], endpoint=False)
y = np.linspace(0, system_size[2], discretization[2], endpoint=False)
x, y = np.meshgrid(x, y, indexing='ij')
center = [0.5, 0.5]
r = np.sqrt((x-center[0])**2 + (y-center[1])**2)
# Span of local data
local_size = model.boundary_shape
local_offset = tm.mpi.local_offset(r.shape)
local_slice = np.s_[local_offset:local_offset+local_size[0], :]
# Sphere
R = args.radius
P = args.load
# Contact area
a = (3*P*R/(4*E_star))**(1/3)
p_0 = 3 * P / (2 * np.pi * a**2)
# Pressure definition
traction = model.traction
hertz_traction = np.zeros(discretization[1:])
hertz_traction[r < a] = p_0 * np.sqrt(1 - (r[r < a] / a)**2)
traction[..., 2] = hertz_traction[local_slice]
# Array references
displacement = model.displacement
stress = model['stress']
gradient = np.zeros_like(stress)
# Getting integral operators
boussinesq = model.operators["boussinesq"]
boussinesq_gradient = model.operators["boussinesq_gradient"]
# Applying operators
boussinesq(traction, displacement)
boussinesq_gradient(traction, gradient)
model.operators["hooke"](gradient, stress)
# Dumper
dumper_helper = Dumper(args.name, 'stress')
model.addDumper(dumper_helper)
model.dump()
# Converting HDF dump to VTK
with tm.mpi.sequential():
if tm.mpi.rank() == 0:
hdf5toVTK(os.path.join('hdf5', f'{args.name}_0000.h5'), args.name)
if args.plots:
import matplotlib.pyplot as plt # noqa
fig, ax = plt.subplots(1, 2)
fig.suptitle("Rank {}".format(tm.mpi.rank()))
ax[0].set_title("Traction")
ax[1].set_title("Normal Stress")
ax[0].imshow(traction[..., 2])
ax[1].imshow(stress[0, ..., 2])
fig.tight_layout()
plt.show()
publications()
diff --git a/python/SConscript b/python/SConscript
index f6467b9..506aeb9 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,177 +1,177 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 SCons.Script import Import, Split, Copy, Dir
Import('main_env')
# Pybind11 wrapper
env_pybind = main_env.Clone(SHLIBPREFIX='')
# Remove pedantic warnings
cxx_flags = env_pybind['CXXFLAGS']
try:
del cxx_flags[cxx_flags.index('-pedantic')]
except ValueError:
pass
env_pybind.Tool(pybind11)
pybind_sources = Split("""
tamaas_module.cpp
wrap/core.cpp
wrap/percolation.cpp
wrap/surface.cpp
wrap/model.cpp
wrap/solvers.cpp
wrap/compute.cpp
wrap/mpi.cpp
wrap/test_features.cpp
""")
# Setting paths to find libTamaas
env_pybind.AppendUnique(LIBPATH=['../src'])
# Link against a static libTamaas
if env_pybind['build_static_lib']:
env_pybind.PrependUnique(LIBS=['Tamaas']) # keep other libs for link
env_pybind['RPATH'] = "" # no need for rpath w/ static lib
# Link against a dynamic libTamaas
else:
env_pybind.AppendUnique(RPATH=[
"'$$$$ORIGIN/../../src'", # path to lib in build_dir
"'$$$$ORIGIN/../../..'", # path to lib in install prefix
])
env_pybind['LIBS'] = ['Tamaas'] # discard other libs for link
# Building the pybind library
tamaas_wrap = env_pybind.Pybind11Module(
target='tamaas/_tamaas',
source=pybind_sources,
)
# For some reason link happens too early
Import('libTamaas')
env_pybind.Depends(tamaas_wrap, libTamaas)
# Copying the __init__.py file with extra python classes
copy_env = env_pybind.Clone()
# Copying additional python files
python_files = """
__main__.py
compute.py
utils.py
dumpers/__init__.py
dumpers/_helper.py
nonlinear_solvers/__init__.py
""".split()
targets = [tamaas_wrap]
targets += [
copy_env.Command(copy_env.File(f, 'tamaas'),
copy_env.File(f, '#python/tamaas'),
Copy("$TARGET", "$SOURCE"))
for f in python_files
]
dist_files = """
MANIFEST.in
pypi.md
setup.py
""".split()
# pyproject.toml causes issues with develop mode
dist_files.append("pyproject.toml")
targets += [
copy_env.Command(copy_env.File(f, ''),
copy_env.File(f, '#python'),
Copy("$TARGET", "$SOURCE"))
for f in dist_files
]
subst_env = env_pybind.Clone(
SUBST_DICT={
'@version@': '$version',
'@authors@': str(copy_env['authors']),
'@author_list@': ', '.join(copy_env['authors']),
'@email@': '$email',
'@description@': '$description',
'@copyright@': '$copyright',
'@maintainer@': '$maintainer',
'@license@': '$license',
}
)
subst_env.Tool('textfile')
targets.append(subst_env.Substfile('tamaas/__init__.py.in'))
targets.append(subst_env.Substfile('setup.cfg.in'))
# Defining alias for python builds
main_env.Alias('build-python', targets)
# Checking if we can use pip to install (more convenient for end-user)
install_env = main_env.Clone()
conf = Configure(install_env,
custom_tests={'CheckPythonModule': CheckPythonModule})
has_pip = conf.CheckPythonModule('pip')
install_env = conf.Finish()
# Current build directory
install_env['PYDIR'] = Dir('.')
# Setting command line for installation
if has_pip:
install_env['PYINSTALLCOM'] = '${py_exec} -m pip install -U $PYOPTIONS .'
install_env['PYDEVELOPCOM'] = \
'${py_exec} -m pip install $PYOPTIONS -e .[all]'
else:
install_env['PYINSTALLCOM'] = '${py_exec} setup.py install $PYOPTIONS'
install_env['PYDEVELOPCOM'] = '${py_exec} setup.py develop $PYOPTIONS'
install_env['py_version'] = get_python_version(install_env)
install_env.PrependENVPath(
'PYTHONPATH',
install_env.subst('${prefix}/lib/python${py_version}/site-packages'))
# Specify install target
PYOPTIONS = ['${"" if verbose else "-q"}']
python_install = install_env.Command(
'.python_install_phony',
targets, install_env['PYINSTALLCOM'],
PYOPTIONS=['--prefix', '${prefix}'] + PYOPTIONS,
chdir=install_env['PYDIR'])
python_install_dev = install_env.Command(
'.python_install_local_phony',
targets, install_env['PYDEVELOPCOM'],
# Temporary fix for https://github.com/pypa/pip/issues/7953
PYOPTIONS=['--prefix=$$HOME/.local/'] + PYOPTIONS,
chdir=install_env['PYDIR'])
# Defining aliases
main_env.Alias('install-python', python_install)
main_env.Alias('dev', python_install_dev)
diff --git a/python/cast.hh b/python/cast.hh
index 7ab3f19..2cd30f4 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,245 +1,245 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef CAST_HH
#define CAST_HH
/* -------------------------------------------------------------------------- */
#include "field_container.hh"
#include "grid.hh"
#include "grid_base.hh"
#include "model_type.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace pybind11 {
// Format descriptor necessary for correct wrap of tamaas complex type
template
struct format_descriptor<
tamaas::complex, detail::enable_if_t::value>> {
static constexpr const char c = format_descriptor::c;
static constexpr const char value[3] = {'Z', c, '\0'};
static std::string format() { return std::string(value); }
};
#ifndef PYBIND11_CPP17
template
constexpr const char format_descriptor<
tamaas::complex,
detail::enable_if_t::value>>::value[3];
#endif
namespace detail {
// declare tamaas complex as a complex type for pybind11
template
struct is_complex> : std::true_type {};
template
struct is_fmt_numeric,
detail::enable_if_t::value>>
: std::true_type {
static constexpr int index = is_fmt_numeric::index + 3;
};
static inline handle policy_switch(return_value_policy policy, handle parent) {
switch (policy) {
case return_value_policy::copy:
case return_value_policy::move:
return handle();
case return_value_policy::automatic_reference: // happens in python-derived
// classes
case return_value_policy::reference:
return none();
case return_value_policy::reference_internal:
return parent;
default:
TAMAAS_EXCEPTION("Policy is not handled");
}
}
template
handle grid_to_python(const tamaas::Grid& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes(dim);
std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
if (grid.getNbComponents() != 1)
sizes.push_back(grid.getNbComponents());
return array(sizes, grid.getInternalData(), parent).release();
}
template
handle grid_to_python(const tamaas::GridBase& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes = {grid.dataSize()};
return array(sizes, grid.getInternalData(), parent).release();
}
/**
* Type caster for grid classes
* inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
* and pybind11/eigen.h
*/
template class G, typename T,
tamaas::UInt dim>
struct type_caster> {
using type = G;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridWrap"));
/**
* Conversion part 1 (Python->C++): convert a PyObject into a grid
* instance or return false upon failure. The second argument
* indicates whether implicit conversions should be applied.
*/
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf)
return false;
value.move(tamaas::wrap::GridNumpy>(buf));
return true;
}
/**
* Conversion part 2 (C++ -> Python): convert a grid instance into
* a Python object. The second and third arguments are used to
* indicate the return value policy and parent object (for
* ``return_value_policy::reference_internal``) and are generally
* ignored by implicit casters.
*/
static handle cast(const type& src, return_value_policy policy,
handle parent) {
return grid_to_python(src, policy, parent);
}
};
/**
* Type caster for GridBase classes
*/
template
struct type_caster> {
using type = tamaas::GridBase;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridBaseWrap"));
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf)
return false;
value.move(tamaas::wrap::GridBaseNumpy(buf));
return true;
}
static handle cast(const type& src, return_value_policy policy,
handle parent) {
using namespace ::tamaas::detail;
return tuple_dispatch_with_default(
[&](auto&& _) {
constexpr auto dim = std::decay_t::value;
const auto* conv = dynamic_cast*>(&src);
if (conv)
return grid_to_python(*conv, policy, parent);
return grid_to_python(src, policy, parent);
},
[&](auto&&) { return grid_to_python(src, policy, parent); },
src.getDimension());
}
};
/**
* Type caster for grid variant
*/
template <>
struct type_caster {
using type = typename tamaas::FieldContainer::Value;
using array_type = array;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridVariant"));
bool load(handle /*src*/, bool /*convert*/) {
return false;
// if (!array_type::check_(src) or !convert)
// return false;
// auto buf = array_type::ensure(src);
// if (!buf)
// return false;
// value.move(tamaas::wrap::GridBaseNumpy(buf));
// return true;
}
static handle cast(const type& src, return_value_policy policy,
handle parent) {
using namespace ::tamaas::detail;
return boost::apply_visitor(
[policy, parent](auto&& grid_ptr) {
return tuple_dispatch_with_default(
[&](auto&& _) {
using T =
typename std::decay_t::value_type;
constexpr auto dim = std::decay_t::value;
const auto* conv =
dynamic_cast*>(grid_ptr.get());
if (conv)
return grid_to_python(*conv, policy, parent);
return grid_to_python(*grid_ptr, policy, parent);
},
[&](auto&&) {
return grid_to_python(*grid_ptr, policy, parent);
},
grid_ptr->getDimension());
},
src);
}
};
} // namespace detail
} // namespace pybind11
#endif // CAST_HH
diff --git a/python/numpy.hh b/python/numpy.hh
index e18b78f..80d3ff9 100644
--- a/python/numpy.hh
+++ b/python/numpy.hh
@@ -1,93 +1,93 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef NUMPY_HH
#define NUMPY_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
#include
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// Type alias for usual numpy array type
template
using numpy = py::array_t;
/// GridBase helper class wrapping numpy array
template
class GridBaseNumpy : public GridBase {
public:
GridBaseNumpy(numpy& buffer) : GridBase() {
this->nb_components = 1;
this->data.wrap(buffer.mutable_data(), buffer.size());
}
};
/// Grid helper class wrapping numpy array
template
class GridNumpy : public Parent {
public:
GridNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
const UInt ndim = buffer.ndim();
if (ndim > Parent::dimension + 1 || ndim < Parent::dimension)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected grid dimensions");
if (ndim == Parent::dimension + 1)
this->nb_components = array_shape[ndim - 1];
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrap(buffer.mutable_data(), this->computeSize());
}
};
/// Surface helper class wrapping numpy array
template
class SurfaceNumpy : public Parent {
public:
SurfaceNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
if (buffer.ndim() != 2)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected surface dimensions");
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrap(buffer.mutable_data(), this->computeSize());
}
};
} // namespace wrap
} // namespace tamaas
#endif // NUMPY_HH
diff --git a/python/tamaas/__init__.py.in b/python/tamaas/__init__.py.in
index 9180365..af2fa62 100644
--- a/python/tamaas/__init__.py.in
+++ b/python/tamaas/__init__.py.in
@@ -1,59 +1,59 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""
@description@
See __author__, __license__, __copyright__ for extra information about Tamaas.
- User documentation: https://tamaas.readthedocs.io
- Bug Tracker: https://gitlab.com/tamaas/tamaas/-/issues
- Source Code: https://gitlab.com/tamaas/tamaas
"""
__author__ = @authors@
__copyright__ = "@copyright@"
__license__ = "@license@"
__maintainer__ = "@maintainer@"
__email__ = "@email@"
try:
from ._tamaas import model_type, TamaasInfo
from ._tamaas import _type_traits as __tt
type_traits = {
model_type.basic_1d: __tt.basic_1d,
model_type.basic_2d: __tt.basic_2d,
model_type.surface_1d: __tt.surface_1d,
model_type.surface_2d: __tt.surface_2d,
model_type.volume_1d: __tt.volume_1d,
model_type.volume_2d: __tt.volume_2d,
}
del __tt
from ._tamaas import * # noqa
__version__ = TamaasInfo.version
except ImportError as e:
print("Error trying to import _tamaas:\n{}".format(e))
raise e
diff --git a/python/tamaas/__main__.py b/python/tamaas/__main__.py
index ca81bab..879fa74 100755
--- a/python/tamaas/__main__.py
+++ b/python/tamaas/__main__.py
@@ -1,257 +1,257 @@
#!/usr/bin/env python3
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""Module entry point."""
import sys
import io
import time
import argparse
import tamaas as tm
import numpy as np
__author__ = "Lucas Frérot"
__copyright__ = (
- "Copyright (©) 2019-2022, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2023, EPFL (École Polytechnique Fédérale de Lausanne),"
"\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
- "\nCopyright (©) 2020-2022 Lucas Frérot"
+ "\nCopyright (©) 2020-2023 Lucas Frérot"
)
__license__ = "SPDX-License-Identifier: AGPL-3.0-or-later"
def load_stream(stream):
"""
Load numpy from binary stream (allows piping).
Code from
https://gist.github.com/CMCDragonkai/3c99fd4aabc8278b9e17f50494fcc30a
"""
np_magic = stream.read(6)
# use the sys.stdin.buffer to read binary data
np_data = stream.read()
# read it all into an io.BytesIO object
return io.BytesIO(np_magic + np_data)
def print_version():
"""Print version information."""
print(
f"""\
Tamaas {tm.__version__}
{tm.__copyright__}
Authors: {', '.join(tm.__author__)}
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Tamaas is the fruit of a research effort. Consider citing 10.21105/joss.02121
and the relevant references therein. Use the function tamaas.utils.publications
at the end of your python scripts for an exhaustive publication list."""
)
def surface(args):
"""Generate a surface."""
if args.generator == "random_phase":
generator = tm.SurfaceGeneratorRandomPhase2D(args.sizes)
elif args.generator == "filter":
generator = tm.SurfaceGeneratorFilter2D(args.sizes)
else:
raise ValueError("Unknown generator method {}".format(args.generator))
generator.spectrum = tm.Isopowerlaw2D()
generator.spectrum.q0 = args.cutoffs[0]
generator.spectrum.q1 = args.cutoffs[1]
generator.spectrum.q2 = args.cutoffs[2]
generator.spectrum.hurst = args.hurst
generator.random_seed = args.seed
surface = (
generator.buildSurface() / generator.spectrum.rmsSlopes() * args.rms
)
output = args.output if args.output is not None else sys.stdout
params = {
"q0": generator.spectrum.q0,
"q1": generator.spectrum.q1,
"q2": generator.spectrum.q2,
"hurst": generator.spectrum.hurst,
"random_seed": generator.random_seed,
"rms_heights": args.rms,
"generator": args.generator,
}
try:
np.savetxt(output, surface, header=str(params))
except BrokenPipeError:
pass
def contact(args):
"""Solve a contact problem."""
from tamaas.dumpers import NumpyDumper
tm.set_log_level(tm.LogLevel.error)
if not args.input:
input = sys.stdin
else:
input = args.input
surface = np.loadtxt(input)
discretization = surface.shape
system_size = [1.0, 1.0]
model = tm.ModelFactory.createModel(
tm.model_type.basic_2d, system_size, discretization
)
solver = tm.PolonskyKeerRey(model, surface, args.tol)
solver.solve(args.load)
dumper = NumpyDumper("numpy", "traction", "displacement")
dumper.dump_to_file(sys.stdout.buffer, model)
def plot(args):
"""Plot displacement and pressure maps."""
try:
import matplotlib.pyplot as plt
except ImportError:
print(
"Please install matplotlib to use the 'plot' command",
file=sys.stderr,
)
sys.exit(1)
fig, (ax_traction, ax_displacement) = plt.subplots(1, 2)
ax_traction.set_title("Traction")
ax_displacement.set_title("Displacement")
with load_stream(sys.stdin.buffer) as f_np:
data = np.load(f_np)
ax_traction.imshow(data["traction"])
ax_displacement.imshow(data["displacement"])
fig.set_size_inches(10, 6)
fig.tight_layout()
plt.show()
def main():
"""Main function."""
parser = argparse.ArgumentParser(
prog="tamaas",
description=(
"The tamaas command is a simple utility for surface"
" generation, contact computation and"
" plotting of contact solutions"
),
)
parser.add_argument(
"--version", action="store_true", help="print version info"
)
subs = parser.add_subparsers(
title="commands", description="utility commands"
)
# Arguments for surface command
parser_surface = subs.add_parser(
"surface", description="Generate a self-affine rough surface"
)
parser_surface.add_argument(
"--cutoffs",
"-K",
nargs=3,
type=int,
help="Long, rolloff, short wavelength cutoffs",
metavar=("k_l", "k_r", "k_s"),
required=True,
)
parser_surface.add_argument(
"--sizes",
nargs=2,
type=int,
help="Number of points",
metavar=("nx", "ny"),
required=True,
)
parser_surface.add_argument(
"--hurst", "-H", type=float, help="Hurst exponent", required=True
)
parser_surface.add_argument(
"--rms", type=float, help="Root-mean-square of slopes", default=1.0
)
parser_surface.add_argument(
"--seed", type=int, help="Random seed", default=int(time.time())
)
parser_surface.add_argument(
"--generator",
help="Generation method",
choices=("random_phase", "filter"),
default="random_phase",
)
parser_surface.add_argument(
"--output", "-o", help="Output file name (compressed if .gz)"
)
parser_surface.set_defaults(func=surface)
# Arguments for contact command
parser_contact = subs.add_parser(
"contact",
description="Compute the elastic contact solution with a given surface",
)
parser_contact.add_argument(
"--input", "-i", help="Rough surface file (default stdin)"
)
parser_contact.add_argument(
"--tol", type=float, default=1e-12, help="Solver tolerance"
)
parser_contact.add_argument(
"load", type=float, help="Applied average pressure"
)
parser_contact.set_defaults(func=contact)
# Arguments for plot command
parser_plot = subs.add_parser("plot", description="Plot contact solution")
parser_plot.set_defaults(func=plot)
args = parser.parse_args()
if args.version:
print_version()
return
try:
args.func(args)
except AttributeError:
parser.print_usage()
if __name__ == "__main__":
main()
diff --git a/python/tamaas/compute.py b/python/tamaas/compute.py
index b6630b5..6739fa4 100644
--- a/python/tamaas/compute.py
+++ b/python/tamaas/compute.py
@@ -1,21 +1,21 @@
# -*- coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 ._tamaas.compute import * # noqa
diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py
index 6e01423..b76c2a3 100644
--- a/python/tamaas/dumpers/__init__.py
+++ b/python/tamaas/dumpers/__init__.py
@@ -1,581 +1,581 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""Dumpers for the class :py:class:`Model `."""
from pathlib import Path
from os import PathLike
import glob
import json
import io
import typing as ts
from collections.abc import Collection
import numpy as np
from .. import (
ModelDumper,
model_type,
mpi,
type_traits,
ModelFactory,
Model,
__version__,
)
from ._helper import (
step_dump,
directory_dump,
local_slice,
_is_surface_field,
_basic_types,
file_handler,
)
__all__ = [
"JSONDumper",
"FieldDumper",
"NumpyDumper",
]
FileType = ts.Union[str, PathLike, io.TextIOBase]
NameType = ts.Union[str, PathLike]
_reverse_trait_map = {
'model_type.' + t.__name__: mtype
for mtype, t in type_traits.items()
}
def _get_attributes(model: Model):
"""Get model attributes."""
return {
'model_type': str(model.type),
'system_size': model.system_size,
'discretization': model.global_shape,
'boundary_fields': model.boundary_fields,
'program': f"Tamaas {__version__}, DOI:10.21105/joss.02121",
}
def _create_model(attrs: ts.MutableMapping):
"""Create a model from attribute dictionary."""
mtype = _reverse_trait_map[attrs['model_type']]
# netcdf4 converts 1-lists attributes to numbers
for attr in ['system_size', 'discretization']:
if not isinstance(attrs[attr], Collection):
attrs[attr] = [attrs[attr]]
return ModelFactory.createModel(mtype, attrs['system_size'],
attrs['discretization'])
class MPIIncompatibilityError(RuntimeError):
"""Raised when code is not meant to be executed in MPI environment."""
class ModelError(ValueError):
"""Raised when unexpected model is passed to a dumper with a state."""
class ComponentsError(ValueError):
"""Raised when an unexpected number of components is encountred."""
class _ModelJSONEncoder(json.JSONEncoder):
"""Encode a model to JSON."""
def default(self, obj):
"""Encode model."""
if isinstance(obj, Model):
model = obj
attrs = _get_attributes(model)
model_dict = {
'attrs': attrs,
'fields': {},
'operators': [],
}
for field in model:
model_dict['fields'][field] = model[field].tolist()
for op in model.operators:
model_dict['operators'].append(op)
return model_dict
return json.JSONEncoder.default(self, obj)
class JSONDumper(ModelDumper):
"""Dumper to JSON."""
def __init__(self, file_descriptor: FileType):
"""Construct with file handle."""
super(JSONDumper, self).__init__()
self.fd = file_descriptor
@file_handler('w')
def _dump_to_file(self, fd: FileType, model: Model):
json.dump(model, fd, cls=_ModelJSONEncoder)
def dump(self, model: Model):
"""Dump model."""
self._dump_to_file(self.fd, model)
@classmethod
@file_handler('r')
def read(cls, fd: FileType):
"""Read model from file."""
properties = json.load(fd)
model = _create_model(properties['attrs'])
for name, field in properties['fields'].items():
v = np.asarray(field)
if model.type in _basic_types:
v = v.reshape(list(v.shape) + [1])
model[name] = v
return model
class FieldDumper(ModelDumper):
"""Abstract dumper for python classes using fields."""
postfix = ""
extension = ""
name_format = "{basename}{postfix}.{extension}"
def __init__(self, basename: NameType, *fields, **kwargs):
"""Construct with desired fields."""
super(FieldDumper, self).__init__()
self.basename = basename
self.fields: ts.List[str] = list(fields)
self.all_fields: bool = kwargs.get('all_fields', False)
def add_field(self, field: str):
"""Add another field to the dump."""
if field not in self.fields:
self.fields.append(field)
def _dump_to_file(self, file_descriptor: FileType, model: Model):
"""Dump to a file (path-like or file handle)."""
raise NotImplementedError()
def get_fields(self, model: Model):
"""Get the desired fields."""
if not self.all_fields:
requested_fields = self.fields
else:
requested_fields = list(model)
return {field: model[field] for field in requested_fields}
def dump(self, model: Model):
"""Dump model."""
self._dump_to_file(self.file_path, model)
@classmethod
def read(cls, file_descriptor: FileType):
"""Read model from file."""
raise NotImplementedError(
f'read() method not implemented in {cls.__name__}')
@classmethod
def read_sequence(cls, glob_pattern):
"""Read models from a file sequence."""
return map(cls.read, glob.iglob(glob_pattern))
@property
def file_path(self):
"""Get the default filename."""
return self.name_format.format(basename=self.basename,
postfix=self.postfix,
extension=self.extension)
@directory_dump('numpys')
@step_dump
class NumpyDumper(FieldDumper):
"""Dumper to compressed numpy files."""
extension = 'npz'
def _dump_to_file(self, file_descriptor: FileType, model: Model):
"""Save to compressed multi-field Numpy format."""
if mpi.size() > 1:
raise MPIIncompatibilityError("NumpyDumper does not function "
"at all in parallel")
np.savez_compressed(file_descriptor,
attrs=json.dumps(_get_attributes(model)),
**self.get_fields(model))
@classmethod
def read(cls, file_descriptor: FileType):
"""Create model from Numpy file."""
data = np.load(file_descriptor, mmap_mode='r')
model = _create_model(json.loads(str(data['attrs'])))
for k, v in filter(lambda k: k[0] != 'attrs', data.items()):
if model.type in _basic_types:
v = v.reshape(list(v.shape) + [1])
model[k] = v
return model
try:
import h5py
__all__.append("H5Dumper")
@directory_dump('hdf5')
@step_dump
class H5Dumper(FieldDumper):
"""Dumper to HDF5 file format."""
extension = 'h5'
@staticmethod
def _hdf5_args():
if mpi.size() > 1:
from mpi4py import MPI # noqa
mpi_args = dict(driver='mpio', comm=MPI.COMM_WORLD)
comp_args = {} # compression does not work in parallel
else:
mpi_args = {}
comp_args = dict(compression='gzip', compression_opts=7)
return mpi_args, comp_args
def _dump_to_file(self, file_descriptor: FileType, model: Model):
"""Save to HDF5 with metadata about the model."""
# Setup for MPI
if not h5py.get_config().mpi and mpi.size() > 1:
raise MPIIncompatibilityError("HDF5 does not have MPI support")
mpi_args, comp_args = self._hdf5_args()
with h5py.File(file_descriptor, 'w', **mpi_args) as handle:
# Writing data
for name, field in self.get_fields(model).items():
shape = list(field.shape)
if mpi.size() > 1:
xdim = 0 if _is_surface_field(field, model) else 1
shape[xdim] = mpi_args['comm'].allreduce(shape[xdim])
dset = handle.create_dataset(name, shape, field.dtype,
**comp_args)
s = local_slice(field, model, name in model.boundary_fields)
dset[s] = field
# Writing metadata
for name, attr in _get_attributes(model).items():
handle.attrs[name] = attr
@classmethod
def read(cls, file_descriptor: FileType):
"""Create model from HDF5 file."""
mpi_args, _ = cls._hdf5_args()
with h5py.File(file_descriptor, 'r', **mpi_args) as handle:
model = _create_model(handle.attrs)
for k, v in handle.items():
v = np.asanyarray(v)
if model.type in _basic_types:
v = v.reshape(list(v.shape) + [1])
surface_field = \
k in handle.attrs.get('boundary_fields', {}) \
or _is_surface_field(v, model)
s = local_slice(v, model, surface_field)
if (surface_field and v.ndim == len(model.boundary_shape)) \
or (not surface_field and v.ndim == len(model.shape)):
s = s + (np.newaxis, )
model[k] = v[s].copy()
return model
except ImportError:
pass
try:
import uvw # noqa
__all__ += [
"UVWDumper",
"UVWGroupDumper",
]
@directory_dump('paraview')
@step_dump
class UVWDumper(FieldDumper):
"""Dumper to VTK files for elasto-plastic calculations."""
extension = 'vtr'
def _dump_to_file(self, file_descriptor: FileType, model: Model):
"""Dump displacements, plastic deformations and stresses."""
if mpi.size() > 1:
raise MPIIncompatibilityError("UVWDumper does not function "
"properly in parallel")
bdim = len(model.boundary_shape)
# Local MPI size
lsize = model.shape
gsize = mpi.global_shape(model.boundary_shape)
gshape = gsize
if len(lsize) > bdim:
gshape = [model.shape[0]] + gshape
# Space coordinates
coordinates = [
np.linspace(0, L, N, endpoint=False)
for L, N in zip(model.system_size, gshape)
]
# If model has subsurfce domain, z-coordinate is always first
dimension_indices = np.arange(bdim)
if len(lsize) > bdim:
dimension_indices += 1
dimension_indices = np.concatenate((dimension_indices, [0]))
coordinates[0] = \
np.linspace(0, model.system_size[0], gshape[0])
offset = np.zeros_like(dimension_indices)
offset[0] = mpi.local_offset(gsize)
rectgrid = uvw.RectilinearGrid if mpi.size() == 1 \
else uvw.parallel.PRectilinearGrid
# Creating rectilinear grid with correct order for components
coordlist = [
coordinates[i][o:o + lsize[i]]
for i, o in zip(dimension_indices, offset)
]
grid = rectgrid(
file_descriptor,
coordlist,
compression=True,
offsets=offset,
)
fields = self.get_fields(model).items()
# Iterator over fields we want to dump on system geometry
if model.type in {model_type.volume_1d, model_type.volume_2d}:
fields_it = filter(lambda t: not t[0] in model.boundary_fields,
fields)
else:
fields_it = iter(fields)
for name, field in fields_it:
array = uvw.DataArray(field, dimension_indices, name)
grid.addPointData(array)
grid.write()
@directory_dump('paraview')
class UVWGroupDumper(FieldDumper):
"""Dumper to ParaViewData files."""
extension = 'pvd'
def __init__(self, basename: NameType, *fields, **kwargs):
"""Construct with desired fields."""
super(UVWGroupDumper, self).__init__(basename, *fields, **kwargs)
subdir = Path('paraview') / f'{basename}-VTR'
subdir.mkdir(parents=True, exist_ok=True)
self.uvw_dumper = UVWDumper(
Path(f'{basename}-VTR') / basename, *fields, **kwargs)
self.group = uvw.ParaViewData(self.file_path, compression=True)
def _dump_to_file(self, file_descriptor, model):
self.group.addFile(
self.uvw_dumper.file_path.replace('paraview/', ''),
timestep=self.uvw_dumper.count,
)
self.group.write()
self.uvw_dumper.dump(model)
except ImportError:
pass
try:
from netCDF4 import Dataset
__all__.append("cls")
@directory_dump('netcdf')
class NetCDFDumper(FieldDumper):
"""Dumper to netCDF4 files."""
extension = "nc"
time_dim = 'frame'
format = 'NETCDF4'
def _file_setup(self, grp, model: Model):
grp.createDimension(self.time_dim, None)
# Attibutes
for k, v in _get_attributes(model).items():
grp.setncattr(k, v)
# Local dimensions
voigt_dim = type_traits[model.type].voigt
components = type_traits[model.type].components
self._vec = grp.createDimension('spatial', components)
self._tens = grp.createDimension('Voigt', voigt_dim)
self.model_info = model.global_shape, model.type
global_boundary_shape = mpi.global_shape(model.boundary_shape)
# Create boundary dimensions
for label, size, length in zip("xy", global_boundary_shape,
model.boundary_system_size):
grp.createDimension(label, size)
coord = grp.createVariable(label, 'f8', (label, ))
coord[:] = np.linspace(0, length, size, endpoint=False)
self._create_variables(grp, model,
lambda f: _is_surface_field(f[1], model),
global_boundary_shape, "xy")
# Create volume dimension
if model.type in {model_type.volume_1d, model_type.volume_2d}:
size = model.shape[0]
grp.createDimension("z", size)
coord = grp.createVariable("z", 'f8', ("z", ))
coord[:] = np.linspace(0, model.system_size[0], size)
self._create_variables(
grp, model, lambda f: not _is_surface_field(f[1], model),
model.global_shape, "zxy")
self.has_setup = True
def _set_collective(self, rootgrp):
if mpi.size() == 1:
return
for v in rootgrp.variables.values():
if self.time_dim in v.dimensions:
v.set_collective(True)
def _dump_to_file(self, file_descriptor: NameType, model: Model):
mode = 'a' if Path(file_descriptor).is_file() \
and getattr(self, 'has_setup', False) else 'w'
try:
with Dataset(file_descriptor,
mode,
format=self.format,
parallel=mpi.size() > 1) as rootgrp:
if rootgrp.dimensions == {}:
self._file_setup(rootgrp, model)
self._set_collective(rootgrp)
if self.model_info != (model.global_shape, model.type):
raise ModelError(f"Unexpected model {mode}")
self._dump_generic(rootgrp, model)
except ValueError:
raise MPIIncompatibilityError("NetCDF4 has no MPI support")
def _create_variables(self, grp, model, predicate, shape, dimensions):
field_dim = len(shape)
fields = list(filter(predicate, self.get_fields(model).items()))
dim_labels = list(dimensions[:field_dim])
for label, data in fields:
local_dim = []
# If we have an extra component
if data.ndim > field_dim:
if data.shape[-1] == self._tens.size:
local_dim = [self._tens.name]
elif data.shape[-1] == self._vec.size:
local_dim = [self._vec.name]
else:
raise ComponentsError(
f"{label} has unexpected number of components "
f"({data.shape[-1]})")
# Downcasting in case of 128 bit float
dtype = data.dtype if data.dtype.str[1:] != 'f16' else 'f8'
grp.createVariable(label,
dtype,
[self.time_dim] + dim_labels + local_dim,
zlib=mpi.size() == 0)
def _dump_generic(self, grp, model):
fields = self.get_fields(model).items()
new_frame = len(grp.dimensions[self.time_dim])
for label, data in fields:
var = grp[label]
slice_in_global = (new_frame, ) + local_slice(data, model)
var[slice_in_global] = np.array(data, dtype=var.dtype)
@classmethod
def _open_read(cls, fd):
return Dataset(fd, 'r', format=cls.format, parallel=mpi.size() > 1)
@staticmethod
def _create_model(rootgrp):
attrs = {k: rootgrp.getncattr(k) for k in rootgrp.ncattrs()}
return _create_model(attrs)
@staticmethod
def _set_model_fields(rootgrp, model, frame):
dims = rootgrp.dimensions.keys()
for k, v in filter(lambda k: k[0] not in dims,
rootgrp.variables.items()):
v = v[frame, :]
if model.type in _basic_types:
v = np.asarray(v).reshape(list(v.shape) + [1])
model[k] = v[local_slice(v, model)].copy()
@classmethod
def read(cls, file_descriptor: NameType):
"""Create model with last frame."""
with cls._open_read(file_descriptor) as rootgrp:
model = cls._create_model(rootgrp)
cls._set_model_fields(rootgrp, model, -1)
return model
@classmethod
def read_sequence(cls, file_descriptor: NameType):
with cls._open_read(file_descriptor) as rootgrp:
model = cls._create_model(rootgrp)
for frame in range(len(rootgrp.dimensions[cls.time_dim])):
cls._set_model_fields(rootgrp, model, frame)
yield model
except ImportError:
pass
diff --git a/python/tamaas/dumpers/_helper.py b/python/tamaas/dumpers/_helper.py
index 692c2c7..c36e499 100644
--- a/python/tamaas/dumpers/_helper.py
+++ b/python/tamaas/dumpers/_helper.py
@@ -1,174 +1,174 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""
Helper functions for dumpers
"""
from os import PathLike
from functools import wraps
from pathlib import Path
import io
import numpy as np
from .. import model_type, type_traits, mpi
__all__ = ["step_dump", "directory_dump", "local_slice"]
_basic_types = [t for t, trait in type_traits.items() if trait.components == 1]
def _is_surface_field(field, model):
def _to_global(shape):
if len(shape) == len(model.boundary_shape) + 1:
return mpi.global_shape(model.boundary_shape) + [shape[-1]]
return mpi.global_shape(shape)
def _compare_shape(a, b):
return a == b
b_shapes = [list(model[name].shape) for name in model.boundary_fields]
shape = list(field.shape)
return any(_compare_shape(shape, s) for s in b_shapes) \
or any(_compare_shape(shape, _to_global(s)) for s in b_shapes)
def local_slice(field, model, surface_field=None):
n = model.shape
bn = model.boundary_shape
gshape = mpi.global_shape(bn)
offsets = np.zeros_like(gshape)
offsets[0] = mpi.local_offset(gshape)
if not surface_field:
surface_field = _is_surface_field(field, model)
if not surface_field and len(n) > len(bn):
gshape = [n[0]] + gshape
offsets = np.concatenate(([0], offsets))
shape = bn if surface_field else n
if len(field.shape) > len(shape):
shape += field.shape[len(shape):]
def sgen(offset, size):
return slice(offset, offset + size, None)
def sgen_basic(offset, size):
return slice(offset, offset + size)
slice_gen = sgen_basic if model_type in _basic_types else sgen
return tuple(map(slice_gen, offsets, shape))
def step_dump(cls):
"""
Decorator for dumper with counter for steps
"""
orig_init = cls.__init__
orig_dump = cls.dump
@wraps(cls.__init__)
def __init__(obj, *args, **kwargs):
orig_init(obj, *args, **kwargs)
obj.count = 0
def postfix(obj):
return "_{:04d}".format(obj.count)
@wraps(cls.dump)
def dump(obj, *args, **kwargs):
orig_dump(obj, *args, **kwargs)
obj.count += 1
cls.__init__ = __init__
cls.dump = dump
cls.postfix = property(postfix)
return cls
def directory_dump(directory=""):
"Decorator for dumper in a directory"
directory = Path(directory)
def actual_decorator(cls):
orig_dump = cls.dump
orig_filepath = cls.file_path.fget
orig_init = cls.__init__
@wraps(cls.__init__)
def init(obj, *args, **kwargs):
orig_init(obj, *args, **kwargs)
obj.mkdir = kwargs.get('mkdir', True)
@wraps(cls.dump)
def dump(obj, *args, **kwargs):
if mpi.rank() == 0 and getattr(obj, 'mkdir'):
directory.mkdir(parents=True, exist_ok=True)
orig_dump(obj, *args, **kwargs)
@wraps(cls.file_path.fget)
def file_path(obj):
if getattr(obj, 'mkdir'):
return str(directory / orig_filepath(obj))
return orig_filepath(obj)
cls.__init__ = init
cls.dump = dump
cls.file_path = property(file_path)
return cls
return actual_decorator
def hdf5toVTK(inpath: PathLike, outname: str):
"""Convert HDF5 dump of a model to VTK."""
from . import UVWDumper, H5Dumper # noqa
UVWDumper(outname, all_fields=True) << H5Dumper.read(inpath)
def netCDFtoParaview(inpath: PathLike, outname: str):
"""Convert NetCDF dump of model sequence to Paraview."""
from . import UVWGroupDumper, NetCDFDumper # noqa
dumper = UVWGroupDumper(outname, all_fields=True)
for model in NetCDFDumper.read_sequence(inpath):
dumper << model
def file_handler(mode):
"""Decorate a function to accept path-like or file handles."""
def _handler(func):
@wraps(func)
def _wrapped(self, fd, *args, **kwargs):
if isinstance(fd, (str, PathLike)):
with open(fd, mode) as fd:
return _wrapped(self, fd, *args, **kwargs)
elif isinstance(fd, io.TextIOBase):
return func(self, fd, *args, **kwargs)
raise TypeError(
f"Expected a path-like or file handle, got {type(fd)}")
return _wrapped
return _handler
diff --git a/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py
index d30b39a..1e28a69 100644
--- a/python/tamaas/nonlinear_solvers/__init__.py
+++ b/python/tamaas/nonlinear_solvers/__init__.py
@@ -1,171 +1,171 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""Nonlinear solvers for plasticity problems.
Solvers in this module use :py:mod:`scipy.optimize` to solve the implicit
non-linear equation for plastic deformations with fixed contact pressures.
"""
from functools import wraps
from scipy.optimize import newton_krylov, root
from scipy.optimize.nonlin import NoConvergence
from .. import EPSolver, Logger, LogLevel, mpi
from .._tamaas import _tolerance_manager
from .._tamaas import _DFSANESolver as DFSANECXXSolver
__all__ = ['NLNoConvergence',
'DFSANESolver',
'DFSANECXXSolver',
'NewtonKrylovSolver',
'ToleranceManager']
class NLNoConvergence(Exception):
"""Convergence not reached exception."""
class ScipySolver(EPSolver):
"""Base class for solvers wrapping SciPy routines."""
def __init__(self, residual, model=None, callback=None):
"""Construct nonlinear solver with residual.
:param residual: plasticity residual object
:param model: Deprecated
:param callback: passed on to the Scipy solver
"""
super(ScipySolver, self).__init__(residual)
if mpi.size() > 1:
raise RuntimeError("Scipy solvers cannot be used with MPI; "
"DFSANECXXSolver can be used instead")
self.callback = callback
self._x = self.getStrainIncrement()
self._residual = self.getResidual()
self.options = {'ftol': 0, 'fatol': 1e-9}
def solve(self):
"""Solve R(δε) = 0 using a Scipy function."""
# For initial guess, compute the strain due to boundary tractions
# self._residual.computeResidual(self._x)
# self._x[...] = self._residual.getVector()
EPSolver.beforeSolve(self)
# Scipy root callback
def compute_residual(vec):
self._residual.computeResidual(vec)
return self._residual.getVector().copy()
# Solve
Logger().get(LogLevel.debug) << \
"Entering non-linear solve\n"
self._x[...] = self._scipy_solve(compute_residual)
Logger().get(LogLevel.debug) << \
"Non-linear solve returned"
# Computing displacements
self._residual.computeResidualDisplacement(self._x)
def reset(self):
"""Set solution vector to zero."""
self._x[...] = 0
def _scipy_solve(self, compute_residual):
"""Actually call the scipy solver.
:param compute_residual: function returning residual for given variable
"""
raise NotImplementedError()
class NewtonKrylovSolver(ScipySolver):
"""Solve using a finite-difference Newton-Krylov method."""
def _scipy_solve(self, compute_residual):
try:
return newton_krylov(compute_residual, self._x,
f_tol=self.tolerance,
verbose=True, callback=self.callback)
except NoConvergence:
raise NLNoConvergence("Newton-Krylov did not converge")
class DFSANESolver(ScipySolver):
"""Solve using a spectral residual jacobianless method.
See :doi:`10.1090/S0025-5718-06-01840-0` for details on method and
the relevant Scipy `documentation
`_ for
details on parameters.
"""
def _scipy_solve(self, compute_residual):
solution = root(compute_residual,
self._x,
method='df-sane',
options={'ftol': 0, 'fatol': self.tolerance},
callback=self.callback)
Logger().get(LogLevel.info) << \
"DF-SANE/Scipy: {} ({} iterations, {})".format(
solution.message,
solution.nit,
self.tolerance)
if not solution.success:
raise NLNoConvergence("DF-SANE/Scipy did not converge")
return solution.x.copy()
def ToleranceManager(start, end, rate):
"""Decorate solver to manage tolerance of non-linear solve step."""
def actual_decorator(cls):
orig_init = cls.__init__
orig_solve = cls.solve
orig_update_state = cls.updateState
@wraps(cls.__init__)
def __init__(obj, *args, **kwargs):
orig_init(obj, *args, **kwargs)
obj.setToleranceManager(_tolerance_manager(start, end, rate))
@wraps(cls.solve)
def new_solve(obj, *args, **kwargs):
ftol = obj.tolerance
ftol *= rate
obj.tolerance = max(ftol, end)
return orig_solve(obj, *args, **kwargs)
@wraps(cls.updateState)
def updateState(obj, *args, **kwargs):
obj.tolerance = start
return orig_update_state(obj, *args, **kwargs)
cls.__init__ = __init__
# cls.solve = new_solve
# cls.updateState = updateState
return cls
return actual_decorator
diff --git a/python/tamaas/utils.py b/python/tamaas/utils.py
index 8147d41..825a167 100644
--- a/python/tamaas/utils.py
+++ b/python/tamaas/utils.py
@@ -1,311 +1,311 @@
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-# Copyright (©) 2020-2022 Lucas Frérot
+# Copyright (©) 2020-2023 Lucas Frérot
#
# 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 .
"""Convenience utilities."""
import inspect
from collections import namedtuple
from contextlib import contextmanager
from typing import Iterable, Union
from functools import partial
from operator import contains
import numpy as np
from . import (
ContactSolver,
Model,
SurfaceGenerator1D,
SurfaceGenerator2D,
dtype,
set_log_level,
get_log_level,
LogLevel,
Logger,
)
__all__ = [
"log_context",
"publications",
"load_path",
"seeded_surfaces",
"hertz_surface",
]
class NoConvergenceError(RuntimeError):
"""Convergence not reached exception."""
@contextmanager
def log_context(log_level: LogLevel):
"""Context manager to easily control Tamaas' logging level."""
current = get_log_level()
set_log_level(log_level)
try:
yield
finally:
set_log_level(current)
def publications(format_str="{pub.citation}\n\t{pub.doi}"):
"""Print publications associated with objects in use."""
Publication = namedtuple("Publication", ["citation", "doi"])
joss = Publication(
citation=(
"Frérot, L., Anciaux, G., Rey, V., Pham-Ba, S. & Molinari, J.-F."
" Tamaas: a library for elastic-plastic contact of periodic rough"
" surfaces. Journal of Open Source Software 5, 2121 (2020)."),
doi="10.21105/joss.02121",
)
zenodo = Publication(
citation=(
"Frérot, L., Anciaux, G., Rey, V., Pham-Ba, S., "
"& Molinari, J.-F. Tamaas, a high-performance "
"library for periodic rough surface contact. Zenodo (2019)."),
doi="10.5281/zenodo.3479236",
)
_publications = {
k: v
for keys, v in [
(
(
"SurfaceGeneratorRandomPhase1D",
"SurfaceGeneratorRandomPhase2D",
),
[
Publication(
citation=(
"Wu, J.-J. Simulation of rough surfaces with FFT."
" Tribology International 33, 47–58 (2000)."),
doi="10.1016/S0301-679X(00)00016-5",
),
],
),
(
("SurfaceGeneratorFilter1D", "SurfaceGeneratorFilter2D"),
[
Publication(
citation=(
"Hu, Y. Z. & Tonder, K. Simulation of 3-D random"
" rough surface by 2-D digital filter and fourier"
" analysis. International Journal of Machine Tools"
" and Manufacture 32, 83–90 (1992)."),
doi="10.1016/0890-6955(92)90064-N",
),
],
),
(
("PolonskyKeerRey", ),
[
Publication(
citation=(
"Polonsky, I. A. & Keer, L. M. A numerical method"
" for solving rough contact problems based on the"
" multi-level multi-summation and conjugate"
" gradient techniques. Wear 231, 206–219 (1999)."),
doi="10.1016/S0043-1648(99)00113-1",
),
Publication(
citation=(
"Rey, V., Anciaux, G. & Molinari, J.-F. Normal"
" adhesive contact on rough surfaces: efficient"
" algorithm for FFT-based BEM resolution. Comput"
" Mech 1–13 (2017)."),
doi="10.1007/s00466-017-1392-5",
),
],
),
(
("DFSANESolver", "DFSANECXXSolver"),
[
Publication(
citation=(
"La Cruz, W., Martínez, J. & Raydan, M. Spectral"
" residual method without gradient information for"
" solving large-scale nonlinear systems of"
" equations. Math. Comp. 75, 1429–1448 (2006)."),
doi="10.1090/S0025-5718-06-01840-0",
),
],
),
(
("Residual", ),
[
Publication(
citation=("Frérot, L., Bonnet, M., Molinari, J.-F. &"
" Anciaux, G. A Fourier-accelerated volume"
" integral method for elastoplastic contact."
" Computer Methods in Applied Mechanics and"
" Engineering 351, 951–976 (2019)."),
doi="10.1016/j.cma.2019.04.006",
),
],
),
(
("EPICSolver", ),
[
Publication(
citation=("Frérot, L., Bonnet, M., Molinari, J.-F. &"
" Anciaux, G. A Fourier-accelerated volume"
" integral method for elastoplastic contact."
" Computer Methods in Applied Mechanics and"
" Engineering 351, 951–976 (2019)."),
doi="10.1016/j.cma.2019.04.006",
),
Publication(
citation=(
"Jacq, C., Nélias, D., Lormand, G. & Girodin, D."
" Development of a Three-Dimensional"
" Semi-Analytical Elastic-Plastic Contact Code."
" Journal of Tribology 124, 653 (2002)."),
doi="10.1115/1.1467920",
),
],
),
(
("Condat", ),
[
Publication(
citation=(
"Condat, L. A Primal–Dual Splitting Method for"
" Convex Optimization Involving Lipschitzian,"
" Proximable and Linear Composite Terms. J Optim"
" Theory Appl 158, 460–479 (2012)."),
doi="10.1007/s10957-012-0245-9",
),
],
),
(
("KatoSaturated", ),
[
Publication(
citation=(
"Stanley, H. M. & Kato, T. An FFT-Based Method for"
" Rough Surface Contact. J. Tribol 119, 481–485"
" (1997)."),
doi="10.1115/1.2833523",
),
Publication(
citation=(
"Almqvist, A., Sahlin, F., Larsson, R. &"
" Glavatskih, S. On the dry elasto-plastic contact"
" of nominally flat surfaces. Tribology"
" International 40, 574–579 (2007)."),
doi="10.1016/j.triboint.2005.11.008",
),
],
),
] for k in keys
}
frame = inspect.stack()[1][0]
caller_vars = (frame.f_globals | frame.f_locals)
citable = filter(
partial(contains, _publications.keys()),
map(lambda x: type(x).__name__, caller_vars.values()),
)
citable = [joss, zenodo] \
+ list({pub for k in citable for pub in _publications[k]})
msg = "Please cite the following publications:\n\n"
msg += "\n".join(format_str.format(pub=pub) for pub in citable)
Logger().get(LogLevel.info) << msg
return citable
def load_path(
solver: ContactSolver,
loads: Iterable[Union[float, np.ndarray]],
verbose: bool = False,
callback=None,
) -> Iterable[Model]:
"""
Generate model objects solutions for a sequence of applied loads.
:param solver: a contact solver object
:param loads: an iterable sequence of loads
:param verbose: print info output of solver
:param callback: a callback executed after the yield
"""
log_level = LogLevel.info if verbose else LogLevel.warning
with log_context(log_level):
for load in loads:
if solver.solve(load) > solver.tolerance:
raise NoConvergenceError("Solver error exceeded tolerance")
yield solver.model
if callback is not None:
callback()
def seeded_surfaces(
generator: Union[SurfaceGenerator1D, SurfaceGenerator2D],
seeds: Iterable[int],
) -> Iterable[np.ndarray]:
"""
Generate rough surfaces with a prescribed seed sequence.
:param generator: surface generator object
:param seeds: random seed sequence
"""
for seed in seeds:
generator.random_seed = seed
yield generator.buildSurface()
def hertz_surface(system_size: Iterable[float], shape: Iterable[int],
radius: float) -> np.ndarray:
"""
Construct a parabolic surface.
:param system_size: size of the domain in each direction
:param shape: number of points in each direction
:param radius: radius of surface
"""
coords = map(
lambda L, N: np.linspace(0, L, N, endpoint=False, dtype=dtype),
system_size,
shape,
)
coords = np.meshgrid(*coords, "ij", sparse=True)
surface = (-1 / (2 * radius)) * sum(
(x - L / 2)**2 for x, L in zip(coords, system_size))
return surface
def radial_average(x, y, values, r, theta, method='linear', endpoint=False):
"""Compute the radial average of a 2D field."""
try:
from scipy.interpolate import RegularGridInterpolator
except ImportError:
raise ImportError("Install scipy to use tamaas.utils.radial_average")
interpolator = RegularGridInterpolator((x, y), values, method=method)
rr, tt = np.meshgrid(r, theta, indexing='ij', sparse=True)
x, y = rr * np.cos(tt), rr * np.sin(tt)
X = np.vstack((x.flatten(), y.flatten())).T
return interpolator(X).reshape(x.shape).sum(axis=1) \
/ (theta.size - int(not endpoint))
diff --git a/python/tamaas_module.cpp b/python/tamaas_module.cpp
index 554dfe7..ce6df99 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,87 +1,87 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "tamaas.hh"
#include "tamaas_info.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace py = pybind11;
namespace detail {
template
struct dtype_helper {
static const py::dtype dtype;
};
template <>
const py::dtype dtype_helper::dtype("=f8");
template <>
const py::dtype dtype_helper::dtype("=f16");
} // namespace detail
/// Creating the tamaas python module
PYBIND11_MODULE(_tamaas, mod) {
mod.doc() = "Compiled component of Tamaas";
// Wrapping the base methods
mod.def("initialize", &initialize, py::arg("num_threads") = 0,
"Initialize tamaas with desired number of threads. Automatically "
"called upon import of the tamaas module, but can be manually called "
"to set the desired number of threads.");
mod.def("finalize", []() {
PyErr_WarnEx(PyExc_DeprecationWarning,
"finalize() is deprecated, it is now automatically called "
"at the end of the program", 1);
});
// Default dtype of numpy arrays
mod.attr("dtype") = detail::dtype_helper::dtype;
// Wrapping release information
auto info = py::class_(mod, "TamaasInfo");
info.attr("version") = TamaasInfo::version;
info.attr("build_type") = TamaasInfo::build_type;
info.attr("branch") = TamaasInfo::branch;
info.attr("commit") = TamaasInfo::commit;
info.attr("diff") = TamaasInfo::diff;
info.attr("remotes") = TamaasInfo::remotes;
info.attr("has_mpi") = TamaasInfo::has_mpi;
info.attr("backend") = TamaasInfo::backend;
// Wrapping tamaas components
wrap::wrapCore(mod);
wrap::wrapPercolation(mod);
wrap::wrapSurface(mod);
wrap::wrapModel(mod);
wrap::wrapSolvers(mod);
wrap::wrapCompute(mod);
wrap::wrapMPI(mod);
/// Wrapping test features
wrap::wrapTestFeatures(mod);
}
} // namespace tamaas
diff --git a/python/wrap.hh b/python/wrap.hh
index 4a18d8f..8998497 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,77 +1,77 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef WRAP_HH
#define WRAP_HH
/* -------------------------------------------------------------------------- */
#include "cast.hh"
#include "numpy.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// For naming classes templated with dimension
inline std::string makeDimensionName(const std::string& name, UInt dim) {
std::stringstream str;
str << name << dim << "D";
return str.str();
}
/* -------------------------------------------------------------------------- */
/* Prototypes for wrapping of tamaas components */
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod);
void wrapPercolation(py::module& mod);
void wrapSurface(py::module& mod);
void wrapModel(py::module& mod);
void wrapSolvers(py::module& mod);
void wrapTestFeatures(py::module& mod);
void wrapCompute(py::module& mod);
void wrapMPI(py::module& mod);
/* -------------------------------------------------------------------------- */
/* Deprecation macro */
/* -------------------------------------------------------------------------- */
#define TAMAAS_DEPRECATE(olds, news) \
do { \
PyErr_WarnEx(PyExc_DeprecationWarning, \
olds " is deprecated, use " news " instead.", 1); \
} while (0)
#define TAMAAS_DEPRECATE_ACCESSOR(acc, type, property) \
#acc, [](const type& m) -> decltype(m.acc()) { \
TAMAAS_DEPRECATE(#acc "()", "the " property " property"); \
return m.acc(); \
}
} // namespace wrap
} // namespace tamaas
#endif // WRAP_HH
diff --git a/python/wrap/compute.cpp b/python/wrap/compute.cpp
index fee8395..8900c2e 100644
--- a/python/wrap/compute.cpp
+++ b/python/wrap/compute.cpp
@@ -1,91 +1,91 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "computes.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
void wrapCompute(py::module& mod) {
auto compute_mod = mod.def_submodule("compute");
compute_mod.def(
"eigenvalues",
[](model_type type, Grid& eigs, const Grid& field) {
eigenvalues(type, eigs, field);
},
"model_type"_a, "eigenvalues_out"_a, "field"_a,
"Compute eigenvalues of a tensor field");
compute_mod.def(
"von_mises",
[](model_type type, Grid& vm, const Grid& field) {
vonMises(type, vm, field);
},
"model_type"_a, "von_mises"_a, "field"_a,
"Compute the Von Mises invariant of a tensor field");
compute_mod.def(
"deviatoric",
[](model_type type, Grid& dev, const Grid& field) {
deviatoric(type, dev, field);
},
"model_type"_a, "deviatoric"_a, "field"_a,
"Compute the deviatoric part of a tensor field");
compute_mod.def(
"to_voigt",
[](const Grid& field) {
if (field.getNbComponents() != 9)
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
Grid voigt(field.sizes(), 6);
Loop::loop(
[] CUDA_LAMBDA(MatrixProxy in,
SymMatrixProxy out) { out.symmetrize(in); },
range>(field),
range>(voigt));
return voigt;
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
compute_mod.def(
"from_voigt",
[](const Grid& field) {
if (field.getNbComponents() != 6)
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
Grid full(field.sizes(), 9);
Loop::loop([] CUDA_LAMBDA(
SymMatrixProxy in,
MatrixProxy out) { out.fromSymmetric(in); },
range>(field),
range>(full));
return full;
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index 2b5f5d7..4b889ed 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,123 +1,123 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "logger.hh"
#include "statistics.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
void wrapStatistics(py::module& mod) {
auto name = makeDimensionName("Statistics", dim);
py::class_>(mod, name.c_str())
.def_static("computePowerSpectrum",
&Statistics::computePowerSpectrum,
py::return_value_policy::move, "Compute PSD of surface")
.def_static(
"computeAutocorrelation", &Statistics::computeAutocorrelation,
py::return_value_policy::move, "Compute autocorrelation of surface")
.def_static("computeMoments", &Statistics::computeMoments,
"Compute spectral moments")
.def_static("computeSpectralRMSSlope",
&Statistics::computeSpectralRMSSlope,
"Compute hrms' in Fourier space")
.def_static("computeFDRMSSlope", &Statistics::computeFDRMSSlope,
"Compute hrms' with finite differences")
.def_static("computeRMSHeights", &Statistics::computeRMSHeights,
"Compute hrms")
.def_static("contact", &Statistics::contact, "tractions"_a,
"perimeter"_a = 0,
"Compute the (corrected) contact area. Permieter is the "
"total contact perimeter in number of segments.")
.def_static("graphArea", &Statistics::graphArea, "zdisplacement"_a,
"Compute area defined by function graph.");
}
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod) {
// Exposing logging facility
py::enum_(mod, "LogLevel")
.value("debug", LogLevel::debug)
.value("info", LogLevel::info)
.value("warning", LogLevel::warning)
.value("error", LogLevel::error);
mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); });
mod.def("get_log_level", Logger::getCurrentLevel);
py::class_(mod, "Logger")
.def(py::init<>())
.def(
"get",
[](Logger& logger, LogLevel level) -> Logger& {
logger.get(level);
return logger;
},
"Get a logger object for a log level")
.def(
"__lshift__",
[](Logger& logger, std::string msg) -> Logger& {
auto level = logger.getWishLevel();
if (level >= Logger::getCurrentLevel() and
not(mpi::rank() != 0 and level == LogLevel::info)) {
py::print(msg,
"file"_a = py::module::import("sys").attr("stderr"));
}
return logger;
},
"Log a message");
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
mod.def(
"to_voigt",
[](const Grid& field) {
TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()");
if (field.getNbComponents() == 9) {
Grid voigt(field.sizes(), 6);
Loop::loop([] CUDA_LAMBDA(
MatrixProxy in,
SymMatrixProxy out) { out.symmetrize(in); },
range>(field),
range>(voigt));
return voigt;
}
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 7d37f60..7731b1b 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,538 +1,538 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "adhesion_functional.hh"
#include "elastic_functional.hh"
#include "field_container.hh"
#include "functional.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_extensions.hh"
#include "model_factory.hh"
#include "numpy.hh"
#include "residual.hh"
#include "wrap.hh"
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
struct model_operator_accessor {
Model& m;
decltype(auto) get(const std::string& name) {
return m.getIntegralOperator(name);
}
};
template
std::unique_ptr> instanciateFromNumpy(numpy& num) {
std::unique_ptr> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique>>(num);
return result;
case 3:
result = std::make_unique>>(num);
return result;
case 4:
result = std::make_unique>>(num);
return result;
default:
TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
"array to be the number of components");
}
}
void wrapFieldContainer(py::module& mod) {
py::class_(mod, "FieldContainer")
.def(py::init<>())
.def("__getitem__", &FieldContainer::at,
py::return_value_policy::reference_internal)
.def(
"__setitem__",
[](FieldContainer& f, const std::string& key, numpy grid) {
f[key] = instanciateFromNumpy(grid);
},
py::keep_alive<1, 3>())
.def(
"__setitem__",
[](FieldContainer& f, const std::string& key, numpy grid) {
f[key] = instanciateFromNumpy(grid);
},
py::keep_alive<1, 3>())
.def(
"__setitem__",
[](FieldContainer& f, const std::string& key, numpy grid) {
f[key] = instanciateFromNumpy(grid);
},
py::keep_alive<1, 3>())
.def(
"__contains__",
[](const FieldContainer& m, const std::string& key) {
const auto fields = m.fields();
return std::find(fields.begin(), fields.end(), key) != fields.end();
},
"Test field existence")
.def(
"__iter__",
[](const FieldContainer& m) {
const auto& fields = m.fields_map();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>(), "Iterator on fields")
/// Deprecated interface
.def(
"registerField",
[](Model& m, const std::string& name, numpy field) {
TAMAAS_DEPRECATE("registerField()", "the [] operator");
m[name] = instanciateFromNumpy(field);
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def(
"getField",
[](const Model& m, const std::string& name) -> decltype(m.at(name)) {
TAMAAS_DEPRECATE("getField()", "the [] operator");
return m.at(name);
},
"field_name"_a, py::return_value_policy::reference_internal)
.def(
"getFields",
[](const Model& m) {
TAMAAS_DEPRECATE("getFields()", "list(model)");
return m.fields();
},
"Return fields list");
}
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_,
functional::wrap::PyFunctional>
func(mod, "Functional");
func.def(py::init<>())
.def("computeF", &functional::Functional::computeF,
"Compute functional value")
.def("computeGradF", &functional::Functional::computeGradF,
"Compute functional gradient");
py::class_(
mod, "ElasticFunctionalPressure", func)
.def(py::init&>());
py::class_(mod, "ElasticFunctionalGap",
func)
.def(py::init&>());
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters,
"Parameters dictionary")
.def("setParameters", [](functional::AdhesionFunctional& f,
const std::map& m) {
TAMAAS_DEPRECATE("setParameters()", "the parameters property");
f.setParameters(m);
});
py::class_(
mod, "ExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-g/ρ)")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh,
"Cohesive zone potential F = H(g - ρ)·γ/ρ")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-0.5·(g/ρ)²)")
.def(py::init&>(), "surface"_a);
}
/// Wrap IntegralOperator
void wrapIntegralOperator(py::module& mod) {
py::class_(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy input, numpy output) {
TAMAAS_DEPRECATE("apply()", "the () operator");
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"),
py::return_value_policy::reference)
.def(TAMAAS_DEPRECATE_ACCESSOR(getKind, IntegralOperator, "kind"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getType, IntegralOperator, "type"))
.def(
"__call__",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
},
"Apply the integral operator")
.def("updateFromModel", &IntegralOperator::updateFromModel,
"Resets internal persistent variables from the model")
.def_property_readonly("kind", &IntegralOperator::getKind)
.def_property_readonly("model", &IntegralOperator::getModel)
.def_property_readonly("type", &IntegralOperator::getType)
.def_property_readonly("shape", &IntegralOperator::matvecShape)
.def(
"matvec",
[](const IntegralOperator& op, numpy X) -> GridBase {
GridBaseNumpy x(X);
auto y = op.matvec(x);
return y;
},
py::return_value_policy::move);
py::enum_(mod, "integration_method",
"Integration method used for the computation "
"of volumetric Fourier operators")
.value("linear", integration_method::linear,
"No approximation error, O(N₁·N₂·N₃) time complexity, may cause "
"float overflow/underflow")
.value("cutoff", integration_method::cutoff,
"Approximation, O(sqrt(N₁²+N₂²)·N₃²) time complexity, no "
"overflow/underflow risk");
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet)
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, BEEngine, "model"),
py::return_value_policy::reference)
.def_property_readonly("model", &BEEngine::getModel);
}
template
void wrapModelTypeTrait(py::module& mod) {
using trait = model_type_traits;
py::class_(mod, trait::repr)
.def_property_readonly_static(
"dimension", [](const py::object&) { return trait::dimension; },
"Dimension of computational domain")
.def_property_readonly_static(
"components", [](const py::object&) { return trait::components; },
"Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
[](const py::object&) { return trait::boundary_dimension; },
"Dimension of boundary of computational domain")
.def_property_readonly_static(
"voigt", [](const py::object&) { return trait::voigt; },
"Number of components of symmetrical tensor fields")
.def_property_readonly_static(
"indices", [](const py::object&) { return trait::indices; });
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
.value("basic_1d", model_type::basic_1d,
"Normal contact with 1D interface")
.value("basic_2d", model_type::basic_2d,
"Normal contact with 2D interface")
.value("surface_1d", model_type::surface_1d,
"Normal & tangential contact with 1D interface")
.value("surface_2d", model_type::surface_2d,
"Normal & tangential contact with 2D interface")
.value("volume_1d", model_type::volume_1d,
"Contact with volumetric representation and 1D interface")
.value("volume_2d", model_type::volume_2d,
"Contact with volumetric representation and 2D interface");
auto trait_mod = mod.def_submodule("_type_traits");
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
py::class_(mod, "_model_operator_acessor")
.def(py::init())
.def(
"__getitem__",
[](model_operator_accessor& acc, std::string name) {
try {
return acc.get(name);
} catch (std::out_of_range&) {
throw py::key_error(name);
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](model_operator_accessor& acc, std::string key) {
const auto ops = acc.m.getIntegralOperators();
return std::find(ops.begin(), ops.end(), key) != ops.end();
})
.def(
"__iter__",
[](const model_operator_accessor& acc) {
const auto& ops = acc.m.getIntegralOperatorsMap();
return py::make_key_iterator(ops.cbegin(), ops.cend());
},
py::keep_alive<0, 1>());
py::class_(mod, "Model")
.def(py::init(py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel)),
R"-(Create a new model of a given type, physical size and *global*
discretization.
:param model_type: the type of desired model
:param system_size: the physical size of the domain in each direction
:param global_discretization: number of points in each direction)-")
.def_property_readonly("type", &Model::getType)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus,
"Young's modulus")
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio,
"Poisson's ratio")
.def_property_readonly("mu", &Model::getShearModulus, "Shear modulus")
.def_property_readonly("E_star", &Model::getHertzModulus,
"Contact (Hertz) modulus")
.def_property_readonly("be_engine", &Model::getBEEngine,
"Boundary element engine")
.def(
"setElasticity",
[](Model& m, Real E, Real nu) {
TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties");
m.setElasticity(E, nu);
},
"E"_a, "nu"_a)
.def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model,
"boundary_system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model,
"boundary_shape"))
.def("solveNeumann", &Model::solveNeumann,
"Solve surface tractions -> displacements")
.def("solveDirichlet", &Model::solveDirichlet,
"Solve surface displacemnts -> tractions")
.def("dump", &Model::dump, "Write model data to registered dumpers")
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>(),
"Register a dumper")
.def(
"getBEEngine",
[](Model& m) -> decltype(m.getBEEngine()) {
TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property");
return m.getBEEngine();
},
py::return_value_policy::reference_internal)
.def(
"getIntegralOperator",
[](const Model& m, std::string name) {
TAMAAS_DEPRECATE("getIntegralOperator()", "the operators property");
return m.getIntegralOperator(name);
},
"operator_name"_a, py::return_value_policy::reference_internal)
.def(
"applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
},
"Apply Hooke's law")
// Python magic functions
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def(
"__copy__",
[](const Model&) {
throw std::runtime_error("__copy__ not implemented");
},
"Shallow copy of model. Not implemented.")
.def(
"__deepcopy__",
[](const Model& m, py::dict) { return ModelFactory::createModel(m); },
"Deep copy of model.")
.def_property_readonly("boundary_fields", &Model::getBoundaryFields)
.def_property_readonly(
"operators", [](Model& m) { return model_operator_accessor{m}; },
"Returns a dict-like object allowing access to the model's "
"integral "
"operators")
// More python-like access to model properties
.def_property_readonly("shape", &Model::getDiscretization,
"Discretization (local in MPI environment)")
.def_property_readonly("global_shape", &Model::getGlobalDiscretization,
"Global discretization (in MPI environement)")
.def_property_readonly("boundary_shape",
&Model::getBoundaryDiscretization,
"Number of points on boundary")
.def_property_readonly("system_size", &Model::getSystemSize,
"Size of physical domain")
.def_property_readonly("boundary_system_size",
&Model::getBoundarySystemSize,
"Physical size of surface")
.def_property_readonly("traction",
(const GridBase& (Model::*)() const) &
Model::getTraction,
"Surface traction field")
.def_property_readonly("displacement",
(const GridBase& (Model::*)() const) &
Model::getDisplacement,
"Displacement field");
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a, "Dump model")
.def(
"__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; },
"Dump model");
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static(
"createModel",
py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel),
"model_type"_a, "system_size"_a, "global_discretization"_a,
R"-(Create a new model of a given type, physical size and *global*
discretization.
:param model_type: the type of desired model
:param system_size: the physical size of the domain in each direction
:param global_discretization: number of points in each direction)-")
.def_static("createModel",
py::overload_cast(&ModelFactory::createModel),
"model"_a, "Create a deep copy of a model.")
.def_static("createResidual", &ModelFactory::createResidual, "model"_a,
"sigma_y"_a, "hardening"_a = 0.,
R"-(Create an isotropic linear hardening residual.
:param model: the model on which to define the residual
:param sigma_y: the (von Mises) yield stress
:param hardening: the hardening modulus)-")
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a,
"Register Boussinesq and Mindlin operators to model.")
.def_static("setIntegrationMethod", &ModelFactory::setIntegrationMethod,
"operator"_a, "method"_a, "cutoff"_a,
"Set the integration method (linear or cutoff) for a volume "
"integral operator");
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def(
"applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
"cutoff"_a = 1e-12)
.def_property("yield_stress", &Residual::getYieldStress,
&Residual::setYieldStress)
.def_property("hardening_modulus", &Residual::getHardeningModulus,
&Residual::setHardeningModulus)
.def_property_readonly("model", &Residual::getModel);
}
void wrapModel(py::module& mod) {
wrapFieldContainer(mod);
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/model_extensions.hh b/python/wrap/model_extensions.hh
index 09debbd..10679e1 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,148 +1,148 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "functional.hh"
#include "model.hh"
#include "model_dumper.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace functional {
namespace wrap {
/// Class for extension of Functional in python
class PyFunctional : public Functional {
public:
using Functional::Functional;
// Overriding pure virtual functions
Real computeF(GridBase& variable, GridBase& dual) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
}
void computeGradF(GridBase& variable,
GridBase& gradient) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
}
};
} // namespace wrap
} // namespace functional
/* -------------------------------------------------------------------------- */
namespace wrap {
/* -------------------------------------------------------------------------- */
class PyModelDumper : public ModelDumper {
public:
using ModelDumper::ModelDumper;
void dump(const Model& model) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, ModelDumper, dump, model);
}
};
/* -------------------------------------------------------------------------- */
class PyResidual : public Residual {
public:
using Residual::Residual;
void computeResidual(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidual, strain_increment);
}
void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, applyTangent, output, input,
current_strain_increment);
}
void computeStress(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeStress, strain_increment);
}
void updateState(GridBase& converged_strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, updateState,
converged_strain_increment);
}
const GridBase& getVector() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getVector);
}
const GridBase& getPlasticStrain() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getPlasticStrain);
}
const GridBase& getStress() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getStress);
}
void computeResidualDisplacement(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidualDisplacement,
strain_increment);
}
void setIntegrationMethod(integration_method method, Real cutoff) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setIntegrationMethod, method,
cutoff);
}
Real getYieldStress() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Residual, getYieldStress);
}
Real getHardeningModulus() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Residual, getHardeningModulus);
}
void setYieldStress(Real val) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setYieldStress, val);
}
void setHardeningModulus(Real val) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setHardeningModulus, val);
}
};
/* -------------------------------------------------------------------------- */
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/mpi.cpp b/python/wrap/mpi.cpp
index b36ad30..18d8cdd 100644
--- a/python/wrap/mpi.cpp
+++ b/python/wrap/mpi.cpp
@@ -1,101 +1,101 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "mpi_interface.hh"
#include "partitioner.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
void wrapMPI(py::module& mod) {
auto mpi_mod = mod.def_submodule("mpi");
py::class_(mpi_mod, "sequential")
.def(py::init<>())
.def("__enter__",
[](mpi::sequential& /*obj*/) { mpi::sequential::enter(); })
.def("__exit__", [](mpi::sequential& /*obj*/, py::object /*unused*/,
py::object /*unused*/,
py::object /*unused*/) { mpi::sequential::exit(); });
mpi_mod.def(
"local_shape",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_size(global);
case 2:
return Partitioner<2>::local_size(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a, "Gives the local size of a 1D/2D global shape");
mpi_mod.def(
"global_shape",
[](std::vector local) {
switch (local.size()) {
case 1:
return Partitioner<1>::global_size(local);
case 2:
return Partitioner<2>::global_size(local);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"local_shape"_a, "Gives the global shape of a 1D/2D local shape");
mpi_mod.def(
"local_offset",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_offset(global);
case 2:
return Partitioner<2>::local_offset(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a, "Gives the local offset of a 1D/2D global shape");
mpi_mod.def("gather", &Partitioner<2>::gather,
py::return_value_policy::move, "Gather 2D surfaces");
mpi_mod.def("scatter", &Partitioner<2>::scatter,
py::return_value_policy::move, "Scatter 2D surfaces");
mpi_mod.def(
"size", []() { return mpi::size(); },
"Returns the number of MPI processes");
mpi_mod.def(
"rank", []() { return mpi::rank(); },
"Returns the rank of the local process");
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/percolation.cpp b/python/wrap/percolation.cpp
index cfff21c..8be4fa5 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,95 +1,95 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "flood_fill.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
/* -------------------------------------------------------------------------- */
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
void wrapCluster(py::module& mod) {
auto name = makeDimensionName("Cluster", dim);
py::class_>(mod, name.c_str())
.def(py::init<>())
.def_property_readonly("area", &Cluster::getArea, "Area of cluster")
.def_property_readonly("points", &Cluster::getPoints,
"Get list of points of cluster")
.def_property_readonly("perimeter", &Cluster::getPerimeter,
"Get perimeter of cluster")
.def_property_readonly("bounding_box", &Cluster::boundingBox,
"Compute the bounding box of a cluster")
.def_property_readonly("extent", &Cluster::extent,
"Compute the extents of a cluster")
.def(TAMAAS_DEPRECATE_ACCESSOR(getArea, Cluster, "area"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPoints, Cluster, "points"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPerimeter, Cluster, "perimeter"))
.def("__str__", [](const Cluster& cluster) {
std::stringstream sstr;
sstr << '{';
auto&& points = cluster.getPoints();
auto print_point = [&sstr](const auto& point) {
sstr << '(';
for (UInt i = 0; i < point.size() - 1; ++i)
sstr << point[i] << ", ";
sstr << point.back() << ')';
};
auto point_it = points.begin();
for (UInt p = 0; p < points.size() - 1; ++p) {
print_point(*(point_it++));
sstr << ", ";
}
print_point(points.back());
sstr << "}";
return sstr.str();
});
}
void wrapPercolation(py::module& mod) {
wrapCluster<1>(mod);
wrapCluster<2>(mod);
wrapCluster<3>(mod);
py::class_(mod, "FloodFill")
.def_static("getSegments", &FloodFill::getSegments,
"Return a list of segments from boolean map", "contact"_a)
.def_static("getClusters", &FloodFill::getClusters,
"Return a list of clusters from boolean map", "contact"_a,
"diagonal"_a)
.def_static("getVolumes", &FloodFill::getVolumes,
"Return a list of volume clusters", "map"_a, "diagonal"_a);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 3f9e3c0..aacdce7 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,239 +1,239 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "beck_teboulle.hh"
#include "condat.hh"
#include "contact_solver.hh"
#include "dfsane_solver.hh"
#include "ep_solver.hh"
#include "epic.hh"
#include "kato.hh"
#include "kato_saturated.hh"
#include "polonsky_keer_rey.hh"
#include "polonsky_keer_tan.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
/* -------------------------------------------------------------------------- */
using namespace py::literals;
class PyEPSolver : public EPSolver {
public:
using EPSolver::EPSolver;
// NOLINTNEXTLINE(readability-else-after-return)
void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
void updateState() override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD(void, EPSolver, updateState);
}
};
class PyContactSolver : public ContactSolver {
public:
using ContactSolver::ContactSolver;
Real solve(std::vector load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
Real solve(Real load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_(mod, "ContactSolver")
.def(py::init&, Real>(),
py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def(
"setMaxIterations",
[](ContactSolver& m, UInt iter) {
TAMAAS_DEPRECATE("setMaxIterations()", "the max_iter property");
m.setMaxIterations(iter);
},
"max_iter"_a)
.def(
"setDumpFrequency",
[](ContactSolver& m, UInt freq) {
TAMAAS_DEPRECATE("setDumpFrequency()", "the dump_freq property");
m.setDumpFrequency(freq);
},
"dump_freq"_a)
.def_property("tolerance", &ContactSolver::getTolerance,
&ContactSolver::setTolerance, "Solver tolerance")
.def_property("max_iter", &ContactSolver::getMaxIterations,
&ContactSolver::setMaxIterations,
"Maximum number of iterations")
.def_property("dump_freq", &ContactSolver::getDumpFrequency,
&ContactSolver::setDumpFrequency,
"Frequency of displaying solver info")
.def_property_readonly("model", &ContactSolver::getModel)
.def_property_readonly("surface", &ContactSolver::getSurface)
.def_property_readonly("functional", &ContactSolver::getFunctional)
.def("addFunctionalTerm", &ContactSolver::addFunctionalTerm,
"Add a term to the contact functional to minimize")
.def("solve", py::overload_cast>(&ContactSolver::solve),
"target_force"_a,
py::call_guard(),
"Solve the contact for a mean traction/gap vector")
.def("solve", py::overload_cast(&ContactSolver::solve),
"target_normal_pressure"_a,
py::call_guard(),
"Solve the contact for a mean normal pressure/gap");
py::class_ pkr(
mod, "PolonskyKeerRey",
"Main solver class for normal elastic contact problems. Its functional "
"can be customized to add an adhesion term, and its primal variable can "
"be set to either the gap or the pressure.");
// Need to export enum values before defining PKR constructor
py::enum_(pkr, "type")
.value("gap", PolonskyKeerRey::gap)
.value("pressure", PolonskyKeerRey::pressure)
.export_values();
pkr.def(py::init&, Real, PolonskyKeerRey::type,
PolonskyKeerRey::type>(),
"model"_a, "surface"_a, "tolerance"_a,
"primal_type"_a = PolonskyKeerRey::type::pressure,
"constraint_type"_a = PolonskyKeerRey::type::pressure,
py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def("computeError", &PolonskyKeerRey::computeError);
py::class_(
mod, "KatoSaturated",
"Solver for pseudo-plasticity problems where the normal pressure is "
"constrained above by a saturation pressure \"pmax\"")
.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "pmax"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def_property("pmax", &KatoSaturated::getPMax, &KatoSaturated::setPMax,
"Saturation normal pressure");
py::class_ kato(mod, "Kato");
kato.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &Kato::solve, "p0"_a, "proj_iter"_a = 50)
.def("solveRelaxed", &Kato::solveRelaxed, "g0"_a)
.def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a = 0.01)
.def("computeCost", &Kato::computeCost, "use_tresca"_a = false);
py::class_ bt(mod, "BeckTeboulle");
bt.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &BeckTeboulle::solve, "p0"_a)
.def("computeCost", &BeckTeboulle::computeCost);
py::class_ cd(
mod, "Condat",
"Main solver for frictional contact problems. It has no restraint on the "
"material properties or friction coefficient values, but solves an "
"associated version of the Coulomb friction law, which differs from the "
"traditional Coulomb friction in that the normal and tangential slip "
"components are coupled.");
cd.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &Condat::solve, "p0"_a, "grad_step"_a = 0.9)
.def("computeCost", &Condat::computeCost);
py::class_ pkt(mod, "PolonskyKeerTan");
pkt.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &PolonskyKeerTan::solve, "p0"_a)
.def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a)
.def("computeCost", &PolonskyKeerTan::computeCost,
"use_tresca"_a = false);
py::class_(
mod, "_tolerance_manager",
"Manager object for the tolereance of nonlinear plasticity solvers. "
"Decreases the solver tolerance by geometric progression.")
.def(py::init(), "start_tol"_a, "end_tol"_a, "rate"_a);
py::class_(
mod, "EPSolver", "Mother class for nonlinear plasticity solvers")
.def(py::init(), "residual"_a, py::keep_alive<1, 2>())
.def("solve", &EPSolver::solve)
.def("getStrainIncrement", &EPSolver::getStrainIncrement,
py::return_value_policy::reference_internal)
.def("getResidual", &EPSolver::getResidual,
py::return_value_policy::reference_internal)
.def("updateState", &EPSolver::updateState)
.def_property("tolerance", &EPSolver::getTolerance,
&EPSolver::setTolerance)
.def("setToleranceManager", &EPSolver::setToleranceManager)
.def("beforeSolve", &EPSolver::beforeSolve);
py::class_(mod, "_DFSANESolver")
.def(py::init(), "residual"_a, py::keep_alive<1, 2>())
.def(py::init([](Residual& res, Model&) {
return std::make_unique(res);
}),
"residual"_a, "model"_a, py::keep_alive<1, 2>(),
"Deprecated constructor: only here for compatibility reasons");
py::class_(
mod, "EPICSolver",
"Main solver class for elastic-plastic contact problems")
.def(py::init(),
"contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10,
"relaxation"_a = 0.3, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def(
"solve",
[](EPICSolver& solver, Real pressure) {
return solver.solve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard(),
"Solves the EP contact with a relaxed fixed-point scheme. Adjust "
"the relaxation parameter to help convergence.")
.def(
"acceleratedSolve",
[](EPICSolver& solver, Real pressure) {
return solver.acceleratedSolve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard(),
"Solves the EP contact with an accelerated fixed-point scheme. May "
"not converge!");
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 2f8a0f4..f1fd979 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,178 +1,178 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- * Copyright (©) 2020-2022 Lucas Frérot
+ * Copyright (©) 2020-2023 Lucas Frérot
*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include "regularized_powerlaw.hh"
#include "surface_generator_filter.hh"
#include "surface_generator_random_phase.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
class PyFilter : public Filter {
public:
using Filter::Filter;
// Overriding pure virtual functions
void
computeFilter(GridHermitian& filter_coefficients) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Filter