diff --git a/INFOS.py b/INFOS.py
index 11d1f5b..267076f 100644
--- a/INFOS.py
+++ b/INFOS.py
@@ -1,34 +1,34 @@
# -*- mode:python; coding: utf-8 -*-
"Defines the information to be used throughout the builds"
from collections import namedtuple
import versioneer
TamaasInfo = namedtuple('TamaasInfo',
['version',
'authors',
'maintainer',
'email',
'copyright',
'description'])
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='lfrerot1@jhu.edu',
copyright=(
- u"Copyright (©) 2016-2021 EPFL "
+ u"Copyright (©) 2016-2022 EPFL "
u"(École Polytechnique Fédérale de Lausanne), "
u"Laboratory (LSMS - Laboratoire de Simulation en "
u"Mécanique des Solides)"
),
description='A high-performance library for periodic rough surface contact',
)
diff --git a/SConstruct b/SConstruct
index a39a6c0..ceff295 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,474 +1,474 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 .
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
from __future__ import print_function
import sys
import os
from subprocess import check_output
# Import below not strictly necessary, but good for pep8
from SCons.Script import (
EnsurePythonVersion,
EnsureSConsVersion,
Help,
Environment,
Variables,
EnumVariable,
PathVariable,
BoolVariable,
ListVariable,
Split,
Export,
Dir,
)
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(2, 7)
EnsureSConsVersion(2, 4)
# ------------------------------------------------------------------------------
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/preprocessor/seq.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', 'omp',
allowed_values=('cpp', 'omp', 'tbb', 'cuda'),
ignorecase=2),
EnumVariable('fftw_threads', 'Threads FFTW library preference', 'omp',
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'),
# 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),
# Dependencies provided as submodule get different default
PathVariable('GTEST_ROOT',
'Googletest custom path',
os.getenv('GTEST_ROOT', '#third-party/googletest/googletest'),
PathVariable.PathAccept),
PathVariable('PYBIND11_ROOT',
'Pybind11 custom path',
os.getenv('PYBIND11_ROOT', '#third-party/pybind11/include'),
PathVariable.PathAccept),
PathVariable('EXPOLIT_ROOT',
'Expolit custom path',
os.getenv('EXPOLIT_ROOT', '#third-party/expolit/include'),
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', False),
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', False),
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]))
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'])
# 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')
main_env.CompilationDatabase(PRINT_CMD_LINE_FUNC=pretty_cmd_print)
# 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 90bce9c..2b70a47 100644
--- a/doc/SConscript
+++ b/doc/SConscript
@@ -1,104 +1,104 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
from __future__ import print_function
from 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)
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 876e1f2..35c6b25 100644
--- a/examples/adhesion.py
+++ b/examples/adhesion.py
@@ -1,120 +1,120 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 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 = 0
# 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()
diff --git a/examples/pipe_tools/contact b/examples/pipe_tools/contact
index eb6a325..634c3a1 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-2021, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2022, 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 d59c2f6..b6c3f6e 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-2021, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2022, 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 8aecd5d..0997600 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-2021, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2022, 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 6d680b3..11ae2b0 100644
--- a/examples/plasticity.py
+++ b/examples/plasticity.py
@@ -1,83 +1,83 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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
# 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)
# Extract local part of the global surface
local_shape = tm.mpi.local_shape(surface.shape)
offset = tm.mpi.local_offset(surface.shape)
local_surface = surface[offset:offset+local_shape[0], :].copy()
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))
diff --git a/examples/rough_contact.py b/examples/rough_contact.py
index ba49fb0..079c2f2 100644
--- a/examples/rough_contact.py
+++ b/examples/rough_contact.py
@@ -1,62 +1,62 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import tamaas as tm
import matplotlib.pyplot as plt
# Initialize threads and fftw
tm.set_log_level(tm.LogLevel.info) # Show progression of solver
# Surface size
n = 512
# Surface generator
sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 0
# 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()
diff --git a/examples/saturated.py b/examples/saturated.py
index c6017e2..a9312a5 100644
--- a/examples/saturated.py
+++ b/examples/saturated.py
@@ -1,64 +1,64 @@
# -*- coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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
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()
diff --git a/examples/statistics.py b/examples/statistics.py
index cefd45e..eebe552 100644
--- a/examples/statistics.py
+++ b/examples/statistics.py
@@ -1,81 +1,81 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import os
import numpy as np
import matplotlib.pyplot as plt
import tamaas as tm
from matplotlib.colors import LogNorm
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 = 0
# 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(psd.real, norm=LogNorm())
plt.gca().set_title('Power Spectrum Density')
plt.gcf().tight_layout()
# Computing autocorrelation and shifting for plot
acf = tm.Statistics2D.computeAutocorrelation(surface)
acf = np.fft.fftshift(acf)
plt.figure()
plt.imshow(acf)
plt.gca().set_title('Autocorrelation')
plt.gcf().tight_layout()
plt.show()
# Write the rough surface in paraview's VTK format
try:
import uvw
try:
os.mkdir('paraview')
except FileExistsError:
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
diff --git a/examples/stresses.py b/examples/stresses.py
index 40d73e2..fc06c67 100644
--- a/examples/stresses.py
+++ b/examples/stresses.py
@@ -1,120 +1,120 @@
#!/usr/bin/env python3
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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
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', 'stress_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()
diff --git a/python/SConscript b/python/SConscript
index beca60a..fe39c55 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,161 +1,161 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
from __future__ import print_function
from 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
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
]
targets.append(copy_env.Command('MANIFEST.in', '#python/MANIFEST.in',
Copy("$TARGET", "$SOURCE")))
subst_env = env_pybind.Clone(
SUBST_DICT={
'@version@': '$version',
'@authors@': str(copy_env['authors']),
'@email@': '$email',
'@description@': '$description',
# TODO change when issue with unicode fixed
# '@copyright@': '$copyright',
# '@maintainer@': '$maintainer',
}
)
subst_env.Tool('textfile')
targets.append(subst_env.Substfile('setup.py.in'))
targets.append(subst_env.Substfile('tamaas/__init__.py.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'],
PYOPTIONS=['--user'] + 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 e8b608c..07b203a 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,183 +1,183 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef CAST_HH
#define CAST_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "grid.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include
#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();
}
/**
* 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) {
#define GRID_BASE_CASE(unused1, unused2, dim) \
case dim: { \
const auto& conv = dynamic_cast&>(src); \
return grid_to_python(conv, policy, parent); \
}
switch (src.getDimension()) {
BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
default:
return nullptr;
}
#undef GRID_BASE_CASE
}
};
} // namespace detail
} // namespace pybind11
#endif // CAST_HH
diff --git a/python/numpy.hh b/python/numpy.hh
index 9cea292..0e8f591 100644
--- a/python/numpy.hh
+++ b/python/numpy.hh
@@ -1,92 +1,92 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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.wrapMemory(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.wrapMemory(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.wrapMemory(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 f0cddcf..bfa1e7a 100644
--- a/python/tamaas/__init__.py.in
+++ b/python/tamaas/__init__.py.in
@@ -1,61 +1,61 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 .
"""
@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@
# TODO Change copyright when is issue with unicode is found
-__copyright__ = u"Copyright (©) 2016-2021 EPFL " \
+__copyright__ = u"Copyright (©) 2016-2022 EPFL " \
+ u"(École Polytechnique Fédérale de Lausanne), " \
+ u"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
__license__ = "SPDX-License-Identifier: AGPL-3.0-or-later"
__maintainer__ = "Lucas Frérot"
__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 e222b44..1a23177 100755
--- a/python/tamaas/__main__.py
+++ b/python/tamaas/__main__.py
@@ -1,227 +1,227 @@
#!/usr/bin/env python3
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 sys
import io
import time
import argparse
import tamaas as tm
import numpy as np
__author__ = "Lucas Frérot"
__copyright__ = (
- "Copyright (©) 2019-2021, EPFL (École Polytechnique Fédérale de Lausanne),"
+ "Copyright (©) 2019-2022, 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)
def print_version():
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.""")
def surface(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
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):
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., 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)
def plot(args):
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():
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.)
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 5523b12..c4dc253 100644
--- a/python/tamaas/compute.py
+++ b/python/tamaas/compute.py
@@ -1,20 +1,20 @@
# -*- coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
from ._tamaas.compute import * # noqa
diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py
index aa86758..38b1009 100644
--- a/python/tamaas/dumpers/__init__.py
+++ b/python/tamaas/dumpers/__init__.py
@@ -1,493 +1,493 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 .
"""
Dumpers for the class tamaas.Model
"""
from pathlib import Path
from os import PathLike
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
from ._helper import (
step_dump, directory_dump, local_slice, _is_surface_field, _basic_types,
)
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,
}
def _create_model(attrs: ts.MutableMapping):
"Create a model from attribute dictionaty"
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):
def default(self, obj):
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):
def __init__(self, file_descriptor: FileType):
super(JSONDumper, self).__init__()
self.fd = file_descriptor
def dump_to_file(self, fd, model: Model):
json.dump(model, fd, cls=ModelJSONEncoder)
def dump(self, model: Model):
if isinstance(self.fd, (str, PathLike)):
with open(self.fd, 'w') as fd:
self.dump_to_file(fd, model)
elif isinstance(self.fd, io.TextIOBase):
self.dump_to_file(self.fd, model)
else:
raise TypeError(
f"Expected a path or file descriptor, got {self.fd}")
def read_from_file(self, fd: io.TextIOBase):
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
def read(self, file_descriptor: FileType):
fd = file_descriptor
if isinstance(fd, (str, PathLike)):
with open(fd, 'r') as fd:
return self.read_from_file(fd)
if isinstance(fd, io.TextIOBase):
return self.read_from_file(fd)
raise TypeError(f"Expected a path or file descriptor, got {type(fd)}")
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 (name or 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)
def read(self, file_descriptor: FileType):
raise NotImplementedError(
f'read() method not implemented in {type(self).__name__}')
@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):
"""Saving 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))
def read(self, 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
@directory_dump('hdf5')
@step_dump
class H5Dumper(FieldDumper):
"""Dumper to HDF5 file format"""
extension = 'h5'
def _hdf5_args(self):
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):
"""Saving 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)
dset[local_slice(field, model)] = field
# Writing metadata
for name, attr in _get_attributes(model).items():
handle.attrs[name] = attr
def read(self, file_descriptor: FileType):
"Create model from HDF5 file"
mpi_args, _ = self._hdf5_args()
with h5py.File(file_descriptor, 'r', **mpi_args) as handle:
model = _create_model(handle.attrs)
for k, v in handle.items():
if model.type in _basic_types:
v = np.asarray(v).reshape(list(v.shape) + [1])
model[k] = v[local_slice(v, model)]
return model
except ImportError:
pass
try:
import uvw # noqa
@directory_dump('paraview')
@step_dump
class UVWDumper(FieldDumper):
"""Dumper to VTK files for elasto-plastic calculations"""
extension = 'vtr'
forbidden_fields = ['traction', 'gap']
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,
)
# Iterator over fields we want to dump
fields_it = filter(lambda t: t[0] not in self.forbidden_fields,
self.get_fields(model).items())
# We make fields periodic for visualization
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
@directory_dump('netcdf')
class NetCDFDumper(FieldDumper):
"""Dumper to netCDF4 files"""
extension = "nc"
boundary_fields = ['traction', 'gap']
def _file_setup(self, grp, model: Model):
grp.createDimension('frame', 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 dump_to_file(self, file_descriptor: NameType, model: Model):
if mpi.size() > 1:
raise MPIIncompatibilityError("NetCDFDumper does not function "
"properly in parallel")
mode = 'a' if Path(file_descriptor).is_file() \
and getattr(self, 'has_setup', False) else 'w'
with Dataset(file_descriptor, mode,
format='NETCDF4_CLASSIC',
parallel=mpi.size() > 1) as rootgrp:
if rootgrp.dimensions == {}:
self._file_setup(rootgrp, model)
if self.model_info != (model.global_shape, model.type):
raise ModelError(f"Unexpected model {mode}")
self._dump_generic(rootgrp, model)
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]})")
grp.createVariable(label, 'f8',
['frame'] + dim_labels + local_dim,
zlib=True)
def _dump_generic(self, grp, model):
fields = self.get_fields(model).items()
new_frame = grp.dimensions['frame'].size
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=np.double)
def read(self, file_descriptor: NameType):
"Create model with last frame"
with Dataset(file_descriptor, 'r',
format='NETCDF4_CLASSIC') as rootgrp:
attrs = {k: rootgrp.getncattr(k) for k in rootgrp.ncattrs()}
model = _create_model(attrs)
dims = rootgrp.dimensions.keys()
for k, v in filter(lambda k: k[0] not in dims,
rootgrp.variables.items()):
v = v[-1, :]
if model.type in _basic_types:
v = np.asarray(v).reshape(list(v.shape) + [1])
model[k] = v
return model
except ImportError:
pass
diff --git a/python/tamaas/dumpers/_helper.py b/python/tamaas/dumpers/_helper.py
index 5a78236..1297dc0 100644
--- a/python/tamaas/dumpers/_helper.py
+++ b/python/tamaas/dumpers/_helper.py
@@ -1,160 +1,160 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 .
"""
Helper functions for dumpers
"""
from functools import wraps
from pathlib import Path
import numpy as np
from .. import model_type, type_traits, mpi
__all__ = ["step_dump", "directory_dump"]
_basic_types = [t for t, trait in type_traits.items() if trait.components == 1]
def _is_surface_field(field, model):
bn = model.boundary_shape
gn = mpi.global_shape(bn)
shape = list(field.shape)
# Also test number of components to weed out cases where shape is the same
# in all directions
if model.type not in _basic_types:
bn.append(type_traits[model.type].components)
gn.append(bn[-1])
# Testing works for both local and global fields
return shape == bn or shape == gn
def local_slice(field, model):
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 _is_surface_field(field, model) and len(n) > len(bn):
gshape = [n[0]] + gshape
offsets = np.concatenate(([0], offsets))
shape = bn if _is_surface_field(field, model) else n
if len(field.shape) > len(shape):
shape += field.shape[len(shape):]
def sgen(pair):
offset, size = pair
return slice(offset, offset + size, None)
def sgen_basic(pair):
offset, size = pair
return slice(offset, offset + size)
slice_gen = sgen_basic if model_type in _basic_types else sgen
return tuple(map(slice_gen, zip(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
@wraps(cls.dump)
def dump(obj, *args, **kwargs):
if mpi.rank() == 0:
directory.mkdir(parents=True, exist_ok=True)
orig_dump(obj, *args, **kwargs)
@wraps(cls.file_path.fget)
def file_path(obj):
return str(directory / orig_filepath(obj))
cls.dump = dump
cls.file_path = property(file_path)
return cls
return actual_decorator
def hdf5toVTK(inpath, outname):
"Convert HDF5 dump of a model to VTK"
import h5py as h5 # noqa
from .. import ModelFactory # noqa
from . import UVWDumper # noqa
type_translate = {
'model_type.basic_1d': model_type.basic_1d,
'model_type.basic_2d': model_type.basic_2d,
'model_type.surface_1d': model_type.surface_1d,
'model_type.surface_2d': model_type.surface_2d,
'model_type.volume_1d': model_type.volume_1d,
'model_type.volume_2d': model_type.volume_2d,
}
with h5.File(inpath, 'r') as h5file:
model_t = h5file.attrs['model_type']
system_size = list(h5file.attrs['system_size'])
n = list(h5file.attrs['discretization'])
model = ModelFactory.createModel(type_translate[model_t],
system_size, n)
fields = []
for field in h5file:
model.registerField(field, h5file[field][:])
fields.append(field)
uvw_dumper = UVWDumper(outname, *fields)
uvw_dumper << model
diff --git a/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py
index afa9b68..d0500b1 100644
--- a/python/tamaas/nonlinear_solvers/__init__.py
+++ b/python/tamaas/nonlinear_solvers/__init__.py
@@ -1,163 +1,163 @@
# -*- mode:python; coding: utf-8 -*-
#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+# Copyright (©) 2016-2022 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 .
"""
Pulling solvers to nonlinear_solvers module
"""
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, _, callback=None):
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 the nonlinear plasticity equation using the scipy_solve routine
"""
# 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
self._x[...] = self.scipy_solve(compute_residual)
# Computing displacements
self._residual.computeResidualDisplacement(self._x)
def reset(self):
"Set solution vector to zero"
self._x[...] = 0
class NewtonKrylovSolver(ScipySolver):
"""
Solve using a finite-difference Newton-Krylov method
"""
def __init__(self, residual, model=None, callback=None):
ScipySolver.__init__(self, residual, model,
callback=callback)
def scipy_solve(self, compute_residual):
"Solve R(delta epsilon) = 0 using a newton-krylov method"
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
"""
def __init__(self, residual, model=None, callback=None):
ScipySolver.__init__(self, residual, model,
callback=callback)
def scipy_solve(self, compute_residual):
"Solve R(delta epsilon) = 0 using a df-sane method"
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):
"Decorator to manage tolerance of non-linear solver"
# start /= rate # just anticipating first multiplication
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_module.cpp b/python/tamaas_module.cpp
index 0e88ec2..e2445b3 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,86 +1,86 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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 c5e8ace..c6ac47d 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,76 +1,76 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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 c323dd9..eb7ea6a 100644
--- a/python/wrap/compute.cpp
+++ b/python/wrap/compute.cpp
@@ -1,72 +1,72 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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) {
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/core.cpp b/python/wrap/core.cpp
index 083afd9..c0f27ec 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,116 +1,116 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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("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.");
}
/* -------------------------------------------------------------------------- */
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); });
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 714e3a6..c29c2b4 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,457 +1,457 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "adhesion_functional.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
/* -------------------------------------------------------------------------- */
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);
}
};
/// 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_ 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);
}
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");
}
}
/// 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);
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",
[](py::object) { return trait::dimension; },
"Dimension of computational domain")
.def_property_readonly_static(
"components", [](py::object) { return trait::components; },
"Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
[](py::object) { return trait::boundary_dimension; },
"Dimension of boundary of computational domain")
.def_property_readonly_static(
"voigt", [](py::object) { return trait::voigt; },
"Number of components of symmetrical tensor fields")
.def_property_readonly_static("indices",
[](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_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("registerField",
[](Model& m, std::string name, numpy field) {
TAMAAS_DEPRECATE("registerField()", "the [] operator");
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def("getField",
[](const Model& m, std::string name) -> decltype(m.getField(name)) {
TAMAAS_DEPRECATE("getField()", "the [] operator");
return m.getField(name);
},
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields",
[](const Model& m) {
TAMAAS_DEPRECATE("getFields()", "list(model)");
return m.getFields();
},
"Return fields list")
.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("__getitem__",
[](const Model& m, std::string key) -> decltype(m[key]) {
try {
return m[key];
} catch (std::out_of_range&) {
throw py::key_error(key);
}
},
py::return_value_policy::reference_internal, "Get field")
.def("__setitem__",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
py::keep_alive<1, 3>(), "Register new field")
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), key) !=
fields.end();
},
py::keep_alive<0, 1>(), "Test field existence")
.def("__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>(), "Iterator on fields")
.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", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "global_discretization"_a,
"Create a new model of a given type, physical size and "
"*global* discretization")
.def_static("createResidual", &ModelFactory::createResidual, "model"_a,
"sigma_y"_a, "hardening"_a = 0.,
"Create an isotropic linear hardening residual")
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a,
"Register Boussinesq and Mindlin operators to model");
}
/// 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) {
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 5cc5972..9ff457d 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,147 +1,147 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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 7ebb7fb..806f90e 100644
--- a/python/wrap/mpi.cpp
+++ b/python/wrap/mpi.cpp
@@ -1,92 +1,92 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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("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 3fa6c13..210f17c 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,92 +1,92 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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(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 3f6d67d..ad74fc2 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,216 +1,216 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_(mod, "ContactSolver")
.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("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("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 8967ff6..07ec191 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,177 +1,177 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Copyright (©) 2016-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#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, computeFilter,
filter_coefficients);
}
};
template
void wrapFilter(py::module& mod) {
auto name = makeDimensionName("Filter", dim);
py::class_, std::shared_ptr>, PyFilter>(
mod, name.c_str(), "Mother class for Fourier filter objects")
.def(py::init<>())
.def("computeFilter",
(void (Filter::*)(GridHermitian&) const) &
Filter::computeFilter,
"Compute the Fourier coefficient of the surface");
}
/* -------------------------------------------------------------------------- */
template
void wrapIsopowerlaw(py::module& mod) {
std::string name = makeDimensionName("Isopowerlaw", dim);
py::class_, Filter, std::shared_ptr>>(
mod, name.c_str(), "Isotropic powerlaw spectrum with a rolloff plateau")
.def(py::init<>())
.def_property("q0", &Isopowerlaw::getQ0, &Isopowerlaw::setQ0,
"Long wavelength cutoff")
.def_property("q1", &Isopowerlaw::getQ1, &Isopowerlaw::setQ1,
"Rolloff wavelength")
.def_property("q2", &Isopowerlaw::getQ2, &Isopowerlaw::setQ2,
"Short wavelength cutoff")
.def_property("hurst", &Isopowerlaw::getHurst,
&Isopowerlaw::setHurst, "Hurst exponent")
.def("rmsHeights", &Isopowerlaw::rmsHeights,
"Theoretical RMS of heights")
.def("moments", &Isopowerlaw::moments,
"Theoretical first 3 moments of spectrum")
.def("alpha", &Isopowerlaw::alpha, "Nayak's bandwidth parameter")
.def("rmsSlopes", &Isopowerlaw::rmsSlopes,
"Theoretical RMS of slopes");
name = makeDimensionName("RegularizedPowerlaw", dim);
py::class_, Filter,
std::shared_ptr>>(
mod, name.c_str(),
"Isotropic regularized powerlaw with a plateau extending to the size of "
"the system")
.def(py::init<>())
.def_property("q1", &RegularizedPowerlaw::getQ1,
&RegularizedPowerlaw::setQ1, "Rolloff wavelength")
.def_property("q2", &RegularizedPowerlaw::getQ2,
&RegularizedPowerlaw::setQ2, "Short wavelength cutoff")
.def_property("hurst", &RegularizedPowerlaw::getHurst,
&RegularizedPowerlaw::setHurst, "Hurst exponent");
}
/* -------------------------------------------------------------------------- */
template
void wrapSurfaceGenerators(py::module& mod) {
std::string generator_name = makeDimensionName("SurfaceGenerator", dim);
py::class_>(mod, generator_name.c_str())
.def("buildSurface", &SurfaceGenerator::buildSurface,
py::return_value_policy::reference_internal,
"Generate a surface and return a reference to it")
.def("setSizes",
[](SurfaceGenerator& m, std::array s) {
TAMAAS_DEPRECATE("setSizes()", "the shape property");
m.setSizes(std::move(s));
})
.def("setRandomSeed",
[](SurfaceGenerator& m, long s) {
TAMAAS_DEPRECATE("setRandomSeed()", "the random_seed property");
m.setRandomSeed(s);
})
.def_property("random_seed", &SurfaceGenerator::getRandomSeed,
&SurfaceGenerator::setRandomSeed,
"Random generator seed")
.def_property("shape", &SurfaceGenerator::getSizes,
py::overload_cast>(
&SurfaceGenerator::setSizes),
"Global shape of surfaces");
std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim);
py::class_, SurfaceGenerator>(
mod, filter_name.c_str(),
"Generates a rough surface with Gaussian noise in the PSD")
.def(py::init<>(), "Default constructor")
.def(py::init>(),
"Initialize with global surface shape")
.def("setFilter",
[](SurfaceGeneratorFilter