diff --git a/SConstruct b/SConstruct
index 0cfb893..3c025be 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,477 +1,477 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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-2020 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)"
)
# ------------------------------------------------------------------------------
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/examples/adhesion.py b/examples/adhesion.py
index b7ca58f..dfc8ea4 100644
--- a/examples/adhesion.py
+++ b/examples/adhesion.py
@@ -1,127 +1,127 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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)
# Spectrum
spectrum = tm.Isopowerlaw2D()
sg.setFilter(spectrum)
# Parameters
spectrum.q0 = 16
spectrum.q1 = 16
spectrum.q2 = 64
spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
surface /= 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()
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
index 8bc9036..64ed1d5 100644
--- a/examples/legacy/adhesion.py
+++ b/examples/legacy/adhesion.py
@@ -1,161 +1,161 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
index a03370b..8b1bca2 100644
--- a/examples/legacy/cluster_statistics.py
+++ b/examples/legacy/cluster_statistics.py
@@ -1,104 +1,104 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
index 25abd57..0adbc84 100644
--- a/examples/legacy/contact.py
+++ b/examples/legacy/contact.py
@@ -1,73 +1,73 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
index b68fb4c..0a07a9c 100644
--- a/examples/legacy/fractal_surface.py
+++ b/examples/legacy/fractal_surface.py
@@ -1,102 +1,102 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
index cb5bfaa..2a30094 100644
--- a/examples/legacy/hertz.py
+++ b/examples/legacy/hertz.py
@@ -1,116 +1,116 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
index 8cdd108..6c48dca 100644
--- a/examples/legacy/surface_1d.py
+++ b/examples/legacy/surface_1d.py
@@ -1,83 +1,83 @@
#!/usr/bin/env python
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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/plasticity.py b/examples/plasticity.py
index 3896bee..31f1fd4 100644
--- a/examples/plasticity.py
+++ b/examples/plasticity.py
@@ -1,86 +1,86 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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)
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 4cef430..bdf967f 100644
--- a/examples/rough_contact.py
+++ b/examples/rough_contact.py
@@ -1,71 +1,71 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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.random_seed = 0
# Spectrum
spectrum = tm.Isopowerlaw2D()
sg.setFilter(spectrum)
# Parameters
spectrum.q0 = 16
spectrum.q1 = 16
spectrum.q2 = 64
spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
surface /= n
# print(spectrum.rmsSlopes())
# print(tm.SurfaceStatistics.computeRMSSlope(surface))
plt.imshow(surface)
# Creating model
model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n])
# Solver
solver = tm.PolonskyKeerRey(model, surface, 1e-12,
tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.pressure)
# Solve for target pressure
p_target = 0.1
solver.solve(p_target)
plt.figure()
plt.imshow(model.getTraction())
plt.title('Contact tractions')
print(model.getTraction().mean())
plt.show()
diff --git a/examples/saturated.py b/examples/saturated.py
index e576bde..81f4a59 100644
--- a/examples/saturated.py
+++ b/examples/saturated.py
@@ -1,67 +1,67 @@
# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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.random_seed = 2
sg.setSizes([grid_size, grid_size])
sg.setSpectrum(iso)
surface = sg.buildSurface()
surface -= np.max(surface)
model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
[1., 1.], [grid_size, grid_size])
model.E = 1.
model.nu = 0
esolver = tm.PolonskyKeerRey(model, surface, 1e-12,
tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.pressure)
esolver.solve(load)
elastic_tractions = model['traction'].copy()
plt.imshow(elastic_tractions)
plt.title('Elastic contact tractions')
plt.colorbar()
model['traction'][:] = 0.
solver = tm.KatoSaturated(model, surface, 1e-12, p_sat)
solver.setDumpFrequency(1)
solver.setMaxIterations(200)
solver.solve(load)
plt.figure()
plt.imshow(model['traction'])
plt.title('Saturated contact tractions')
plt.colorbar()
plt.show()
diff --git a/examples/statistics.py b/examples/statistics.py
index e1e37c1..798a591 100644
--- a/examples/statistics.py
+++ b/examples/statistics.py
@@ -1,86 +1,86 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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.random_seed = 0
# Spectrum
spectrum = tm.Isopowerlaw2D()
sg.setFilter(spectrum)
# Parameters
spectrum.q0 = 16
spectrum.q1 = 16
spectrum.q2 = 64
spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
# Computing PSD and shifting for plot
psd = tm.Statistics2D.computePowerSpectrum(surface)
psd = np.fft.fftshift(psd, axes=0)
plt.imshow(psd.real, norm=LogNorm())
plt.gca().set_title('Power Spectrum Density')
plt.gcf().tight_layout()
# Computing autocorrelation and shifting for plot
acf = tm.Statistics2D.computeAutocorrelation(surface)
acf = np.fft.fftshift(acf)
plt.figure()
plt.imshow(acf)
plt.gca().set_title('Autocorrelation')
plt.gcf().tight_layout()
plt.show()
# Write the rough surface in paraview's VTK format
try:
import uvw
try:
os.mkdir('paraview')
except FileExistsError:
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 53151ff..ca51412 100644
--- a/examples/stresses.py
+++ b/examples/stresses.py
@@ -1,119 +1,119 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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()
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()
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)
# 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 123a023..6859b58 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,162 +1,162 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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
])
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 c672217..b97750c 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,211 +1,211 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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__
/* -------------------------------------------------------------------------- */
#include "grid_base.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/numpy.hh b/python/numpy.hh
index 1f196ea..80a840f 100644
--- a/python/numpy.hh
+++ b/python/numpy.hh
@@ -1,93 +1,93 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 __NUMPY_HH__
#define __NUMPY_HH__
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
#include
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// Type alias for usual numpy array type
template
using numpy = py::array_t;
/// GridBase helper class wrapping numpy array
template
class GridBaseNumpy : public GridBase {
public:
GridBaseNumpy(numpy& buffer) : GridBase() {
this->nb_components = 1;
this->data.wrapMemory(buffer.mutable_data(), buffer.size());
}
};
/// Grid helper class wrapping numpy array
template
class GridNumpy : public Parent {
public:
GridNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
const UInt ndim = buffer.ndim();
if (ndim > Parent::dimension + 1 || ndim < Parent::dimension)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected grid dimensions");
if (ndim == Parent::dimension + 1)
this->nb_components = array_shape[ndim - 1];
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
/// Surface helper class wrapping numpy array
template
class SurfaceNumpy : public Parent {
public:
SurfaceNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
if (buffer.ndim() != 2)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected surface dimensions");
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
} // namespace wrap
} // namespace tamaas
#endif // __NUMPY_HH__
diff --git a/python/tamaas/__init__.py.in b/python/tamaas/__init__.py.in
index 646df6a..0605d7d 100644
--- a/python/tamaas/__init__.py.in
+++ b/python/tamaas/__init__.py.in
@@ -1,83 +1,83 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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 __authors__, __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 " \
+ 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/compute.py b/python/tamaas/compute.py
index 765987c..b300a90 100644
--- a/python/tamaas/compute.py
+++ b/python/tamaas/compute.py
@@ -1,22 +1,22 @@
# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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.compute import * # noqa
diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py
index 2d29c6e..4b57833 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-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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()
return {field: model.getField(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())
# Local MPI size
lsize = model.getDiscretization()
gsize = mpi.global_shape(model.getBoundaryDiscretization())
gshape = gsize
if len(lsize) > bdim:
gshape = [model.getDiscretization()[0]] + gshape
# Space coordinates
coordinates = [np.linspace(0, L, N, endpoint=False)
for L, N in zip(model.getSystemSize(), 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])
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())
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
# Create boundary dimensions
for label, size, length in zip(
"xy",
model.getBoundaryDiscretization(),
model.getBoundarySystemSize()
):
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"
)
# Create volume dimension
if model.type in {model_type.volume_1d, model_type.volume_2d}:
size = model.getDiscretization()[0]
grp.createDimension("z", size)
coord = grp.createVariable("z", 'f8', ("z",))
coord[:] = np.linspace(0, model.getSystemSize()[0], size)
self._create_variables(
grp, model,
lambda f: not _is_surface_field(f[1], model),
model.getDiscretization(), "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):
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 66bbbcb..4dcbcc8 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-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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()
return list(field.shape[:len(bn)]) == bn
def local_slice(field, model):
n = model.getDiscretization()
bn = model.getBoundaryDiscretization()
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/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py
index c911dc3..e4e6368 100644
--- a/python/tamaas/nonlinear_solvers/__init__.py
+++ b/python/tamaas/nonlinear_solvers/__init__.py
@@ -1,168 +1,168 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
-# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+# 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 .
"""
Pulling solvers to nonlinear_solvers module
"""
from functools import wraps
import numpy as np
from scipy.sparse.linalg import LinearOperator, lgmres
from scipy.linalg import norm
from scipy.optimize import newton_krylov, root
from scipy.optimize.nonlin import NoConvergence
from .. import EPSolver, Logger, LogLevel, mpi
from .._tamaas import _tolerance_manager
from .._tamaas import _DFSANESolver as DFSANECXXSolver
__all__ = ['NLNoConvergence',
'DFSANESolver',
'DFSANECXXSolver',
'NewtonKrylovSolver',
'ToleranceManager']
class NLNoConvergence(Exception):
"""Convergence not reached exception"""
class ScipySolver(EPSolver):
"""
Base class for solvers wrapping SciPy routines
"""
def __init__(self, residual, _, callback=None):
super(ScipySolver, self).__init__(residual)
if mpi.size() > 1:
raise RuntimeError("Scipy solvers cannot be used with MPI; "
"DFSANECXXSolver can be used instead")
self.callback = callback
self._x = self.getStrainIncrement()
self._residual = self.getResidual()
self.options = {'ftol': 0, 'fatol': 1e-9}
def solve(self):
"""
Solve the nonlinear plasticity equation using the scipy_solve routine
"""
# For initial guess, compute the strain due to boundary tractions
# self._residual.computeResidual(self._x)
# self._x[...] = self._residual.getVector()
EPSolver.beforeSolve(self)
# Scipy root callback
def compute_residual(vec):
self._residual.computeResidual(vec)
return self._residual.getVector().copy()
# Solve
self._x[...] = self.scipy_solve(compute_residual)
# Computing displacements
self._residual.computeResidualDisplacement(self._x)
def reset(self):
"Set solution vector to zero"
self._x[...] = 0
class NewtonKrylovSolver(ScipySolver):
"""
Solve using a finite-difference Newton-Krylov method
"""
def __init__(self, residual, model=None, callback=None):
ScipySolver.__init__(self, residual, model,
callback=callback)
def scipy_solve(self, compute_residual):
"Solve R(delta epsilon) = 0 using a newton-krylov method"
try:
return newton_krylov(compute_residual, self._x,
f_tol=self.tolerance,
verbose=True, callback=self.callback)
except NoConvergence:
raise NLNoConvergence("Newton-Krylov did not converge")
class DFSANESolver(ScipySolver):
"""
Solve using a spectral residual jacobianless method
"""
def __init__(self, residual, model=None, callback=None):
ScipySolver.__init__(self, residual, model,
callback=callback)
def scipy_solve(self, compute_residual):
"Solve R(delta epsilon) = 0 using a df-sane method"
solution = root(compute_residual,
self._x,
method='df-sane',
options={'ftol': 0, 'fatol': self.tolerance},
callback=self.callback)
Logger().get(LogLevel.info) << \
"DF-SANE/Scipy: {} ({} iterations, {})".format(
solution.message,
solution.nit,
self.tolerance)
if not solution.success:
raise NLNoConvergence("DF-SANE/Scipy did not converge")
return solution.x.copy()
def ToleranceManager(start, end, rate):
"Decorator to manage tolerance of non-linear solver"
# start /= rate # just anticipating first multiplication
def actual_decorator(cls):
orig_init = cls.__init__
orig_solve = cls.solve
orig_update_state = cls.updateState
@wraps(cls.__init__)
def __init__(obj, *args, **kwargs):
orig_init(obj, *args, **kwargs)
obj.setToleranceManager(_tolerance_manager(start, end, rate))
@wraps(cls.solve)
def new_solve(obj, *args, **kwargs):
ftol = obj.tolerance
ftol *= rate
obj.tolerance = max(ftol, end)
return orig_solve(obj, *args, **kwargs)
@wraps(cls.updateState)
def updateState(obj, *args, **kwargs):
obj.tolerance = start
return orig_update_state(obj, *args, **kwargs)
cls.__init__ = __init__
# cls.solve = new_solve
# cls.updateState = updateState
return cls
return actual_decorator
diff --git a/python/tamaas_module.cpp b/python/tamaas_module.cpp
index 60d2fc6..b7b4885 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,92 +1,92 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 820fa33..5756359 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,75 +1,75 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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
/* -------------------------------------------------------------------------- */
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
} // namespace wrap
} // namespace tamaas
#endif // __WRAP_HH__
diff --git a/python/wrap/bem.cpp b/python/wrap/bem.cpp
index 1b83d3f..4f0d56c 100644
--- a/python/wrap/bem.cpp
+++ b/python/wrap/bem.cpp
@@ -1,115 +1,115 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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/compute.cpp b/python/wrap/compute.cpp
index b7444cb..126fd69 100644
--- a/python/wrap/compute.cpp
+++ b/python/wrap/compute.cpp
@@ -1,67 +1,67 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 "computes.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
void wrapCompute(py::module& mod) {
auto compute_mod = mod.def_submodule("compute");
compute_mod.def(
"eigenvalues",
[](model_type type, Grid& eigs, const Grid& field) {
eigenvalues(type, eigs, field);
},
"model_type"_a, "eigenvalues_out"_a, "field"_a);
compute_mod.def("von_mises",
[](model_type type, Grid& vm,
const Grid& field) { vonMises(type, vm, field); },
"model_type"_a, "von_mises"_a, "field"_a);
compute_mod.def(
"deviatoric",
[](model_type type, Grid& dev, const Grid& field) {
deviatoric(type, dev, field);
},
"model_type"_a, "deviatoric"_a, "field"_a);
compute_mod.def("to_voigt",
[](const Grid& field) {
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/core.cpp b/python/wrap/core.cpp
index d0f8562..feec9ad 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,138 +1,138 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 2afa2f6..da53cc2 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,320 +1,320 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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;
/// 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);
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) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def("updateFromModel", &IntegralOperator::updateFromModel)
.def("getModel", &IntegralOperator::getModel,
py::return_value_policy::reference)
.def("getKind", &IntegralOperator::getKind)
.def("getType", &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);
}
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")
.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,
py::return_value_policy::reference_internal)
.def("getDisplacement",
(GridBase & (Model::*)()) & Model::getDisplacement,
py::return_value_policy::reference_internal)
.def("getSystemSize", &Model::getSystemSize)
.def("getDiscretization", &Model::getDiscretization)
.def("getBoundarySystemSize", &Model::getBoundarySystemSize)
.def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
.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,
py::return_value_policy::reference_internal)
.def("getIntegralOperator", &Model::getIntegralOperator,
"operator_name"_a, py::return_value_policy::reference_internal)
.def("registerField",
[](Model& m, std::string name, numpy field) {
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,
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields", &Model::getFields)
.def("applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
})
// 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& {
try {
return m.getField(key);
} catch (std::out_of_range&) {
throw py::key_error();
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), key) !=
fields.end();
},
py::keep_alive<0, 1>())
.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)
.def_property_readonly("boundary_shape",
&Model::getBoundaryDiscretization)
.def_property_readonly("system_size", &Model::getSystemSize);
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a)
.def("__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << 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.)
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a);
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
"cutoff"_a = 1e-12)
.def_property("yield_stress", &Residual::getYieldStress,
&Residual::setYieldStress)
.def_property("hardening_modulus", &Residual::getHardeningModulus,
&Residual::setHardeningModulus)
.def_property_readonly("model", &Residual::getModel);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/model_extensions.hh b/python/wrap/model_extensions.hh
index 6b9d39b..7da97c6 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,132 +1,132 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 "functional.hh"
#include "model.hh"
#include "model_dumper.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace functional {
namespace wrap {
/// Class for extension of Functional in python
class PyFunctional : public Functional {
public:
using Functional::Functional;
// Overriding pure virtual functions
Real computeF(GridBase& variable, GridBase& dual) const override {
PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
}
void computeGradF(GridBase& variable,
GridBase& gradient) const override {
PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
}
};
} // namespace wrap
} // namespace functional
/* -------------------------------------------------------------------------- */
namespace wrap {
/* -------------------------------------------------------------------------- */
class PyModelDumper : public ModelDumper {
public:
using ModelDumper::ModelDumper;
void dump(const Model& model) override {
PYBIND11_OVERLOAD_PURE(void, ModelDumper, dump, model);
}
};
/* -------------------------------------------------------------------------- */
class PyResidual : public Residual {
public:
using Residual::Residual;
void computeResidual(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidual, strain_increment);
}
void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, applyTangent, output, input,
current_strain_increment);
}
void computeStress(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeStress, strain_increment);
}
void updateState(GridBase& converged_strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, updateState,
converged_strain_increment);
}
const GridBase& getVector() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getVector);
}
const GridBase& getPlasticStrain() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getPlasticStrain);
}
const GridBase& getStress() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getStress);
}
void computeResidualDisplacement(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidualDisplacement,
strain_increment);
}
void setIntegrationMethod(integration_method method, Real cutoff) override {
PYBIND11_OVERLOAD_PURE(void, Residual, setIntegrationMethod, method,
cutoff);
}
Real getYieldStress() const override {
PYBIND11_OVERLOAD_PURE(Real, Residual, getYieldStress);
}
Real getHardeningModulus() const override {
PYBIND11_OVERLOAD_PURE(Real, Residual, getHardeningModulus);
}
void setYieldStress(Real val) override {
PYBIND11_OVERLOAD_PURE(void, Residual, setYieldStress, val);
}
void setHardeningModulus(Real val) override {
PYBIND11_OVERLOAD_PURE(void, Residual, setHardeningModulus, val);
}
};
/* -------------------------------------------------------------------------- */
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/mpi.cpp b/python/wrap/mpi.cpp
index 2b03c36..96c1262 100644
--- a/python/wrap/mpi.cpp
+++ b/python/wrap/mpi.cpp
@@ -1,89 +1,89 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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 "mpi_interface.hh"
#include "partitioner.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
void wrapMPI(py::module& mod) {
auto mpi_mod = mod.def_submodule("mpi");
py::class_(mpi_mod, "sequential")
.def(py::init<>())
.def("__enter__", [](mpi::sequential& obj) { obj.enter(); })
.def("__exit__", [](mpi::sequential& obj, py::object, py::object,
py::object) { obj.exit(); });
mpi_mod.def("local_shape",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_size(global);
case 2:
return Partitioner<2>::local_size(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a, "Gives the local size of a 1D/2D global shape");
mpi_mod.def("global_shape",
[](std::vector local) {
switch (local.size()) {
case 1:
return Partitioner<1>::global_size(local);
case 2:
return Partitioner<2>::global_size(local);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"local_shape"_a, "Gives the global shape of a 1D/2D local shape");
mpi_mod.def("local_offset",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_offset(global);
case 2:
return Partitioner<2>::local_offset(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a,
"Gives the local offset of a 1D/2D global shape");
mpi_mod.def("size", []() { return mpi::size(); });
mpi_mod.def("rank", []() { return mpi::rank(); });
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/percolation.cpp b/python/wrap/percolation.cpp
index d9948db..a3d770b 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,86 +1,86 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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("__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 6fb149f..7a85807 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,171 +1,171 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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("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 2243ebe..4bfa259 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,202 +1,202 @@
/**
* @file
* @section LICENSE
*
- * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * 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("hurst", &Isopowerlaw::getHurst,
&Isopowerlaw::setHurst)
.def("rmsHeights", &Isopowerlaw::rmsHeights)
.def("moments", &Isopowerlaw::moments)
.def("alpha", &Isopowerlaw::alpha)
.def("rmsSlopes", &Isopowerlaw::rmsSlopes)
// legacy wrapper code
.def("getQ0", &Isopowerlaw::getQ0,
py::return_value_policy::reference)
.def("getQ1", &Isopowerlaw::getQ1,
py::return_value_policy::reference)
.def("getQ2", &Isopowerlaw::getQ2,
py::return_value_policy::reference)
.def("setQ0", &Isopowerlaw::setQ0,
py::return_value_policy::reference)
.def("setQ1", &Isopowerlaw::setQ1,
py::return_value_policy::reference)
.def("setQ2", &Isopowerlaw::setQ2,
py::return_value_policy::reference)
.def("setHurst", &Isopowerlaw::setHurst)
.def("getHurst", &Isopowerlaw