diff --git a/Jenkinsfile b/Jenkinsfile
index d3b08fe..3293b72 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -1,120 +1,119 @@
pipeline {
parameters {string(defaultValue: '', description: 'api-token', name: 'API_TOKEN')
string(defaultValue: '', description: 'Token for readthedocs', name: 'RTD_TOKEN')
string(defaultValue: '', description: 'buildable phid', name: 'BUILD_TARGET_PHID')
string(defaultValue: '', description: 'Commit id', name: 'COMMIT_ID')
string(defaultValue: '', description: 'Diff id', name: 'DIFF_ID')
string(defaultValue: 'PHID-PROJ-gbo56hpf2y5bi7t5jusk', description: 'ID of the project', name: 'PROJECT_ID')
}
environment {
PHABRICATOR_HOST = 'https://c4science.ch/api/'
PYTHONPATH = sh returnStdout: true, script: 'echo ${WORKSPACE}/tests/ci/script/'
}
agent {
dockerfile { additionalBuildArgs '--tag tamaas-environment'}
}
stages {
stage('SCM Checkout') {
steps {
checkout scm: [
$class: 'GitSCM',
branches: scm.branches,
extensions: [[
$class: 'SubmoduleOption',
recursiveSubmodules: true,
]],
userRemoteConfigs: scm.userRemoteConfigs
]
}
}
stage('Configure') {
steps {
sh '''#!/usr/bin/env bash
echo "py_exec = \'python3\'" > build-setup.conf
echo "build_python = \'true\'" >> build-setup.conf
echo "build_tests = \'true\'" >> build-setup.conf
echo "use_googletest = \'true\'" >> build-setup.conf
- echo "legacy_bem = \'true\'" >> build-setup.conf
echo "use_mpi = \'true\'" >> build-setup.conf'''
}
}
stage('Compile') {
steps {
sh '''#!/usr/bin/env bash
set -o pipefail
rm -rf .sconf_temp .sconsign.dblite build-release
scons | tee compilation.txt'''
}
post {
failure {
uploadArtifact('compilation.txt', 'Compilation')
}
}
}
stage('Run tests') {
steps {
sh 'PYTHONPATH=build-release/python/ python3 -m pytest --durations=0 --junitxml=results.xml build-release'
}
}
}
post {
always {
createArtifact("results.xml")
junit 'results.xml'
step([$class: 'XUnitBuilder',
thresholds: [
[$class: 'SkippedThreshold', failureThreshold: '0'],
[$class: 'FailedThreshold', failureThreshold: '0']],
tools: [[$class: 'GoogleTestType', pattern: 'build-release/tests/gtest_results.xml']]])
}
success {
trigger_rtd()
passed_hbm()
}
failure {
failed_hbm()
emailext(
body: '''${SCRIPT, template="groovy-html.template"}''',
mimeType: 'text/html',
subject: "[Jenkins] ${currentBuild.fullDisplayName} Failed",
recipientProviders: [[$class: 'CulpritsRecipientProvider']],
to: 'lucas.frerot@protonmail.com',
attachLog: true,
compressLog: true)
}
}
}
def createArtifact(filename) {
sh "./tests/ci/scripts/hbm send-uri -k 'Jenkins URI' -u ${BUILD_URL} -l 'View Jenkins result'"
sh "./tests/ci/scripts/hbm send-junit-results -f ${filename}"
}
def uploadArtifact(artifact, name) {
sh "./test/ci/scripts/hbm upload-file -f ${artifact} -n \"${name}\" -v ${PROJECT_ID}"
}
def failed_hbm() {
sh "./tests/ci/scripts/hbm failed"
}
def passed_hbm() {
sh "./tests/ci/scripts/hbm passed"
}
def trigger_rtd() {
sh """
set -x
curl -X POST -d "token=${RTD_TOKEN}" https://readthedocs.org/api/v2/webhook/tamaas/106141/
"""
}
diff --git a/SConstruct b/SConstruct
index 3c025be..2a560b7 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,477 +1,469 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
from __future__ import print_function
import os
from subprocess import check_output
import SCons
from os.path import join, abspath, basename
# Import below not strictly necessary, but good for pep8
from SCons.Script import (
EnsurePythonVersion,
EnsureSConsVersion,
Help,
Environment,
Variables,
EnumVariable,
PathVariable,
BoolVariable,
Split,
SConscript,
Export,
Dir
)
from version import get_git_subst
from detect import (
FindFFTW,
FindBoost,
FindThrust,
FindCuda,
FindExpolit,
FindPybind11
)
# ------------------------------------------------------------------------------
EnsurePythonVersion(2, 7)
EnsureSConsVersion(2, 4)
# ------------------------------------------------------------------------------
tamaas_info = dict(
version="2.2.0",
authors=[
u'Lucas Frérot',
'Guillaume Anciaux',
'Valentine Rey',
'Son Pham-Ba',
u'Jean-François Molinari'
],
maintainer=u'Lucas Frérot',
email='lucas.frerot@protonmail.com',
copyright=u"Copyright (©) 2016-2021 EPFL "
+ u"(École Polytechnique Fédérale de Lausanne), "
+ u"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
# ------------------------------------------------------------------------------
def detect_dependencies(env):
"""Detect all dependencies"""
fftw_comp = {
'omp': ['omp'],
'cpp': [],
'tbb': ['threads']
}
fftw_components = fftw_comp[env['backend']]
if main_env['use_mpi']:
fftw_components.append('mpi')
FindFFTW(env, fftw_components, precision=env['real_type'])
FindBoost(env, ['boost/preprocessor/seq.hpp'])
FindExpolit(env)
thrust_var = 'THRUST_ROOT'
# Take cuda version of thrust if available
if 'CUDA_ROOT' in env['ENV']:
thrust_var = 'CUDA_ROOT'
FindThrust(env, env['backend'], thrust_var)
# Activate cuda if needed
if env['backend'] == 'cuda':
FindCuda(env)
if env['build_python']:
FindPybind11(env)
# ------------------------------------------------------------------------------
# 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_info.items():
main_env[k] = v
main_env['COLOR_DICT'] = colors
# 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'),
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', "")),
- # Compile legacy code
- BoolVariable('legacy_bem', 'Compile legacy BEM code', False),
# Cosmetic
BoolVariable('verbose', 'Activate verbosity', False),
BoolVariable('color', 'Color the non-verbose compilation output', False),
# Tamaas components
BoolVariable('build_doc', 'Build documentation', False),
BoolVariable('build_tests', 'Build test suite', False),
BoolVariable('build_python', 'Build python wrapper', True),
# Dependencies
BoolVariable('use_googletest', 'Build tests using GTest', False),
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 + "\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-' + main_env['build_type']
main_env['build_dir'] = build_dir
if main_env['should_configure']:
print("Building in " + build_dir)
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
# Include paths
main_env.AppendUnique(CPPPATH=['#/src',
'#/src/core',
'#/src/mpi',
'#/src/bem',
'#/src/surface',
'#/src/python',
'#/src/percolation',
'#/src/model',
'#/src/model/elasto_plastic',
'#/src/solvers',
'#/src/gpu',
'#/python'])
# Changing the shared object extension
main_env['SHOBJSUFFIX'] = '.o'
# Back to gcc if cuda is activated
if main_env['backend'] == "cuda" and "g++" not in main_env['CXX']:
raise SCons.Errors.StopError('GCC should be used when compiling with CUDA')
# 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 SCons.Errors.StopError('Unsupported compiler: ' + cxx)
cxx = cxx_alias(main_env['CXX'])
# Setting main compilation flags
main_env['CXXFLAGS'] = Split(main_env['CXXFLAGS'])
main_env.AppendUnique(
CXXFLAGS=Split('-std=c++14 -Wall -Wextra -pedantic'),
CPPDEFINES=['TAMAAS_BACKEND=TAMAAS_BACKEND_{}'.format(
main_env['backend'].upper()
)],
)
# 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'])
-# Force deactivate legacy when using intel
-# Intel does not like old openmp loops
-if cxx == 'icpc' and main_env['legacy_bem']:
- print("Using intel compiler => deactivating legacy code")
- main_env['legacy_bem'] = False
-
# 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 SCons.Errors.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']:
detect_dependencies(main_env)
# Setting up the python name with version
if main_env['build_python']:
args = (main_env.subst("${py_exec} -c").split()
+ ["from distutils.sysconfig import get_python_version;"
"print(get_python_version())"])
main_env['py_version'] = bytes(check_output(args)).decode()
# 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@': abspath(build_dir),
'@build_version@': '$version',
})
# 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(join(build_dir, 'tamaas_environment.sh'),
env_content)
# Building sub-directories
def subdir(dir):
return SConscript(join(dir, 'SConscript'),
variant_dir=join(build_dir, dir),
duplicate=True)
# Building Tamaas library
Export('main_env')
subdir('src')
build_targets = ['build-cpp', env_file]
install_targets = ['install-lib']
# Building Tamaas extra components
for dir in ['python', 'tests']:
if main_env['build_{}'.format(dir)] and not main_env.GetOption('help'):
subdir(dir)
build_targets.append('build-{}'.format(dir))
# Building API + Sphinx documentation if requested
if main_env['build_doc']:
doc_env = main_env.Clone()
doc = doc_env.Command('#.phony_doc', '',
'make -C {doc} clean && make -C {doc}'
.format(doc=Dir('#/doc')))
if doc_env['build_python']:
doc_env.PrependENVPath('PYTHONPATH',
doc_env.subst('../${build_dir}/python'))
doc_env.Depends(doc, 'build-python')
main_env.Alias('doc', doc)
else:
dummy_command(main_env, 'doc', 'Command "doc" does not do anything'
' without documentation activated ("build_doc=True")')
# 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-{}.tar.gz'.format(main_env['version']),
'',
('tar --exclude-vcs --exclude-vcs-ignores '
'--exclude=third-party/googletest '
'--exclude=third-party/pybind11 '
'--exclude=joss '
'--exclude=".*" '
'-czf $TARGET {}'.format(basename(os.getcwd()))),
chdir='..',
)
main_env.Alias('archive', archive)
diff --git a/doc/sphinx/source/contact.rst b/doc/sphinx/source/contact.rst
index c72b05d..70dde2f 100644
--- a/doc/sphinx/source/contact.rst
+++ b/doc/sphinx/source/contact.rst
@@ -1,137 +1,137 @@
Solving contact
===============
The resolution of contact problems is done with classes that inherit from :cpp:class:`tamaas::ContactSolver`. These usually take as argument a :cpp:class:`tamaas::Model` object, a surface described by a :cpp:class:`tamaas::Grid` or a 2D numpy array, and a tolerance. We will see the specificities of the different solver objects below.
Elastic contact
---------------
The most common case is normal elastic contact, and is most efficiently solved with :cpp:class:`tamaas::PolonskyKeerRey`. The advantage of this class is that it combines two algorithms into one. By default, it considers that the contact pressure field is the unknown, and tries to minimize the complementary energy of the system under the constraint that the mean pressure should be equal to the value supplied by the user, for example::
# ...
solver = tm.PolonskyKeerRey(model, surface, 1e-12)
solver.solve(1e-2)
Here the average pressure is ``1e-2``. The solver can also be instanciated by specifying the the constraint should be on the mean gap instead of mean pressure::
solver = tm.PolonskyKeerRey(model, surface, 1e-12, constraint_type=tm.PolonskyKeerRey.gap)
solver.solve(1e-2)
The contact problem is now solved for a mean gap of ``1e-2``. Note that the choice of constraint affects the performance of the algorithm.
Contact with adhesion
---------------------
The second algorithm hidden in :cpp:class:`tamaas::PolonskyKeerRey` considers the **gap** to be the unknown. The constaint on the mean value can be chosen as above::
solver = tm.PolonskyKeerRey(model, surface, 1e-12,
primal_type=tm.PolonskyKeerRey.gap,
constraint_type=tm.PolonskyKeerRey.gap)
solver.solve(1e-2)
The advantage of this formulation is to be able to solve adhesion problems (`Rey et al. `_). This is done by adding a term to the potential energy functional that the solver tries to minimize::
adhesion_params = {
"rho": 2e-3,
"surface_energy": 2e-5
}
adhesion = tm.ExponentialAdhesionFunctional(surface)
adhesion.setParameters(adhesion)
solver.addFunctionalterm(adhesion)
solver.solve(1e-2)
Custom classes can be used in place of the example term here. One has to inherit from :cpp:class:`tamaas::Functional`::
class AdhesionPython(tm.Functional):
"""
Functional class that extends a C++ class and implements the virtual
methods
"""
def __init__(self, rho, gamma):
super().__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
This example is actually equivalent to :cpp:class:`tamaas::functional::ExponentialAdhesionFunctional`.
Tangential contact
------------------
For frictional contact, several solver classes are available. Among them, :cpp:class:`tamaas::Condat` is able to solve a coupled normal/tangential contact problem regardless of the material properties. It however solves an associated version of the Coulomb friction law. In general, the Coulomb friction used in contact makes the problem ill-posed.
Solving a tangential contact problem is not much different from normal contact::
mu = 0.3 # Friction coefficient
solver = tm.Condat(model, surface, 1e-12, mu)
- solver.setMaxIterations(5000) # The default of 1000 may be too little
+ solver.max_iter = 5000 # The default of 1000 may be too little
solver.solve([1e-2, 0, 1e-2]) # 3D components of applied mean pressure
Elasto-plastic contact
----------------------
For elastic-plastic contact, one needs three different solvers: an elastic contact solver like the ones described above, a non-linear solver and a coupling solver. The non-linear solvers available in Tamaas are implemented in python and inherit from the C++ class :cpp:class:`tamaas::EPSolver`. They make use of the non-linear solvers available in scipy:
:py:class:`DFSANESolver `
Implements an interface to :py:func:`scipy.optimize.root` wiht the DFSANE method.
:py:class:`NewtonKrylovSolver `
Implements an interface to :py:func:`scipy.optimize.newton_krylov`.
These solvers require a residual vector to cancel, the creation of which is handled with :cpp:class:`tamaas::ModelFactory`. Then, an instance of :cpp:class:`tamaas::EPICSolver` is needed to couple the contact and non-linear solvers for the elastic-plastic contact problem::
import tamaas as tm
from tamaas.nonlinear_solvers import DFSANESolver
# Definition of modeled domain
model_type = tm.model_type.volume_2d
discretization = [32, 51, 51] # Order: [z, x, y]
flat_domain = [1, 1]
system_size = [0.5] + flat_domain
# Setup for plasticity
residual = tm.ModelFactory.createResidual(model,
sigma_y=0.1 * model.E,
hardening=0.01 * model.E)
epsolver = DFSANESolver(residual)
# ...
csolver = tm.PolonskyKeerRey(model, surface, 1e-12)
epic = tm.EPICSolver(csolver, epsolver, 1e-7, relaxation=0.3)
epic.solve(1e-3)
# or use an accelerated scheme (which can sometimes not converge)
epic.acceleratedSolve(1e-3)
By default, :cpp:func:`tamaas::EPICSolver::solve` uses a relaxed fixed point. It can be tricky to make it converge: you need to decrease the relaxation parameter passed as argument of the constructor, but this also hinders the convergence rate. The function :cpp:func:`tamaas::EPICSolver::acceleratedSolve` does not require the tweaking of a relaxation parameter, so it can be faster if the latter does not have an optimal value. However, it is not guaranteed to converge.
Finally, during the first iterations, the fixed point error will be large compared to the error of the non-linear solver. It can improve performance to have the tolerance of the non-linear solver (which is the most expensive step of the fixed point solve) decrease over the iterations. This can be achieved with the decorator class :py:class:`ToleranceManager `::
from tamaas.nonlinear_solvers import ToleranceManager, DFSANESolver
@ToleranceManager(start=1e-2,
end=1e-9,
rate=1 / 3)
class Solver(DFSANESolver):
pass
# ...
epsolver = Solver(residual)
# or
epsolver = ToleranceManager(1e-2, 1e-9, 1/3)(DFSANESolver)(residual)
diff --git a/doc/sphinx/source/model.rst b/doc/sphinx/source/model.rst
index 4be5543..9ffbd20 100644
--- a/doc/sphinx/source/model.rst
+++ b/doc/sphinx/source/model.rst
@@ -1,103 +1,156 @@
Model and integral operators
============================
The class :cpp:class:`tamaas::Model` (and its counterpart :py:class:`Model `) is both a central class in Tamaas and one of the simplest. It mostly serves as holder for material properties, fields and integral operators, and apart from a linear elastic behavior does not perform any computation on its own.
Model types
-----------
:cpp:class:`tamaas::Model` has a concrete subclass :cpp:class:`tamaas::ModelTemplate` which implements the model function for a given :cpp:type:`tamaas::model_type`:
:cpp:enumerator:`tamaas::basic_2d`
Model type used in normal frictionless contact: traction and displacement are 2D fields with only one component.
:cpp:enumerator:`tamaas::surface_2d`
Model type used in frictional contact: traction and displacement are 2D fields with three components.
:cpp:enumerator:`tamaas::volume_2d`
Model type used in elastoplastic contact: tractions are the same as with :cpp:enumerator:`tamaas::surface_2d` but the displacement is a 3D field.
The enumeration values suffixed with ``_1d`` are the one dimensional (line contact) counterparts of the above model types. The domain physical dimension and number of components are encoded in the class :cpp:class:`tamaas::model_type_traits`.
Model creation and basic functionality
--------------------------------------
The instanciation of a :py:class:`Model ` is done with the :py:class:`ModelFactory ` class and its :py:func:`createModel ` function::
physical_size = [1., 1.]
discretization = [512, 512]
model = tm.ModelFactory.createModel(tm.model_type.basic_2d, physical_size, discretization)
.. warning::
For models of type ``volume_*d``, the first component of the ``physical_size`` and ``discretization`` arrays corresponds to the depth dimension (:math:`z` in most cases). For example::
tm.ModelFactory.createModel(tm.model_type.basic_2d, [0.3, 1, 1], [64, 81, 81])
creates a model of depth 0.3 and surface size 1\ :superscript:`2`, while the number of points is 64 in depth and 81\ :sup:`2` on the surface. This is done for data contiguity reasons, as we do discrete Fourier transforms in the horizontal plane.
.. note::
If ran in an MPI context, the method :py:meth:`createModel
` expects the *global* system sizes
and discretization of the model.
The properties ``E`` and ``nu`` can be used to set the Young's modulus and Poisson ratio respectively::
model.E = 1
model.nu = 0.3
Fields can be easlily accessed with the ``[]`` operator, similar to Python's dictionaries::
surface_traction = model['traction']
+ # surface_traction is a numpy array
-To know what fields are available, you can call the :py:func:`getFields ` method. You can add new fields to a model object with :py:func:`registerField `, which is convenient for dumping. A model can also be used to compute stresses from a strain field::
+To know what fields are available, you can call the :py:class:`list` function on
+a model (``list(model)``). You can add new fields to a model object with the
+``[]`` operator:``model['new_field'] = new_field``, which is convenient for
+dumping fields that are computed outside of Tamaas.
+
+.. note::
+ The fields ``traction`` and ``displacement`` are always registered in models,
+ and are accessible via :py:attr:`model.traction
+ ` and :py:attr:`model.displacement
+ `.
+
+A model can also be used to compute stresses from a strain field::
import numpy as np
- strain = np.zeros(model.getDiscretization() + [6]) # Mandel--Voigt notation
+ strain = np.zeros(model.shape + [6]) # Mandel--Voigt notation
stress = np.zeros_like(strain)
model.applyElasticity(stress, strain)
+.. tip::
+ ``print(model)`` gives a lot of information about the model: the model type,
+ shape, registered fields, and more!
+
+
+Integral operators
+------------------
+
+Integral operators are a central part of Tamaas: they are carefully designed for
+performance in periodic system. When a :py:class:`Model `
+object is used with contact solvers or with a residual object (for plasticty),
+the objects using the model register integral operators with the model, so the
+user typically does not have to worry about creating integral operators.
+
+Integral operators are accessed through the :py:attr:`operators
+` property of a model object. The ``[]``
+operator gives access to the operators, and ``list(model.operators)`` gives the
+list of registered operators::
+
+ # Accessing operator
+ elasticity = model.operators['hooke']
+
+ # Applying operator
+ elasticity(strain, stress)
+
+ # Print all registered operators
+ print(list(model.operators))
+
+.. note::
+ At model creation, these operators are automatically registered:
+
+ - ``hooke``: Hooke's elasticity law
+ - ``von_mises``: computes Von Mises stresses
+ - ``deviatoric``: computes the deviatoric part of a stress tensor
+ - ``eigenvalues``: computes the eigenvalues of a symetric tensor field
+
+ :cpp:class:`Westergaard ` operators are automatically
+ registered when :py:meth:`solveNeumann
+ ` or :py:meth:`solveDirichlet
+ ` are called.
+
Model dumpers
-------------
The submodule `tamaas.dumpers` contains a number of classes to save model data into different formats:
:py:class:`UVWDumper `
Dumps a model to `VTK `_ format. Requires the `UVW `_ python package which you can install with pip::
pip install uvw
This dumper is made for visualization with VTK based software like `Paraview `_.
:py:class:`NumpyDumper `
Dumps a model to a compressed Numpy file.
:py:class:`H5Dumper `
Dumps a model to a compressed `HDF5 `_ file. Requires the `h5py `_ package.
The dumpers are initialzed with a basename and the fields that you wish to write to file (optionally you can set ``all_fields`` to ``True`` to dump all fields in the model). By default, each write operation creates a new file in a separate directory (e.g. :py:class:`UVWDumper ` creates a ``paraview`` directory). To write to a specific file you can use the `dump_to_file` method. Here is a usage example::
from tamaas.dumpers import UVWDumper, H5Dumper
# Create dumper
uvw_dumper = UVWDumper('rough_contact_example', 'stress', 'plastic_strain')
# Dump model
uvw_dumper << model
# Or alternatively
model.addDumper(H5Dumper('rough_contact_archive', all_fields=True))
model.addDumper(uvw_dumper)
model.dump()
The last ``model.dump()`` call will call both dumpers. The resulting files will have the following hierachy::
./paraview/rough_contact_example_0000.vtr
./paraview/rough_contact_example_0001.vtr
./hdf5/rough_contact_archive_0000.h5
-
+
.. note::
Currently, only :py:class:`H5Dumper ` supports
parallel output with MPI.
diff --git a/doc/sphinx/source/rough_surfaces.rst b/doc/sphinx/source/rough_surfaces.rst
index 7a4a802..fd82e80 100644
--- a/doc/sphinx/source/rough_surfaces.rst
+++ b/doc/sphinx/source/rough_surfaces.rst
@@ -1,109 +1,111 @@
Random rough surfaces
=====================
The generation of stochatsticly rough surfaces is controlled in Tamaas by two abstract classes: :cpp:class:`tamaas::SurfaceGenerator` and :cpp:class:`tamaas::Filter`. The former provides access lets the user set the surface sizes and random seed, while the latter encodes the information of the spectrum of the surface. Two surface generation methods are provided:
- :cpp:class:`tamaas::SurfaceGeneratorFilter` implements a Fourier filtering algorithm (see `Hu & Tonder `_),
- :cpp:class:`tamaas::SurfaceGeneratorRandomPhase` implements a random phase filter.
Both of these rely on a :cpp:class:`tamaas::Filter` object to provided the filtering information (usually power spectrum density coefficients). Tamaas provides two such classes and allows for Python subclassing:
- :cpp:class:`tamaas::Isopowerlaw` provides a roll-off powerlaw,
- :cpp:class:`tamaas::RegularizedPowerlaw` provides a powerlaw with a regularized rolloff.
Tamaas also provided routines for surface statistics.
Generating a surface
--------------------
Let us now see how to generate a surface. Frist create a filter object and set the surface sizes::
import tamaas as tm
# Create spectrum object
spectrum = tm.Isopowerlaw2D()
# Set spectrum parameters
spectrum.q0 = 4
spectrum.q1 = 4
spectrum.q2 = 32
spectrum.hurst = 0.8
The ``spectrum`` object can be queried for information, such as the root-mean-square of heights, the various statistical moments, the spectrum bandwidth, etc. Then we create a generator object and build the random surface::
- generator = tm.SurfaceGeneratorFilter2D()
- generator.setSizes([128, 128])
- generator.setSpectrum(spectrum)
+ generator = tm.SurfaceGeneratorFilter2D([128, 128])
+ generator.spectrum = spectrum
generator.random_seed = 0
surface = generator.buildSurface()
.. important::
The ``surface`` object is a :py:class:`numpy.ndarray` wrapped
around internal memory in the ``generator`` object, so a subsequent call to
:py:func:`buildSurface ` may
change its content. Note that if ``generator`` goes out of scope its memory
will not be freed if there is still a live reference to the surface data.
-.. note::
- If ran in an MPI context, the method :py:meth:`setSizes
- ` expects the *global* shape of
- the surface.
+.. important::
+ If ran in an MPI context, the constructor of
+ :py:class:`SurfaceGeneratorFilter2D
+ ` (and others) expects the *global*
+ shape of the surface. The shape can also be changed with ``generator.shape =
+ [64, 64]``.
Custom filter
-------------
Tamaas provides several classes that can be derived directly with Python classes, and :cpp:class:`tamaas::Filter` is one of them. Since it provides a single pure virtual method :cpp:func:`computeFilter `, it is easy to write a sub-class. Here we implement a class that takes a user-defined auto-correlation function and implements the :cpp:func:`computeFilter ` virtual function::
class AutocorrelationFilter(tm.Filter2D):
def __init__(self, autocorrelation):
tm.Filter2D.__init__(self)
self.autocorrelation = autocorrelation.copy()
def computeFilter(self, filter_coefficients):
shifted_ac = np.fft.ifftshift(self.autocorrelation)
filter_coefficients[...] = np.sqrt(np.fft.rfft2(shifted_ac))
Here ``filter_coefficients`` is also a :py:class:`numpy.ndarray` and is therefore easily manipulated. The creation of the surface then follows the same pattern as previously::
# Create spectrum object
autocorrelation = ... # set your desired autocorrelation
spectrum = AutocorrelationFilter(autocorrelation)
generator = tm.SurfaceGenerator2D()
- generator.setSpectrum(spectrum)
+ generator.shape = autocorrelation.shape
+ generator.spectrum = spectrum
surface = generator.buildSurface()
The lifetime of the ``spectrum`` object is associated to the ``generator`` when :py:func:`setSpectrum ` is called.
Surface Satistics
-----------------
Tamaas provides the C++ class :cpp:class:`tamaas::Statistics` and its wrapper :py:class:`Statistics2D ` to compute statistics on surfaces, including:
- power spectrum density
- autocorrelation
- spectrum moments
- root-mean-square of heights :math:`\sqrt{\langle h^2 \rangle}`
- root-mean-square of slopes (computed in Fourier domain) :math:`\sqrt{\langle |\nabla h|^2\rangle}`
All these quantities are computed in a discretization-independent manner: increasing the number of points in the surface should not drastically change the computed values (for a given spectrum). This allows to refine the discretization as much as possible to approximate a continuum. Note that the autocorrelation and PSD are fft-shifted. Here is a sample code plotting the PSD and autocorrelation::
psd = tm.Statistics2D.computePowerSpectrum(surface)
psd = np.fft.fftshift(psd, axes=0) # Shifting only once axis because of R2C transform
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
plt.imshow(psd.real, norm=LogNorm())
acf = tm.Statistics2D.computeAutocorrelation(surface)
acf = np.fft.fftshift(acf)
plt.figure()
plt.imshow(acf)
plt.show()
See ``examples/statistics.py`` for more usage examples of statistics.
diff --git a/examples/adhesion.py b/examples/adhesion.py
index dfc8ea4..fbc4973 100644
--- a/examples/adhesion.py
+++ b/examples/adhesion.py
@@ -1,127 +1,122 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
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()
-sg.setSizes([n, n])
-sg.setRandomSeed(0)
+sg = tm.SurfaceGeneratorFilter2D([n, n])
+sg.random_seed = 0
# Spectrum
-spectrum = tm.Isopowerlaw2D()
-sg.setFilter(spectrum)
+sg.spectrum = tm.Isopowerlaw2D()
# Parameters
-spectrum.q0 = 16
-spectrum.q1 = 16
-spectrum.q2 = 64
-spectrum.hurst = 0.8
+sg.spectrum.q0 = 16
+sg.spectrum.q1 = 16
+sg.spectrum.q2 = 64
+sg.spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
-# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
surface /= n
-# print(spectrum.rmsSlopes())
-# print(tm.SurfaceStatistics.computeRMSSlope(surface))
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.setParameters(adhesion_params)
solver.addFunctionalTerm(adhesion)
# Solve for target pressure
g_target = 5e-2
solver.solve(g_target)
-tractions = model.getTraction()
+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/legacy/adhesion.py b/examples/legacy/adhesion.py
deleted file mode 100644
index 64ed1d5..0000000
--- a/examples/legacy/adhesion.py
+++ /dev/null
@@ -1,161 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-import sys
-import tamaas as tm
-import numpy as np
-import matplotlib.pyplot as plt
-
-def plotSurface(surface):
- fig = plt.figure()
- ax = fig.add_subplot(111)
- img = ax.imshow(surface)
- #fig.colorbar(img)
-
-def constructHertzProfile(size, curvature):
- radius = 1. / curvature
- x = np.linspace(-0.5, 0.5, size)
- y = np.linspace(-0.5, 0.5, size)
- x, y = np.meshgrid(x, y)
- surface = np.sqrt(radius**2 - x**2 - y**2)
- surface -= surface.min()
- return surface.copy()
-
-def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
- x = np.linspace(-0.5, 0.5, size)
- y = np.linspace(-0.5, 0.5, size)
- x, y = np.meshgrid(x, y)
- disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - (x**2 + y**2))
- return disp.copy()
-
-def main():
- grid_size = 128
- curvature = 0.5
- effective_modulus = 1.
- load = 0.1
-
- surface_energy = 0
- rho = 2.071e-7
-
- surface = constructHertzProfile(grid_size, curvature)
- # SG = tm.SurfaceGeneratorFilterFFT()
- # SG.getGridSize().assign(grid_size)
- # SG.getHurst().assign(0.8)
- # SG.getRMS().assign(0.002);
- # SG.getQ0().assign(8);
- # SG.getQ1().assign(8);
- # SG.getQ2().assign(16);
- # SG.getRandomSeed().assign(156);
- # SG.Init()
- # surface = SG.buildSurface()
- print "Max height {}".format(surface.max())
- print "Min height {}".format(surface.min())
-
- bem = tm.BemGigi(surface)
- bem.setDumpFreq(1)
- functional = tm.ExponentialAdhesionFunctional(bem)
- functional.setParameter('rho', rho)
- functional.setParameter('surface_energy', surface_energy)
- bem.setEffectiveModulus(effective_modulus)
- bem.addFunctional(functional)
- bem.computeEquilibrium(1e-6, load)
-
- tractions = bem.getTractions()
- print "Average pressure = {}".format(tractions.mean())
- # bem.computeTrueDisplacements()
- t_displacements = bem.getTrueDisplacements()
- t_gap = bem.getGap()
-
- plotSurface(tractions)
-
- plt.figure()
- plt.plot(surface[grid_size/2, :])
- plt.title("Surface")
- plt.figure()
- plt.plot(tractions[grid_size/2, :])
- plt.title("Pressure")
- plt.figure()
- plt.plot(t_gap[grid_size/2, :])
- plt.title("Gap")
- plt.figure()
- plt.plot(t_displacements[grid_size/2, :])
- plt.title("Displacement")
- plt.figure()
- plt.plot(t_displacements[grid_size/2, :] - surface[grid_size/2, :])
- plt.title("Disp-surf")
- plotSurface(t_displacements)
-
- plt.show()
- return 0
-
- # Testing contact area against Hertz solution for solids of revolution
- contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
- hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
- hertz_area = np.pi * hertz_contact_size**2
- area_error = np.abs(hertz_area - contact_area) / hertz_area
- print "Area error: {}".format(area_error)
-
- # Testing maximum pressure
- max_pressure = tractions.max()
- hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
- pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
- print "Max pressure error: {}".format(pressure_error)
-
- # Testing displacements
- hertz_displacements = computeHertzDisplacement(effective_modulus,
- hertz_contact_size,
- hertz_max_pressure,
- grid_size)
-
- # Selecing only the points that are in contact
- contact_indexes = [(i, j, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)]
- contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes))
-
- # Displacements of bem are centered around the mean of the whole surface
- # and Hertz displacements are not centered, so we need to compute mean
- # on the contact zone for both arrays
- bem_mean = 0.
- hertz_mean = 0.
- for index in contact_indexes:
- bem_mean += displacements[index]
- hertz_mean += hertz_displacements[index]
-
- bem_mean /= len(contact_indexes)
- hertz_mean /= len(contact_indexes)
- # Correction applied when computing error
- correction = hertz_mean - bem_mean
-
- # Computation of error of displacement in contact zone
- error = 0.
- hertz_norm = 0.
- for index in contact_indexes:
- error += (hertz_displacements[index] - displacements[index] - correction)**2
- hertz_norm += (hertz_displacements[index] - hertz_mean)**2
-
- displacement_error = np.sqrt(error / hertz_norm)
- print "Displacement error (in contact zone): {}".format(displacement_error)
-
- if area_error > 1e-2 or pressure_error > 1e-2 or displacement_error > 1e-4:
- return 1
-
- return 0
-
-if __name__ == "__main__":
- sys.exit(main())
diff --git a/examples/legacy/cluster_statistics.py b/examples/legacy/cluster_statistics.py
deleted file mode 100644
index 8b1bca2..0000000
--- a/examples/legacy/cluster_statistics.py
+++ /dev/null
@@ -1,104 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-import numpy as np
-
-
-def plotSurface(surf):
- fig = plt.figure()
- axe = fig.add_subplot(111)
- img = axe.imshow(surf)
- cbar = fig.colorbar(img)
-
-
-################################################################
-# surface generation
-################################################################
-
-from tamaas import *
-
-n = 128
-SG = SurfaceGeneratorFilterFFT()
-SG.getGridSize().assign(n)
-SG.getHurst().assign(0.8)
-SG.getRMS().assign(1.);
-SG.getQ0().assign(4);
-SG.getQ1().assign(4);
-SG.getQ2().assign(64);
-SG.getRandomSeed().assign(20);
-SG.Init()
-s = SG.buildSurface()
-rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(s)
-s *= 1./rms_slopes_spectral
-s -= s.max()
-
-################################################################
-# surface plot
-################################################################
-
-import matplotlib.pyplot as plt
-
-#fig1 = plotSurface(s)
-
-################################################################
-# contact contact
-################################################################
-
-bem = BemPolonski(s)
-bem.setEffectiveModulus(1.)
-
-load = 0.1
-bem.computeEquilibrium(1e-13,load)
-
-tractions = bem.getTractions()
-displacements = bem.getDisplacements()
-
-#fig2 = plotSurface(tractions)
-#fig3 = plotSurface(displacements)
-
-################################################################
-# cluster statistics
-################################################################
-
-
-area = ContactArea(tractions.shape[0],1.);
-area.getSurface()[tractions > 0.] = 1
-# Cluster detection
-print "==============================="
-print "Detect contact clusters: "
-area.detectContactClusters()
-print "Found " + str(area.getNumClusters())
-clusters = ContactClusterCollection(area);
-cluster_list = area.getContactClusters();
-cluster_areas = [float(c.getA())/float(n)**2 for c in cluster_list]
-cluster_perimeters = [c.getP() for c in cluster_list]
-print "cluster_areas",cluster_areas
-print "cluster_perimeters",cluster_perimeters
-print "cluster_total_area",clusters.getTotalArea()
-print "cluster_total_perimeter",clusters.getTotalPerimeter()
-print "nb_clusters_with_hole",clusters.getNbClustersWithHole()
-print "nb_clusters",clusters.getNbClusters()
-################################################################
-value,bins = np.histogram(cluster_areas,bins=50)
-fig = plt.figure()
-axe = fig.add_subplot(111)
-axe.loglog(bins[:-1],value,'-o')
-################################################################
-plt.show()
diff --git a/examples/legacy/contact.py b/examples/legacy/contact.py
deleted file mode 100644
index 0adbc84..0000000
--- a/examples/legacy/contact.py
+++ /dev/null
@@ -1,73 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-
-def plotSurface(surf):
- fig = plt.figure()
- axe = fig.add_subplot(111)
- img = axe.imshow(surf)
- cbar = fig.colorbar(img)
-
-
-################################################################
-# surface generation
-################################################################
-
-from tamaas import *
-
-
-SG = SurfaceGeneratorFilterFFT()
-SG.getGridSize().assign(128)
-SG.getHurst().assign(0.8)
-SG.getRMS().assign(1.);
-SG.getQ0().assign(4);
-SG.getQ1().assign(4);
-SG.getQ2().assign(64);
-SG.getRandomSeed().assign(20);
-SG.Init()
-s = SG.buildSurface()
-rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(s)
-s *= 1./rms_slopes_spectral
-s -= s.max()
-
-################################################################
-# surface plot
-################################################################
-
-import matplotlib.pyplot as plt
-
-fig1 = plotSurface(s)
-
-################################################################
-# contact contact
-################################################################
-
-bem = BemPolonski(s)
-bem.setEffectiveModulus(1.)
-
-load = 0.1
-bem.computeEquilibrium(1e-13,load)
-
-tractions = bem.getTractions()
-displacements = bem.getDisplacements()
-
-fig2 = plotSurface(tractions)
-fig3 = plotSurface(displacements)
-plt.show()
diff --git a/examples/legacy/fractal_surface.py b/examples/legacy/fractal_surface.py
deleted file mode 100644
index 0a07a9c..0000000
--- a/examples/legacy/fractal_surface.py
+++ /dev/null
@@ -1,102 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-from tamaas import *
-import argparse
-
-parser = argparse.ArgumentParser()
-parser.add_argument("--rescale", help="Rescale surface for RMS(slopes) = 1", action='store_true')
-parser.add_argument("--N", help="Number of points", type=int, default=512)
-parser.add_argument("--k0", help="Roll-off wave number", type=int, default=4)
-parser.add_argument("--k1", help="Low cutoff wave number", type=int, default=4)
-parser.add_argument("--k2", help="High cutoff wave number", type=int, default=32)
-parser.add_argument("--rms", help="RMS(heights)", default=1.)
-parser.add_argument("--H", help="Hurst exponent", default=0.8)
-
-args = parser.parse_args()
-
-#generate surface
-SG = SurfaceGeneratorFilterFFT()
-SG.getGridSize().assign(args.N)
-SG.getHurst().assign(args.H)
-SG.getRMS().assign(args.rms);
-SG.getQ0().assign(args.k0);
-SG.getQ1().assign(args.k1);
-SG.getQ2().assign(args.k2);
-SG.getRandomSeed().assign(156);
-SG.Init()
-a = SG.buildSurface()
-
-if args.rescale:
- rms_slopes = SurfaceStatistics.computeSpectralRMSSlope(a)
- a /= rms_slopes
-
-#compute and print surface statistics
-class Stats:
-
- def __getitem__(self,key):
- return self.__dict__[key]
-
-
-stats = Stats()
-
-stats.size = SG.getGridSize()
-stats.hurst = SG.getHurst().value()
-stats.rms = SG.getRMS()
-stats.k0 = SG.getQ0()
-stats.k1 = SG.getQ1().value()
-stats.k2 = SG.getQ2().value()
-stats.seed = SG.getRandomSeed()
-stats.rms_spectral = SurfaceStatistics.computeSpectralStdev(a);
-stats.rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(a);
-stats.rms_geometric = a.std(ddof=1)
-stats.rms_slopes_geometric = SurfaceStatistics.computeRMSSlope(a);
-stats.moments = SurfaceStatistics.computeMoments(a);
-stats.m0 = stats['rms_spectral']**2
-stats.m2 = stats.moments[0]
-stats.m4 = stats.moments[1]
-stats.alpha = stats['m0']*stats['m4']/(stats['m2']**2)
-stats.L = 1.
-stats.m0prime = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,1. , stats.L)
-stats.moment_A = stats.m0/stats.m0prime
-stats.analytic_m0 = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
-stats.analytic_m2 = SurfaceStatistics.computeAnalyticFractalMoment(2,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
-stats.analytic_m4 = SurfaceStatistics.computeAnalyticFractalMoment(4,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
-stats.analytic_alpha = stats.analytic_m0*stats.analytic_m4/(stats.analytic_m2*stats.analytic_m2);
-
-print """
-[N] {size}
-[rms] {rms}
-[rmsSpectral] {rms_spectral}
-[rmsSlopeSpectral] {rms_slopes_spectral}
-[rmsSlopeGeometric] {rms_slopes_geometric}
-[Hurst] {hurst}
-[k1] {k1}
-[k2] {k2}
-[moment A] {moment_A}
-[m0] {m0}
-[analytic m0] {analytic_m0}
-[m2] {m2}
-[analytic m2] {analytic_m2}
-[m4] {m4}
-[analytic m4] {analytic_m4}
-[alpha] {alpha}
-[analytic_alpha] {analytic_alpha}
-[seed] {seed}
diff --git a/examples/legacy/hertz.py b/examples/legacy/hertz.py
deleted file mode 100644
index 2a30094..0000000
--- a/examples/legacy/hertz.py
+++ /dev/null
@@ -1,116 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-import tamaas
-
-
-def dumpHertzPrediction(file, bem, pressure, r, E ):
- A = tamaas.SurfaceStatistics.computeContactArea(bem.getTractions())
- Aratio = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
- radius = np.sqrt(A/np.pi)
- #L = bem.getSurface().getL()
- L = 1.
- load = pressure * L*L
- radius_hertz = pow(0.75*load*r/E,1./3.)
- p0_hertz = 1./np.pi*pow(6.*E*E*load/r/r,1./3.)
- p0 = tamaas.SurfaceStatistics.computeMaximum(bem.getTractions())
- n = bem.getTractions().shape[0]
- computed_load = bem.getTractions().sum()/n/n*L*L
- file.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\t{9}\n".format(pressure,load,computed_load,Aratio,A,radius,radius_hertz,radius_hertz/radius,p0,p0_hertz))
-
-
-################################################################
-
-def main(argv):
-
- parser = argparse.ArgumentParser(description='Hertz test for the Boundary element method of Stanley and Kato')
-
- parser.add_argument("--N",type=int, help="Surface size", required=True)
- parser.add_argument("--r",type=float, help="radius of hertz sphere", required=True)
- parser.add_argument("--steps",type=int, help="number of steps within which the pressure is applied.", required=True)
- parser.add_argument("--precision",type=float, help="relative precision, convergence if | (f_i - f_{i-1})/f_i | < Precision.", required=True)
- parser.add_argument("--e_star",type=float, help="effective elastic modulus" , required=True)
- parser.add_argument("--L",type=float, help="size of the surface" , required=True)
- parser.add_argument("--max_pressure",type=float, help="maximal load requested" , required=True)
- parser.add_argument("--plot_surface",type=bool, help="request output of text files containing the contact pressures on the surface" , default=False)
- parser.add_argument("--nthreads",type=int, help="request a number of threads to use via openmp to compute" , default=1)
-
- args = parser.parse_args()
- arguments = vars(args)
-
- n = arguments["N"]
- r = arguments["r"]
- max_pressure = arguments["max_pressure"]
- Ninc = arguments["steps"]
- epsilon = arguments["precision"]
- Estar = arguments["e_star"]
- L = arguments["L"]
- plot_surface = arguments["plot_surface"]
- nthreads = arguments["nthreads"]
-
- pressure = 0.
- dp = max_pressure/float(Ninc)
-
- s = np.zeros((n,n))
- print s.shape
- for i in range(0,n):
- for j in range(0,n):
- x = 1.*i - n/2
- y = 1.*j - n/2
- d = (x*x + y*y)*1./n/n*L*L
-
- if d < r*r: s[i,j] = - r + np.sqrt(r*r-d)
- else: s[i,j] = - r
-
-
- print "\n::: DATA ::::::::::::::::::::::::::::::\n"
- print " [N] {0}\n".format(n)
- print " [r] {0}\n".format(r)
- print " [Pext] {0}\n".format(max_pressure)
- print " [Steps] {0}\n".format(Ninc)
- print " [Precision] {0}\n".format(epsilon)
-
- bem = tamaas.BemPolonski(s)
- bem.setEffectiveModulus(Estar)
-
- file = open("hertz-prediction",'w')
-
- file.write("#pressure\tload\tcomputed-load\tarea_ratio\tarea\tradius\tradius-hertz\tradius/radius-hertz\tp0\tp0-hertz\n")
-
- for i in range(0,Ninc):
-
- pressure += dp
-
- bem.computeEquilibrium(epsilon,pressure)
- A = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
-
- dumpHertzPrediction(file,bem,pressure,r,Estar)
-
- if A == 1.0:
- file.close()
- break
-
- tamaas.dumpTimes()
-
-################################################################
-import sys
-import argparse
-import numpy as np
-main(sys.argv)
diff --git a/examples/legacy/surface_1d.py b/examples/legacy/surface_1d.py
deleted file mode 100644
index 6c48dca..0000000
--- a/examples/legacy/surface_1d.py
+++ /dev/null
@@ -1,83 +0,0 @@
-#!/usr/bin/env python
-# @file
-# @section LICENSE
-#
-# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
-# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-#
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see .
-
-from __future__ import division
-import tamaas as tm
-from numpy.fft import fft, fftshift, fftfreq
-import matplotlib.pyplot as plt
-
-tm.initialize()
-
-
-def draw_triangle(pos, w, ax, style, alpha=1.5, base=10):
- h = alpha * w
- ax.loglog([base**pos[0], base**(pos[0]+w)],
- [base**pos[1], base**pos[1]], style)
- ax.loglog([base**pos[0], base**pos[0]],
- [base**pos[1], base**(pos[1]+h)], style)
- ax.loglog([base**(pos[0]+w), base**pos[0]],
- [base**pos[1], base**(pos[1]+h)], style)
- # ax.text(base**(pos[0]+w/3), base**(pos[1]+h/3), "$-{}$".format(alpha),
- # horizontalalignment='center',
- # verticalalignment='center')
- ax.text(base**(pos[0]+w/2), base**(pos[1]-0.2), '$1$',
- horizontalalignment='center',
- verticalalignment='top')
- ax.text(base**(pos[0]-0.15), base**(pos[1]+h/2), '${}$'.format(alpha),
- horizontalalignment='right',
- verticalalignment='center')
-
-
-iso = tm.Isopowerlaw1D()
-helper = tm.ParamHelper(iso)
-
-params = {
- "Q0": 16,
- "Q1": 16,
- "Q2": 512,
- "Hurst": 0.8
-}
-
-helper.set(params)
-
-gen = tm.SurfaceGeneratorFilter1D()
-gen.setFilter(iso)
-gen.setSizes([2048])
-# gen.setRandomSeed(0) # uncomment for fixed seed
-
-surf = gen.buildSurface()
-
-plt.plot(surf)
-
-surf_fft = fft(surf)
-psd = surf_fft*surf_fft.conj()
-psd = fftshift(psd)
-
-plt.figure()
-freqs = fftshift(fftfreq(psd.size, d=1/psd.size))
-plt.plot(freqs[psd.size/2+1:], psd.real[psd.size/2+1:])
-draw_triangle([5, -4], 2, plt.gca(), '-k', 2*(params["Hurst"]+1), 2)
-plt.yscale('log')
-plt.xscale('log', basex=2)
-plt.ylim([10**(-2), 10**8])
-
-plt.show()
-
-tm.finalize()
diff --git a/examples/pipe_tools/surface b/examples/pipe_tools/surface
index a0d77f7..18bdfe5 100755
--- a/examples/pipe_tools/surface
+++ b/examples/pipe_tools/surface
@@ -1,81 +1,78 @@
#!/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-2020, EPFL (École Polytechnique Fédérale de Lausanne),"
"\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "AGPL"
__email__ = "lucas.frerot@gmail.com"
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()
-spectrum = tm.Isopowerlaw2D()
-spectrum.q0 = args.cutoffs[0]
-spectrum.q1 = args.cutoffs[1]
-spectrum.q2 = args.cutoffs[2]
-spectrum.hurst = args.hurst
-
if args.generator == 'random_phase':
- generator = tm.SurfaceGeneratorRandomPhase2D()
+ generator = tm.SurfaceGeneratorRandomPhase2D(args.sizes)
elif args.generator == 'filter':
- generator = tm.SurfaceGeneratorFilter2D()
+ generator = tm.SurfaceGeneratorFilter2D(args.sizes)
else:
raise ValueError('Unknown generator method {}'.format(args.generator))
-generator.setSizes(args.sizes)
+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
-generator.setSpectrum(spectrum)
-surface = generator.buildSurface() / spectrum.rmsSlopes() * args.rms
+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 31f1fd4..139973e 100644
--- a/examples/plasticity.py
+++ b/examples/plasticity.py
@@ -1,86 +1,85 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import numpy as np
import tamaas as tm
from tamaas.dumpers import H5Dumper as Dumper
from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver, ToleranceManager
-tm.initialize(1)
-
# 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 bdf967f..bc0ca84 100644
--- a/examples/rough_contact.py
+++ b/examples/rough_contact.py
@@ -1,71 +1,64 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
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()
-sg.setSizes([n, n])
+sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 0
# Spectrum
-spectrum = tm.Isopowerlaw2D()
-sg.setFilter(spectrum)
+sg.spectrum = tm.Isopowerlaw2D()
# Parameters
-spectrum.q0 = 16
-spectrum.q1 = 16
-spectrum.q2 = 64
-spectrum.hurst = 0.8
+sg.spectrum.q0 = 16
+sg.spectrum.q1 = 16
+sg.spectrum.q2 = 64
+sg.spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
-# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
-surface /= n
-# print(spectrum.rmsSlopes())
-# print(tm.SurfaceStatistics.computeRMSSlope(surface))
+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,
- tm.PolonskyKeerRey.pressure,
- tm.PolonskyKeerRey.pressure)
+solver = tm.PolonskyKeerRey(model, surface, 1e-12)
# Solve for target pressure
p_target = 0.1
solver.solve(p_target)
plt.figure()
-plt.imshow(model.getTraction())
+plt.imshow(model.traction)
plt.title('Contact tractions')
-print(model.getTraction().mean())
+print(model.traction.mean())
plt.show()
diff --git a/examples/saturated.py b/examples/saturated.py
index 81f4a59..dfad863 100644
--- a/examples/saturated.py
+++ b/examples/saturated.py
@@ -1,67 +1,66 @@
# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
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()
+sg = tm.SurfaceGeneratorFilter2D([grid_size, grid_size])
sg.random_seed = 2
-sg.setSizes([grid_size, grid_size])
-sg.setSpectrum(iso)
+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()
+elastic_tractions = model.traction.copy()
plt.imshow(elastic_tractions)
plt.title('Elastic contact tractions')
plt.colorbar()
-model['traction'][:] = 0.
+model.traction[:] = 0.
solver = tm.KatoSaturated(model, surface, 1e-12, p_sat)
-solver.setDumpFrequency(1)
-solver.setMaxIterations(200)
+solver.dump_freq = 1
+solver.max_iter = 200
solver.solve(load)
plt.figure()
-plt.imshow(model['traction'])
+plt.imshow(model.traction)
plt.title('Saturated contact tractions')
plt.colorbar()
plt.show()
diff --git a/examples/statistics.py b/examples/statistics.py
index 798a591..f47b220 100644
--- a/examples/statistics.py
+++ b/examples/statistics.py
@@ -1,86 +1,83 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import os
import numpy as np
import matplotlib.pyplot as plt
import tamaas as tm
from matplotlib.colors import LogNorm
-# Initialize threads and fftw
tm.set_log_level(tm.LogLevel.info) # Show progression of solver
# Surface size
n = 512
# Surface generator
-sg = tm.SurfaceGeneratorFilter2D()
-sg.setSizes([n, n])
+sg = tm.SurfaceGeneratorFilter2D([n, n])
sg.random_seed = 0
# Spectrum
-spectrum = tm.Isopowerlaw2D()
-sg.setFilter(spectrum)
+sg.spectrum = tm.Isopowerlaw2D()
# Parameters
-spectrum.q0 = 16
-spectrum.q1 = 16
-spectrum.q2 = 64
-spectrum.hurst = 0.8
+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 ca51412..9fbf92e 100644
--- a/examples/stresses.py
+++ b/examples/stresses.py
@@ -1,119 +1,122 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
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)
-local_size = model.getBoundaryDiscretization()
+# 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.getTraction()
+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.getDisplacement()
-stress = residual.getStress()
-gradient = residual.getVector()
-
-# Applying operator
-boussinesq = model.getIntegralOperator("boussinesq")
-boussinesq_gradient = model.getIntegralOperator("boussinesq_gradient")
-boussinesq.apply(traction, displacement)
-boussinesq_gradient.apply(traction, gradient)
-model.applyElasticity(stress, gradient)
+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 6859b58..088b8b9 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,162 +1,157 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
from __future__ import print_function
from os.path import abspath, join
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
""")
-# Adding legacy wrap code
-if env_pybind['legacy_bem']:
- env_pybind.AppendUnique(CPPDEFINES=['TAMAAS_LEGACY_BEM'])
- pybind_sources += ["wrap/bem.cpp"]
-
# Setting paths to find libTamaas
env_pybind.AppendUnique(LIBPATH=[Dir(join('#${build_dir}', '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
+ "'$$$$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 = """
compute.py
dumpers/__init__.py
dumpers/_helper.py
nonlinear_solvers/__init__.py
""".split()
targets = [tamaas_wrap]
targets += [
copy_env.Command(join('tamaas', f),
join(abspath(str(Dir('#python/tamaas'))), f),
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',
# 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)
conf = Configure(main_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
python_install = install_env.Command(
join('$prefix', 'dummy_target'),
targets, install_env['PYINSTALLCOM'], PYOPTIONS='--prefix ${prefix}',
chdir=install_env['PYDIR'])
python_install_dev = install_env.Command(
join('$prefix', 'dummy_target_local'),
targets, install_env['PYDEVELOPCOM'], PYOPTIONS='--user',
chdir=install_env['PYDIR'])
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 b97750c..b861add 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,211 +1,182 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __CAST_HH__
-#define __CAST_HH__
+#ifndef CAST_HH
+#define CAST_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
+#include "grid.hh"
#include "numpy.hh"
-#include "surface.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:
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:
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
}
};
-
-/**
- * Type caster for surface class
- */
-template
-struct type_caster> {
- using type = tamaas::Surface;
- using array_type =
- array_t;
-
-public:
- PYBIND11_TYPE_CASTER(type, _("SurfaceWrap"));
-
- 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::SurfaceNumpy(buf));
-
- return true;
- }
-
- static handle cast(const type& src, return_value_policy policy,
- handle parent) {
- return grid_to_python(src, policy, parent);
- }
-};
} // namespace detail
} // namespace pybind11
#endif // __CAST_HH__
diff --git a/python/tamaas/__init__.py.in b/python/tamaas/__init__.py.in
index 0605d7d..ee05ad3 100644
--- a/python/tamaas/__init__.py.in
+++ b/python/tamaas/__init__.py.in
@@ -1,83 +1,59 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
"""
Tamaas is a library dedicated to the fast treatment of rough contact problems.
-See README.md and examples for guidance as how to use Tamaas.
+See https://tamaas.readthedocs.io for documentation on how to use Tamaas.
-See __authors__, __license__, __copyright__ for extra information about Tamaas.
+See __author__, __license__, __copyright__ for extra information about Tamaas.
"""
__version__ = "@version@"
__author__ = @authors@
# TODO Change copyright when is issue with unicode is found
-__copyright__ = u"Copyright (©) 2016-20 EPFL " \
+__copyright__ = u"Copyright (©) 2016-2021 EPFL " \
+ u"(École Polytechnique Fédérale de Lausanne), " \
+ u"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
__license__ = "AGPLv3"
__maintainer__ = "Lucas Frérot"
__email__ = "@email@"
try:
from ._tamaas import model_type
from ._tamaas import _type_traits as __tt
- # Redefinition of model_type constants (for compatibility)
- # mark as deprecated
- 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
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
except ImportError as e:
print("Error trying to import _tamaas:\n{}".format(e))
raise e
-
-
-class ParamHelper:
- """Legacy class to manage parameters/setters/getters"""
- def __init__(self, obj):
- self.obj = obj
-
- def set(self, params):
- for key in params:
- setter_name = 'set' + key
- try:
- accessor = getattr(self.obj, setter_name)
- accessor(params[key])
- except Exception:
- print("Setter '{}({})' does not exist for object {}"
- .format(setter_name, type(params[key]), self.obj))
diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py
index 4b57833..49ee545 100644
--- a/python/tamaas/dumpers/__init__.py
+++ b/python/tamaas/dumpers/__init__.py
@@ -1,342 +1,342 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
"""
Dumpers for the class tamaas.Model
"""
from __future__ import print_function
from sys import stderr
from os import makedirs
import os.path
import numpy as np
from .. import ModelDumper, model_type, mpi, type_traits
from ._helper import step_dump, directory_dump, local_slice, _is_surface_field
def _get_attributes(model):
"Get model attributes"
return {
'model_type': str(model.type),
'system_size': model.system_size,
'discretization': model.global_shape,
}
class FieldDumper(ModelDumper):
"""Abstract dumper for python classes using fields"""
postfix = ''
extension = ''
name_format = "{basename}{postfix}.{extension}"
def __init__(self, basename, *fields, **kwargs):
"""Construct with desired fields"""
super(FieldDumper, self).__init__()
self.basename = basename
self.fields = list(fields)
self.all_fields = kwargs.get('all_fields', False)
def add_field(self, field):
"""Add another field to the dump"""
if field not in self.fields:
self.fields.append(field)
def dump_to_file(self, file_descriptor, model):
"""Dump to a file (name or handle)"""
def get_fields(self, model):
"""Get the desired fields"""
if not self.all_fields:
requested_fields = self.fields
else:
- requested_fields = model.getFields()
+ requested_fields = list(model)
- return {field: model.getField(field) for field in requested_fields}
+ return {field: model[field] for field in requested_fields}
def dump(self, model):
"Dump model"
self.dump_to_file(self.file_path, model)
@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, model):
"""Saving to compressed multi-field Numpy format"""
if mpi.size() > 1:
raise RuntimeError("NumpyDumper does not function "
"at all in parallel")
np.savez_compressed(file_descriptor, attrs=_get_attributes(model),
**self.get_fields(model))
try:
import h5py
@directory_dump('hdf5')
@step_dump
class H5Dumper(FieldDumper):
"""Dumper to HDF5 file format"""
extension = 'h5'
def dump_to_file(self, file_descriptor, model):
"""Saving to HDF5 with metadata about the model"""
# Setup for MPI
if not h5py.get_config().mpi and mpi.size() > 1:
raise RuntimeError("HDF5 does not have MPI support")
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)
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.COMM_WORLD.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
except ImportError:
pass
try:
import uvw # noqa
import uvw.parallel
@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, model):
"""Dump displacements, plastic deformations and stresses"""
if mpi.size() > 1:
raise RuntimeError("UVWDumper does not function "
"properly in parallel")
- bdim = len(model.getBoundaryDiscretization())
+ bdim = len(model.boundary_shape)
# Local MPI size
- lsize = model.getDiscretization()
- gsize = mpi.global_shape(model.getBoundaryDiscretization())
+ lsize = model.shape
+ gsize = mpi.global_shape(model.boundary_shape)
gshape = gsize
if len(lsize) > bdim:
- gshape = [model.getDiscretization()[0]] + gshape
+ gshape = [model.shape[0]] + gshape
# Space coordinates
coordinates = [np.linspace(0, L, N, endpoint=False)
- for L, N in zip(model.getSystemSize(), gshape)]
+ 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.getSystemSize()[0], gshape[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, *fields, **kwargs):
"""Construct with desired fields"""
super(UVWGroupDumper, self).__init__(basename, *fields, **kwargs)
subdir = os.path.join('paraview', basename + '-VTR')
if not os.path.exists(subdir):
makedirs(subdir)
self.uvw_dumper = UVWDumper(
os.path.join(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 as error:
print(error, file=stderr)
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):
grp.createDimension('frame', None)
# Local dimensions
- model_dim = len(model.getDiscretization())
+ model_dim = len(model.shape)
voigt_dim = type_traits[model.type].voigt
self._vec = grp.createDimension('spatial', model_dim)
self._tens = grp.createDimension('Voigt', voigt_dim)
- self.model_info = model.getDiscretization(), model.type
+ self.model_info = model.shape, model.type
# Create boundary dimensions
for label, size, length in zip(
"xy",
- model.getBoundaryDiscretization(),
- model.getBoundarySystemSize()
+ model.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),
- model.getBoundaryDiscretization(), "xy"
+ model.boundary_shape, "xy"
)
# Create volume dimension
if model.type in {model_type.volume_1d, model_type.volume_2d}:
- size = model.getDiscretization()[0]
+ size = model.shape[0]
grp.createDimension("z", size)
coord = grp.createVariable("z", 'f8', ("z",))
- coord[:] = np.linspace(0, model.getSystemSize()[0], size)
+ coord[:] = np.linspace(0, model.system_size[0], size)
self._create_variables(
grp, model,
lambda f: not _is_surface_field(f[1], model),
- model.getDiscretization(), "zxy"
+ model.shape, "zxy"
)
self.has_setup = True
def dump_to_file(self, file_descriptor, model):
if mpi.size() > 1:
raise RuntimeError("NetCDFDumper does not function "
"properly in parallel")
mode = 'a' if os.path.isfile(file_descriptor) \
and getattr(self, 'has_setup', False) else 'w'
with Dataset(file_descriptor, mode,
format='NETCDF4_CLASSIC') as rootgrp:
if rootgrp.dimensions == {}:
self._file_setup(rootgrp, model)
- if self.model_info != (model.getDiscretization(), model.type):
+ if self.model_info != (model.shape, model.type):
raise Exception("Unexpected model {}".format(model))
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]
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)
except ImportError:
pass
diff --git a/python/tamaas/dumpers/_helper.py b/python/tamaas/dumpers/_helper.py
index 4dcbcc8..60e3e0f 100644
--- a/python/tamaas/dumpers/_helper.py
+++ b/python/tamaas/dumpers/_helper.py
@@ -1,141 +1,141 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
"""
Helper functions for dumpers
"""
from functools import wraps
from collections import defaultdict
import os
import numpy as np
from .. import model_type, type_traits, mpi
__all__ = ["step_dump", "directory_dump"]
def _is_surface_field(field, model):
- bn = model.getBoundaryDiscretization()
+ bn = model.boundary_shape
return list(field.shape[:len(bn)]) == bn
def local_slice(field, model):
- n = model.getDiscretization()
- bn = model.getBoundaryDiscretization()
+ 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))
return tuple((slice(offset, offset+size, None)
for offset, size in zip(offsets, field.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"
def actual_decorator(cls):
orig_dump = cls.dump
orig_filepath = cls.file_path.fget
@wraps(cls.dump)
def dump(obj, *args, **kwargs):
if not os.path.exists(directory) and mpi.rank() == 0:
os.mkdir(directory)
if not os.path.isdir(directory) and mpi.rank() == 0:
raise Exception('{} is not a directory'.format(directory))
orig_dump(obj, *args, **kwargs)
@wraps(cls.file_path.fget)
def file_path(obj):
return os.path.join(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_module.cpp b/python/tamaas_module.cpp
index b7b4885..eebb496 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,92 +1,80 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#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() = "Tamaas module for python";
// Wrapping the base methods
mod.def("initialize", &initialize, py::arg("num_threads") = 0,
"Initialize tamaas with desired number of threads");
mod.def("finalize", &finalize, "Final cleanup");
// 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;
// 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);
-
-#if defined(TAMAAS_LEGACY_BEM)
- // Legacy wrap
- py::class_>(mod, "smart_pointer_UInt")
- .def("assign", &wrap::smart_pointer::assign);
- py::class_>(mod, "smart_pointer_Real")
- .def("assign", &wrap::smart_pointer::assign);
- py::class_>(mod, "smart_pointer_long")
- .def("assign", &wrap::smart_pointer::assign);
-
- wrap::wrapBEM(mod);
-#endif
}
} // namespace tamaas
diff --git a/python/wrap.hh b/python/wrap.hh
index 5756359..6c54053 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,75 +1,77 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __WRAP_HH__
#define __WRAP_HH__
/* -------------------------------------------------------------------------- */
-#include "tamaas.hh"
-#include "numpy.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();
}
-template class smart_pointer {
-public:
- smart_pointer(T* ptr) : ptr(ptr) {}
- void assign(T val) { *ptr = val; }
-
-private:
- T* ptr;
-};
-
/* -------------------------------------------------------------------------- */
/* 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);
-#if defined(TAMAAS_LEGACY_BEM)
-void wrapBEM(py::module& mod);
-#endif
+/* -------------------------------------------------------------------------- */
+/* 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/bem.cpp b/python/wrap/bem.cpp
deleted file mode 100644
index 4f0d56c..0000000
--- a/python/wrap/bem.cpp
+++ /dev/null
@@ -1,115 +0,0 @@
-/**
- * @file
- * @section LICENSE
- *
- * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
- * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Affero General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU Affero General Public License for more details.
- *
- * You should have received a copy of the GNU Affero General Public License
- * along with this program. If not, see .
- *
- */
-/* -------------------------------------------------------------------------- */
-#include "bem_fft_base.hh"
-#include "bem_gigipol.hh"
-#include "bem_grid.hh"
-#include "bem_grid_condat.hh"
-#include "bem_grid_polonski.hh"
-#include "bem_kato.hh"
-#include "bem_polonski.hh"
-#include "wrap.hh"
-/* -------------------------------------------------------------------------- */
-
-namespace tamaas {
-
-namespace wrap {
-
-using namespace py::literals;
-
-/* -------------------------------------------------------------------------- */
-void wrapBEM(py::module& mod) {
- py::class_(mod, "BemFFTBase")
- .def("computeDisplacementsFromTractions",
- &BemFFTBase::computeDisplacementsFromTractions)
- .def("computeTractionsFromDisplacements",
- &BemFFTBase::computeTractionsFromDisplacements)
- .def("applyInfluenceFunctions", &BemFFTBase::applyInfluenceFunctions)
- .def("applyInverseInfluenceFunctions",
- &BemFFTBase::applyInverseInfluenceFunctions)
- .def("computeSpectralInfluenceOverDisplacements",
- &BemFFTBase::computeSpectralInfluenceOverDisplacement)
- .def("computeTrueDisplacements", &BemFFTBase::computeTrueDisplacements)
- .def("computeGaps", &BemFFTBase::computeGaps)
- .def("getSpectralInfluenceOverDisplacement",
- &BemFFTBase::getSpectralInfluenceOverDisplacement)
- .def("getTractions", &BemFFTBase::getTractions,
- py::return_value_policy::reference_internal)
- .def("getDisplacements", &BemFFTBase::getDisplacements,
- py::return_value_policy::reference_internal)
- .def("getTrueDisplacements", &BemFFTBase::getTrueDisplacements,
- py::return_value_policy::reference_internal)
- .def("getGap", &BemFFTBase::getGap,
- py::return_value_policy::reference_internal)
- .def("setEffectiveModulus", &BemFFTBase::setEffectiveModulus, "E_star"_a)
- .def("getEffectiveModulus", &BemFFTBase::getEffectiveModulus)
- .def("setMaxIterations", &BemFFTBase::setMaxIterations)
- .def("getMaxIterations", &BemFFTBase::getMaxIterations)
- .def("setDumpFreq", &BemFFTBase::setDumpFreq)
- .def("getDumpFreq", &BemFFTBase::getDumpFreq);
-
- py::class_(mod, "BemPolonski")
- .def(py::init&>())
- .def("computeEquilibrium", &BemPolonski::computeEquilibrium)
- .def("computeEquilibriuminit", &BemPolonski::computeEquilibriuminit)
- .def("computeMeanGapsInContact", &BemPolonski::computeMeanGapsInContact)
- .def("setMaxPressure", &BemPolonski::setMaxPressure)
- .def("getNbIterations", &BemPolonski::getNbIterations);
-
- py::class_(mod, "BemKato")
- .def(py::init&>())
- .def("computeEquilibrium", &BemKato::computeEquilibrium)
- .def("setMaxPressure", &BemKato::setMaxPressure);
-
- py::class_(mod, "BemGigiPol")
- .def(py::init&>())
- .def("computeEquilibrium", &BemGigipol::computeEquilibrium)
- .def("computeEquilibrium2", &BemGigipol::computeEquilibrium2)
- .def("computeEquilibriuminit", &BemGigipol::computeEquilibriuminit)
- .def("computeEquilibrium2init", &BemGigipol::computeEquilibrium2init)
- .def("computeTractionsFromDisplacements",
- &BemGigipol::computeTractionsFromDisplacements);
-
- py::class_(mod, "BemGrid")
- .def("computeDisplacementsFromTractions",
- &BemGrid::computeDisplacementsFromTractions)
- .def("computeSurfaceNormals", &BemGrid::computeSurfaceNormals)
- .def("computeInfluence", &BemGrid::computeInfluence)
- .def("setElasticity", &BemGrid::setElasticity)
- .def("getTractions", &BemGrid::getTractions,
- py::return_value_policy::reference_internal)
- .def("getDisplacements", &BemGrid::getDisplacements,
- py::return_value_policy::reference_internal)
- .def("getSurfaceNormals", &BemGrid::getSurfaceNormals,
- py::return_value_policy::reference_internal);
-
- py::class_(mod, "BemGridPolonski")
- .def(py::init&>());
-
- py::class_(mod, "BemGridCondat")
- .def(py::init&>())
- .def("computeEquilibrium", &BemGridCondat::computeEquilibrium);
-}
-
-} // namespace wrap
-
-} // namespace tamaas
diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index feec9ad..775e840 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,138 +1,105 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "logger.hh"
#include "statistics.hh"
-#include "surface_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)
.def_static("computeAutocorrelation",
&Statistics::computeAutocorrelation,
py::return_value_policy::move)
.def_static("computeMoments", &Statistics::computeMoments)
.def_static("computeSpectralRMSSlope",
&Statistics::computeSpectralRMSSlope)
.def_static("computeRMSHeights", &Statistics::computeRMSHeights)
.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;
})
.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;
});
- // Exposing SurfaceStatistics (legacy)
-#if defined(TAMAAS_LEGACY_BEM)
- py::class_(mod, "SurfaceStatistics")
- .def_static("computeMaximum", &SurfaceStatistics::computeMaximum)
- .def_static("computeMinimum", &SurfaceStatistics::computeMinimum)
- .def_static("computeSpectralRMSSlope",
- &SurfaceStatistics::computeSpectralRMSSlope)
- .def_static("computeRMSSlope", &SurfaceStatistics::computeRMSSlope)
- .def_static("computeMoments", &SurfaceStatistics::computeMoments)
- .def_static("computeSkewness", &SurfaceStatistics::computeSkewness)
- .def_static("computeKurtosis", &SurfaceStatistics::computeKurtosis)
- .def_static("computeSpectralMeanCurvature",
- &SurfaceStatistics::computeSpectralMeanCurvature)
- .def_static("computeSpectralStdev",
- &SurfaceStatistics::computeSpectralStdev)
- .def_static("computeAnalyticFractalMoment",
- &SurfaceStatistics::computeAnalyticFractalMoment)
- .def_static("computePerimeter", &SurfaceStatistics::computePerimeter)
- .def_static("computeContactArea", &SurfaceStatistics::computeContactArea)
- .def_static("computeContactAreaRatio",
- &SurfaceStatistics::computeContactAreaRatio)
- .def_static("computeSpectralDistribution",
- &SurfaceStatistics::computeSpectralDistribution)
- .def_static("computeSum", &SurfaceStatistics::computeSum)
- .def_static("computeAutocorrelation",
- &SurfaceStatistics::computeAutocorrelation,
- py::return_value_policy::copy)
- .def_static("computePowerSpectrum",
- &SurfaceStatistics::computePowerSpectrum,
- py::return_value_policy::copy);
-#endif
-
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
mod.def("to_voigt",
[](const Grid& field) {
Logger().get(LogLevel::warning)
<< "tamaas.to_voigt deprecated. Use tamaas.compute.to_voigt\n";
if (field.getNbComponents() == 9) {
Grid voigt(field.sizes(), 6);
Loop::loop([](auto in, auto out) { out.symmetrize(in); },
range>(field),
range>(voigt));
return voigt;
} else
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 da53cc2..0ae985c 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,320 +1,431 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#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(std::string name) {
+ return m.getIntegralOperator(std::move(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)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
- // legacy wrapper code
- .def("setParameters", &functional::AdhesionFunctional::setParameters);
+ .def("setParameters",
+ [](functional::AdhesionFunctional& f,
+ const std::map& m) {
+ TAMAAS_DEPRECATE("setParameters()", "the parameters property");
+ f.setParameters(m);
+ });
py::class_(
mod, "ExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh)
.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("updateFromModel", &IntegralOperator::updateFromModel)
- .def("getModel", &IntegralOperator::getModel,
+ .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"),
py::return_value_policy::reference)
- .def("getKind", &IntegralOperator::getKind)
- .def("getType", &IntegralOperator::getType);
+ .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);
+ })
+ .def("updateFromModel", &IntegralOperator::updateFromModel)
+
+ .def_property_readonly("kind", &IntegralOperator::getKind)
+ .def_property_readonly("model", &IntegralOperator::getModel)
+ .def_property_readonly("type", &IntegralOperator::getType);
py::enum_(mod, "integration_method")
.value("linear", integration_method::linear)
.value("cutoff", integration_method::cutoff);
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
- .def("getModel", &BEEngine::getModel, py::return_value_policy::reference)
.def("registerNeumann", &BEEngine::registerNeumann)
- .def("registerDirichlet", &BEEngine::registerDirichlet);
+ .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; })
.def_property_readonly_static(
"components", [](py::object) { return trait::components; })
.def_property_readonly_static(
"boundary_dimension",
[](py::object) { return trait::boundary_dimension; })
.def_property_readonly_static("voigt",
[](py::object) { return trait::voigt; })
.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)
.value("basic_2d", model_type::basic_2d)
.value("surface_1d", model_type::surface_1d)
.value("surface_2d", model_type::surface_2d)
.value("volume_1d", model_type::volume_1d)
.value("volume_2d", model_type::volume_2d);
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("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
- .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus)
- .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio)
- .def("getHertzModulus", &Model::getHertzModulus)
- .def("getYoungModulus", &Model::getYoungModulus)
- .def("getShearModulus", &Model::getShearModulus)
- .def("getPoissonRatio", &Model::getPoissonRatio)
- .def("getTraction", (GridBase & (Model::*)()) & Model::getTraction,
+ .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)
+
+ .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("getDisplacement",
- (GridBase & (Model::*)()) & Model::getDisplacement,
+ .def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
py::return_value_policy::reference_internal)
- .def("getSystemSize", &Model::getSystemSize)
- .def("getDiscretization", &Model::getDiscretization)
- .def("getBoundarySystemSize", &Model::getBoundarySystemSize)
- .def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
+ .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)
.def("solveDirichlet", &Model::solveDirichlet)
.def("dump", &Model::dump)
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
- .def("setDumper",
- [](Model&, const ModelDumper&) {
- TAMAAS_EXCEPTION("setDumper() is not a member of Model; use "
- "addDumper() instead");
- })
- .def("getBEEngine", &Model::getBEEngine,
+
+ .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", &Model::getIntegralOperator,
+ .def("getIntegralOperator",
+ [](const Model& m, std::string name) {
+ TAMAAS_DEPRECATE("getIntegralOperator()",
+ "the operators property");
+ return m.getIntegralOperator(std::move(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",
- (GridBase & (Model::*)(const std::string&)) & Model::getField,
+ [](const Model& m, std::string name) -> decltype(m.getField(name)) {
+ TAMAAS_DEPRECATE("getField()", "the [] operator");
+ return m.getField(std::move(name));
+ },
"field_name"_a, py::return_value_policy::reference_internal)
- .def("getFields", &Model::getFields)
+ .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
+ // Python magic functions
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def("__getitem__",
- [](const Model& m, std::string key) -> const GridBase& {
+ [](const Model& m, std::string key) -> decltype(m[key]) {
try {
- return m.getField(key);
+ return m[key];
} catch (std::out_of_range&) {
- throw py::key_error();
+ throw py::key_error(key);
}
},
- py::return_value_policy::reference_internal)
+ 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) !=
+ return std::find(fields.begin(), fields.end(), std::move(key)) !=
fields.end();
},
- py::keep_alive<0, 1>())
+ 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>())
-
- // More python-like access to model properties
- .def_property_readonly("shape", &Model::getDiscretization)
- .def_property_readonly("global_shape", &Model::getGlobalDiscretization)
+ 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)
- .def_property_readonly("system_size", &Model::getSystemSize);
+ &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)
+ .def("dump", &ModelDumper::dump, "model"_a, "Dump model")
.def("__lshift__",
- [](ModelDumper& dumper, Model& model) { dumper << model; });
+ [](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)
- .def_static("createResidual", &ModelFactory::createResidual,
- "model"_a, "sigma_y"_a, "hardening"_a = 0.)
+ "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);
+ &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/percolation.cpp b/python/wrap/percolation.cpp
index a3d770b..25924a4 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,86 +1,93 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#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("getArea", &Cluster::getArea, "Area of cluster")
- .def("getPoints", &Cluster::getPoints,
- "Get list of points of cluster")
- .def("getPerimeter", &Cluster::getPerimeter,
- "Get perimeter of cluster")
+ .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 7a85807..c7a9187 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,171 +1,188 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#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;
void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
void updateState() override {
PYBIND11_OVERLOAD(void, EPSolver, updateState);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_(mod, "ContactSolver")
- // .def(py::init&, Real>(), "model"_a,
- // "surface"_a, "tolerance"_a)
- .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a)
- .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a)
+ .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)
.def("solve", py::overload_cast>(&ContactSolver::solve),
"target_force"_a,
py::call_guard())
.def("solve", py::overload_cast(&ContactSolver::solve),
"target_normal_pressure"_a,
py::call_guard());
py::class_ pkr(mod, "PolonskyKeerRey");
// 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")
.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);
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");
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")
.def(py::init());
py::class_(mod, "EPSolver")
.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")
.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())
.def("acceleratedSolve",
[](EPICSolver& solver, Real pressure) {
return solver.acceleratedSolve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard());
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 4bfa259..6898229 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,202 +1,170 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include "regularized_powerlaw.hh"
#include "surface_generator_filter.hh"
-#include "surface_generator_filter_fft.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 {
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())
.def(py::init<>())
.def("computeFilter",
(void (Filter::*)(GridHermitian&) const) &
Filter::computeFilter);
}
/* -------------------------------------------------------------------------- */
template
void wrapIsopowerlaw(py::module& mod) {
std::string name = makeDimensionName("Isopowerlaw", dim);
py::class_, Filter, std::shared_ptr>>(
mod, name.c_str())
.def(py::init<>())
- .def_property("q0", &Isopowerlaw::getQ0, &Isopowerlaw::setQ0)
- .def_property("q1", &Isopowerlaw::getQ1, &Isopowerlaw::setQ1)
- .def_property("q2", &Isopowerlaw::getQ2, &Isopowerlaw::setQ2)
+ .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)
- .def("rmsHeights", &Isopowerlaw::rmsHeights)
- .def("moments", &Isopowerlaw::moments)
- .def("alpha", &Isopowerlaw::alpha)
- .def("rmsSlopes", &Isopowerlaw