diff --git a/Jenkinsfile b/Jenkinsfile index d3b08fe..3293b72 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -1,120 +1,119 @@ pipeline { parameters {string(defaultValue: '', description: 'api-token', name: 'API_TOKEN') string(defaultValue: '', description: 'Token for readthedocs', name: 'RTD_TOKEN') string(defaultValue: '', description: 'buildable phid', name: 'BUILD_TARGET_PHID') string(defaultValue: '', description: 'Commit id', name: 'COMMIT_ID') string(defaultValue: '', description: 'Diff id', name: 'DIFF_ID') string(defaultValue: 'PHID-PROJ-gbo56hpf2y5bi7t5jusk', description: 'ID of the project', name: 'PROJECT_ID') } environment { PHABRICATOR_HOST = 'https://c4science.ch/api/' PYTHONPATH = sh returnStdout: true, script: 'echo ${WORKSPACE}/tests/ci/script/' } agent { dockerfile { additionalBuildArgs '--tag tamaas-environment'} } stages { stage('SCM Checkout') { steps { checkout scm: [ $class: 'GitSCM', branches: scm.branches, extensions: [[ $class: 'SubmoduleOption', recursiveSubmodules: true, ]], userRemoteConfigs: scm.userRemoteConfigs ] } } stage('Configure') { steps { sh '''#!/usr/bin/env bash echo "py_exec = \'python3\'" > build-setup.conf echo "build_python = \'true\'" >> build-setup.conf echo "build_tests = \'true\'" >> build-setup.conf echo "use_googletest = \'true\'" >> build-setup.conf - echo "legacy_bem = \'true\'" >> build-setup.conf echo "use_mpi = \'true\'" >> build-setup.conf''' } } stage('Compile') { steps { sh '''#!/usr/bin/env bash set -o pipefail rm -rf .sconf_temp .sconsign.dblite build-release scons | tee compilation.txt''' } post { failure { uploadArtifact('compilation.txt', 'Compilation') } } } stage('Run tests') { steps { sh 'PYTHONPATH=build-release/python/ python3 -m pytest --durations=0 --junitxml=results.xml build-release' } } } post { always { createArtifact("results.xml") junit 'results.xml' step([$class: 'XUnitBuilder', thresholds: [ [$class: 'SkippedThreshold', failureThreshold: '0'], [$class: 'FailedThreshold', failureThreshold: '0']], tools: [[$class: 'GoogleTestType', pattern: 'build-release/tests/gtest_results.xml']]]) } success { trigger_rtd() passed_hbm() } failure { failed_hbm() emailext( body: '''${SCRIPT, template="groovy-html.template"}''', mimeType: 'text/html', subject: "[Jenkins] ${currentBuild.fullDisplayName} Failed", recipientProviders: [[$class: 'CulpritsRecipientProvider']], to: 'lucas.frerot@protonmail.com', attachLog: true, compressLog: true) } } } def createArtifact(filename) { sh "./tests/ci/scripts/hbm send-uri -k 'Jenkins URI' -u ${BUILD_URL} -l 'View Jenkins result'" sh "./tests/ci/scripts/hbm send-junit-results -f ${filename}" } def uploadArtifact(artifact, name) { sh "./test/ci/scripts/hbm upload-file -f ${artifact} -n \"${name}\" -v ${PROJECT_ID}" } def failed_hbm() { sh "./tests/ci/scripts/hbm failed" } def passed_hbm() { sh "./tests/ci/scripts/hbm passed" } def trigger_rtd() { sh """ set -x curl -X POST -d "token=${RTD_TOKEN}" https://readthedocs.org/api/v2/webhook/tamaas/106141/ """ } diff --git a/SConstruct b/SConstruct index 3c025be..2a560b7 100644 --- a/SConstruct +++ b/SConstruct @@ -1,477 +1,469 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . # ------------------------------------------------------------------------------ # Imports # ------------------------------------------------------------------------------ from __future__ import print_function import os from subprocess import check_output import SCons from os.path import join, abspath, basename # Import below not strictly necessary, but good for pep8 from SCons.Script import ( EnsurePythonVersion, EnsureSConsVersion, Help, Environment, Variables, EnumVariable, PathVariable, BoolVariable, Split, SConscript, Export, Dir ) from version import get_git_subst from detect import ( FindFFTW, FindBoost, FindThrust, FindCuda, FindExpolit, FindPybind11 ) # ------------------------------------------------------------------------------ EnsurePythonVersion(2, 7) EnsureSConsVersion(2, 4) # ------------------------------------------------------------------------------ tamaas_info = dict( version="2.2.0", authors=[ u'Lucas Frérot', 'Guillaume Anciaux', 'Valentine Rey', 'Son Pham-Ba', u'Jean-François Molinari' ], maintainer=u'Lucas Frérot', email='lucas.frerot@protonmail.com', copyright=u"Copyright (©) 2016-2021 EPFL " + u"(École Polytechnique Fédérale de Lausanne), " + u"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" ) # ------------------------------------------------------------------------------ def detect_dependencies(env): """Detect all dependencies""" fftw_comp = { 'omp': ['omp'], 'cpp': [], 'tbb': ['threads'] } fftw_components = fftw_comp[env['backend']] if main_env['use_mpi']: fftw_components.append('mpi') FindFFTW(env, fftw_components, precision=env['real_type']) FindBoost(env, ['boost/preprocessor/seq.hpp']) FindExpolit(env) thrust_var = 'THRUST_ROOT' # Take cuda version of thrust if available if 'CUDA_ROOT' in env['ENV']: thrust_var = 'CUDA_ROOT' FindThrust(env, env['backend'], thrust_var) # Activate cuda if needed if env['backend'] == 'cuda': FindCuda(env) if env['build_python']: FindPybind11(env) # ------------------------------------------------------------------------------ # Main compilation # ------------------------------------------------------------------------------ # Compilation colors colors = { 'cyan': '\033[96m', 'purple': '\033[95m', 'blue': '\033[94m', 'green': '\033[92m', 'yellow': '\033[93m', 'gray': '\033[38;5;8m', 'orange': '\033[38;5;208m', 'red': '\033[91m', 'end': '\033[0m' } # Inherit all environment variables (for CXX detection, etc.) main_env = Environment( ENV=os.environ, ) # Set tamaas information for k, v in tamaas_info.items(): main_env[k] = v main_env['COLOR_DICT'] = colors # Build variables vars = Variables('build-setup.conf') vars.AddVariables( EnumVariable('build_type', 'Build type', 'release', allowed_values=('release', 'profiling', 'debug'), ignorecase=2), EnumVariable('backend', 'Thrust backend', 'omp', allowed_values=('cpp', 'omp', 'tbb'), ignorecase=2), EnumVariable('sanitizer', 'Sanitizer type', 'none', allowed_values=('none', 'memory', 'leaks', 'address'), ignorecase=2), PathVariable('prefix', 'Prefix where to install', '/usr/local'), # Dependencies paths PathVariable('FFTW_ROOT', 'FFTW custom path', os.getenv('FFTW_ROOT', ''), PathVariable.PathAccept), PathVariable('THRUST_ROOT', 'Thrust custom path', os.getenv('THRUST_ROOT', ''), PathVariable.PathAccept), PathVariable('BOOST_ROOT', 'Boost custom path', os.getenv('BOOST_ROOT', ''), PathVariable.PathAccept), PathVariable('CUDA_ROOT', 'Cuda custom path', os.getenv('CUDA_ROOT', ''), PathVariable.PathAccept), # Dependencies provided as submodule get different default PathVariable('GTEST_ROOT', 'Googletest custom path', os.getenv('GTEST_ROOT', '#third-party/googletest/googletest'), PathVariable.PathAccept), PathVariable('PYBIND11_ROOT', 'Pybind11 custom path', os.getenv('PYBIND11_ROOT', '#third-party/pybind11/include'), PathVariable.PathAccept), PathVariable('EXPOLIT_ROOT', 'Expolit custom path', os.getenv('EXPOLIT_ROOT', '#third-party/expolit/include'), PathVariable.PathAccept), # Executables ('CXX', 'Compiler', os.getenv('CXX', 'g++')), ('MPICXX', 'MPI Compiler wrapper', os.getenv('MPICXX', 'mpicxx')), ('py_exec', 'Python executable', 'python3'), # Compiler flags ('CXXFLAGS', 'C++ compiler flags', os.getenv('CXXFLAGS', "")), - # Compile legacy code - BoolVariable('legacy_bem', 'Compile legacy BEM code', False), # Cosmetic BoolVariable('verbose', 'Activate verbosity', False), BoolVariable('color', 'Color the non-verbose compilation output', False), # Tamaas components BoolVariable('build_doc', 'Build documentation', False), BoolVariable('build_tests', 'Build test suite', False), BoolVariable('build_python', 'Build python wrapper', True), # Dependencies BoolVariable('use_googletest', 'Build tests using GTest', False), BoolVariable('use_mpi', 'Builds multi-process parallelism', False), # Distribution options BoolVariable('strip_info', 'Strip binary of added information', False), BoolVariable('build_static_lib', "Build a static libTamaas", False), # Type variables EnumVariable('real_type', 'Type for real precision variables', 'double', allowed_values=('double', 'long double')), EnumVariable('integer_type', 'Type for integer variables', 'int', allowed_values=('int', 'long')), ) # Set variables of environment vars.Update(main_env) help_text = vars.GenerateHelpText(main_env) help_text += """ Commands: scons [build] [options]... Compile Tamaas (and additional modules/tests) scons install [prefix=/your/prefix] [options]... Install Tamaas to prefix scons dev Install symlink to Tamaas python module (useful to development purposes) scons test Run tests with pytest scons doc Compile documentation with Doxygen and Sphinx+Breathe scons archive Create a gzipped archive from source """ # noqa Help(help_text) # Save all options, not just those that differ from default with open('build-setup.conf', 'w') as setup: for option in vars.options: setup.write("# " + option.help + "\n") setup.write("{} = '{}'\n".format(option.key, main_env[option.key])) main_env['should_configure'] = \ not main_env.GetOption('clean') and not main_env.GetOption('help') build_type = main_env['build_type'] build_dir = 'build-' + main_env['build_type'] main_env['build_dir'] = build_dir if main_env['should_configure']: print("Building in " + build_dir) verbose = main_env['verbose'] # Remove colors if not set if not main_env['color']: for key in colors: colors[key] = '' if not verbose: main_env['CXXCOMSTR'] = main_env['SHCXXCOMSTR'] = \ u'{0}[Compiling ($SHCXX)] {1}$SOURCE'.format(colors['green'], colors['end']) main_env['LINKCOMSTR'] = main_env['SHLINKCOMSTR'] = \ u'{0}[Linking] {1}$TARGET'.format(colors['purple'], colors['end']) main_env['ARCOMSTR'] = u'{}[Ar]{} $TARGET'.format(colors['purple'], colors['end']) main_env['RANLIBCOMSTR'] = \ u'{}[Randlib]{} $TARGET'.format(colors['purple'], colors['end']) main_env['PRINT_CMD_LINE_FUNC'] = pretty_cmd_print # Include paths main_env.AppendUnique(CPPPATH=['#/src', '#/src/core', '#/src/mpi', '#/src/bem', '#/src/surface', '#/src/python', '#/src/percolation', '#/src/model', '#/src/model/elasto_plastic', '#/src/solvers', '#/src/gpu', '#/python']) # Changing the shared object extension main_env['SHOBJSUFFIX'] = '.o' # Back to gcc if cuda is activated if main_env['backend'] == "cuda" and "g++" not in main_env['CXX']: raise SCons.Errors.StopError('GCC should be used when compiling with CUDA') # OpenMP flags - compiler dependent omp_flags = { "g++": ["-fopenmp"], "clang++": ["-fopenmp"], "icpc": ["-qopenmp"] } def cxx_alias(cxx): for k in omp_flags.keys(): if k in cxx: return k raise SCons.Errors.StopError('Unsupported compiler: ' + cxx) cxx = cxx_alias(main_env['CXX']) # Setting main compilation flags main_env['CXXFLAGS'] = Split(main_env['CXXFLAGS']) main_env.AppendUnique( CXXFLAGS=Split('-std=c++14 -Wall -Wextra -pedantic'), CPPDEFINES=['TAMAAS_BACKEND=TAMAAS_BACKEND_{}'.format( main_env['backend'].upper() )], ) # Adding OpenMP flags if main_env['backend'] == 'omp': main_env.AppendUnique(CXXFLAGS=omp_flags[cxx]) main_env.AppendUnique(LINKFLAGS=omp_flags[cxx]) else: main_env.AppendUnique(CXXFLAGS=['-Wno-unknown-pragmas']) # Correct bug in clang? if main_env['backend'] == 'omp' and cxx == "clang++": main_env.AppendUnique(LIBS=["atomic"]) elif main_env['backend'] == 'tbb': main_env.AppendUnique(LIBS=['tbb']) -# Force deactivate legacy when using intel -# Intel does not like old openmp loops -if cxx == 'icpc' and main_env['legacy_bem']: - print("Using intel compiler => deactivating legacy code") - main_env['legacy_bem'] = False - # Manage MPI compiler if main_env['use_mpi']: main_env['CXX'] = '$MPICXX' main_env.AppendUnique(CPPDEFINES=['TAMAAS_USE_MPI']) main_env.AppendUnique(CXXFLAGS=['-Wno-cast-function-type']) # Flags and options if main_env['build_type'] == 'debug': main_env.AppendUnique(CPPDEFINES=['TAMAAS_DEBUG']) # Define the scalar types main_env.AppendUnique(CPPDEFINES={'TAMAAS_REAL_TYPE': '${real_type}', 'TAMAAS_INT_TYPE': '${integer_type}'}) # Compilation flags cxxflags_dict = { "debug": Split("-g -O0"), "profiling": Split("-g -O3 -fno-omit-frame-pointer"), "release": Split("-O3") } if main_env['sanitizer'] != 'none': if main_env['backend'] == 'cuda': raise SCons.Errors.StopError( "Sanitizers with cuda are not yet supported!") cxxflags_dict[build_type].append('-fsanitize=${sanitizer}') main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type]) main_env.AppendUnique(SHLINKFLAGS=cxxflags_dict[build_type]) main_env.AppendUnique(LINKFLAGS=cxxflags_dict[build_type]) if main_env['should_configure']: detect_dependencies(main_env) # Setting up the python name with version if main_env['build_python']: args = (main_env.subst("${py_exec} -c").split() + ["from distutils.sysconfig import get_python_version;" "print(get_python_version())"]) main_env['py_version'] = bytes(check_output(args)).decode() # Writing information file main_env.Tool('textfile') main_env['SUBST_DICT'] = get_git_subst() # Empty values if requested if main_env['strip_info']: for k in main_env['SUBST_DICT']: main_env['SUBST_DICT'][k] = "" # Substitution of environment file main_env['SUBST_DICT'].update({ '@build_type@': '$build_type', '@build_dir@': abspath(build_dir), '@build_version@': '$version', }) # Environment file content env_content = """export PYTHONPATH=@build_dir@/python:$$PYTHONPATH export LD_LIBRARY_PATH=@build_dir@/src:$$LD_LIBRARY_PATH """ # Writing environment file env_file = main_env.Textfile(join(build_dir, 'tamaas_environment.sh'), env_content) # Building sub-directories def subdir(dir): return SConscript(join(dir, 'SConscript'), variant_dir=join(build_dir, dir), duplicate=True) # Building Tamaas library Export('main_env') subdir('src') build_targets = ['build-cpp', env_file] install_targets = ['install-lib'] # Building Tamaas extra components for dir in ['python', 'tests']: if main_env['build_{}'.format(dir)] and not main_env.GetOption('help'): subdir(dir) build_targets.append('build-{}'.format(dir)) # Building API + Sphinx documentation if requested if main_env['build_doc']: doc_env = main_env.Clone() doc = doc_env.Command('#.phony_doc', '', 'make -C {doc} clean && make -C {doc}' .format(doc=Dir('#/doc'))) if doc_env['build_python']: doc_env.PrependENVPath('PYTHONPATH', doc_env.subst('../${build_dir}/python')) doc_env.Depends(doc, 'build-python') main_env.Alias('doc', doc) else: dummy_command(main_env, 'doc', 'Command "doc" does not do anything' ' without documentation activated ("build_doc=True")') # Define dummy dev command when python is deactivated if not main_env['build_python']: dummy_command(main_env, 'dev', 'Command "dev" does not do anything' + ' without python activated ("build_python=True")') else: install_targets.append('install-python') # Define dummy test command when tests are deactivated if not main_env['build_tests']: dummy_command(main_env, 'test', 'Command "test" does not do anything' + ' without tests activated ("build_tests=True")') # Definition of target aliases, a.k.a. sub-commands main_env.Alias('build', build_targets) # Define proper install targets main_env.Alias('install', install_targets) # Default target is to build stuff main_env.Default('build') # Building a tar archive archive = main_env.Command( 'tamaas-{}.tar.gz'.format(main_env['version']), '', ('tar --exclude-vcs --exclude-vcs-ignores ' '--exclude=third-party/googletest ' '--exclude=third-party/pybind11 ' '--exclude=joss ' '--exclude=".*" ' '-czf $TARGET {}'.format(basename(os.getcwd()))), chdir='..', ) main_env.Alias('archive', archive) diff --git a/doc/sphinx/source/contact.rst b/doc/sphinx/source/contact.rst index c72b05d..70dde2f 100644 --- a/doc/sphinx/source/contact.rst +++ b/doc/sphinx/source/contact.rst @@ -1,137 +1,137 @@ Solving contact =============== The resolution of contact problems is done with classes that inherit from :cpp:class:`tamaas::ContactSolver`. These usually take as argument a :cpp:class:`tamaas::Model` object, a surface described by a :cpp:class:`tamaas::Grid` or a 2D numpy array, and a tolerance. We will see the specificities of the different solver objects below. Elastic contact --------------- The most common case is normal elastic contact, and is most efficiently solved with :cpp:class:`tamaas::PolonskyKeerRey`. The advantage of this class is that it combines two algorithms into one. By default, it considers that the contact pressure field is the unknown, and tries to minimize the complementary energy of the system under the constraint that the mean pressure should be equal to the value supplied by the user, for example:: # ... solver = tm.PolonskyKeerRey(model, surface, 1e-12) solver.solve(1e-2) Here the average pressure is ``1e-2``. The solver can also be instanciated by specifying the the constraint should be on the mean gap instead of mean pressure:: solver = tm.PolonskyKeerRey(model, surface, 1e-12, constraint_type=tm.PolonskyKeerRey.gap) solver.solve(1e-2) The contact problem is now solved for a mean gap of ``1e-2``. Note that the choice of constraint affects the performance of the algorithm. Contact with adhesion --------------------- The second algorithm hidden in :cpp:class:`tamaas::PolonskyKeerRey` considers the **gap** to be the unknown. The constaint on the mean value can be chosen as above:: solver = tm.PolonskyKeerRey(model, surface, 1e-12, primal_type=tm.PolonskyKeerRey.gap, constraint_type=tm.PolonskyKeerRey.gap) solver.solve(1e-2) The advantage of this formulation is to be able to solve adhesion problems (`Rey et al. `_). This is done by adding a term to the potential energy functional that the solver tries to minimize:: adhesion_params = { "rho": 2e-3, "surface_energy": 2e-5 } adhesion = tm.ExponentialAdhesionFunctional(surface) adhesion.setParameters(adhesion) solver.addFunctionalterm(adhesion) solver.solve(1e-2) Custom classes can be used in place of the example term here. One has to inherit from :cpp:class:`tamaas::Functional`:: class AdhesionPython(tm.Functional): """ Functional class that extends a C++ class and implements the virtual methods """ def __init__(self, rho, gamma): super().__init__(self) self.rho = rho self.gamma = gamma def computeF(self, gap, pressure): return -self.gamma * np.sum(np.exp(-gap / self.rho)) def computeGradF(self, gap, gradient): gradient += self.gamma * np.exp(-gap / self.rho) / self.rho This example is actually equivalent to :cpp:class:`tamaas::functional::ExponentialAdhesionFunctional`. Tangential contact ------------------ For frictional contact, several solver classes are available. Among them, :cpp:class:`tamaas::Condat` is able to solve a coupled normal/tangential contact problem regardless of the material properties. It however solves an associated version of the Coulomb friction law. In general, the Coulomb friction used in contact makes the problem ill-posed. Solving a tangential contact problem is not much different from normal contact:: mu = 0.3 # Friction coefficient solver = tm.Condat(model, surface, 1e-12, mu) - solver.setMaxIterations(5000) # The default of 1000 may be too little + solver.max_iter = 5000 # The default of 1000 may be too little solver.solve([1e-2, 0, 1e-2]) # 3D components of applied mean pressure Elasto-plastic contact ---------------------- For elastic-plastic contact, one needs three different solvers: an elastic contact solver like the ones described above, a non-linear solver and a coupling solver. The non-linear solvers available in Tamaas are implemented in python and inherit from the C++ class :cpp:class:`tamaas::EPSolver`. They make use of the non-linear solvers available in scipy: :py:class:`DFSANESolver ` Implements an interface to :py:func:`scipy.optimize.root` wiht the DFSANE method. :py:class:`NewtonKrylovSolver ` Implements an interface to :py:func:`scipy.optimize.newton_krylov`. These solvers require a residual vector to cancel, the creation of which is handled with :cpp:class:`tamaas::ModelFactory`. Then, an instance of :cpp:class:`tamaas::EPICSolver` is needed to couple the contact and non-linear solvers for the elastic-plastic contact problem:: import tamaas as tm from tamaas.nonlinear_solvers import DFSANESolver # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 51, 51] # Order: [z, x, y] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Setup for plasticity residual = tm.ModelFactory.createResidual(model, sigma_y=0.1 * model.E, hardening=0.01 * model.E) epsolver = DFSANESolver(residual) # ... csolver = tm.PolonskyKeerRey(model, surface, 1e-12) epic = tm.EPICSolver(csolver, epsolver, 1e-7, relaxation=0.3) epic.solve(1e-3) # or use an accelerated scheme (which can sometimes not converge) epic.acceleratedSolve(1e-3) By default, :cpp:func:`tamaas::EPICSolver::solve` uses a relaxed fixed point. It can be tricky to make it converge: you need to decrease the relaxation parameter passed as argument of the constructor, but this also hinders the convergence rate. The function :cpp:func:`tamaas::EPICSolver::acceleratedSolve` does not require the tweaking of a relaxation parameter, so it can be faster if the latter does not have an optimal value. However, it is not guaranteed to converge. Finally, during the first iterations, the fixed point error will be large compared to the error of the non-linear solver. It can improve performance to have the tolerance of the non-linear solver (which is the most expensive step of the fixed point solve) decrease over the iterations. This can be achieved with the decorator class :py:class:`ToleranceManager `:: from tamaas.nonlinear_solvers import ToleranceManager, DFSANESolver @ToleranceManager(start=1e-2, end=1e-9, rate=1 / 3) class Solver(DFSANESolver): pass # ... epsolver = Solver(residual) # or epsolver = ToleranceManager(1e-2, 1e-9, 1/3)(DFSANESolver)(residual) diff --git a/doc/sphinx/source/model.rst b/doc/sphinx/source/model.rst index 4be5543..9ffbd20 100644 --- a/doc/sphinx/source/model.rst +++ b/doc/sphinx/source/model.rst @@ -1,103 +1,156 @@ Model and integral operators ============================ The class :cpp:class:`tamaas::Model` (and its counterpart :py:class:`Model `) is both a central class in Tamaas and one of the simplest. It mostly serves as holder for material properties, fields and integral operators, and apart from a linear elastic behavior does not perform any computation on its own. Model types ----------- :cpp:class:`tamaas::Model` has a concrete subclass :cpp:class:`tamaas::ModelTemplate` which implements the model function for a given :cpp:type:`tamaas::model_type`: :cpp:enumerator:`tamaas::basic_2d` Model type used in normal frictionless contact: traction and displacement are 2D fields with only one component. :cpp:enumerator:`tamaas::surface_2d` Model type used in frictional contact: traction and displacement are 2D fields with three components. :cpp:enumerator:`tamaas::volume_2d` Model type used in elastoplastic contact: tractions are the same as with :cpp:enumerator:`tamaas::surface_2d` but the displacement is a 3D field. The enumeration values suffixed with ``_1d`` are the one dimensional (line contact) counterparts of the above model types. The domain physical dimension and number of components are encoded in the class :cpp:class:`tamaas::model_type_traits`. Model creation and basic functionality -------------------------------------- The instanciation of a :py:class:`Model ` is done with the :py:class:`ModelFactory ` class and its :py:func:`createModel ` function:: physical_size = [1., 1.] discretization = [512, 512] model = tm.ModelFactory.createModel(tm.model_type.basic_2d, physical_size, discretization) .. warning:: For models of type ``volume_*d``, the first component of the ``physical_size`` and ``discretization`` arrays corresponds to the depth dimension (:math:`z` in most cases). For example:: tm.ModelFactory.createModel(tm.model_type.basic_2d, [0.3, 1, 1], [64, 81, 81]) creates a model of depth 0.3 and surface size 1\ :superscript:`2`, while the number of points is 64 in depth and 81\ :sup:`2` on the surface. This is done for data contiguity reasons, as we do discrete Fourier transforms in the horizontal plane. .. note:: If ran in an MPI context, the method :py:meth:`createModel ` expects the *global* system sizes and discretization of the model. The properties ``E`` and ``nu`` can be used to set the Young's modulus and Poisson ratio respectively:: model.E = 1 model.nu = 0.3 Fields can be easlily accessed with the ``[]`` operator, similar to Python's dictionaries:: surface_traction = model['traction'] + # surface_traction is a numpy array -To know what fields are available, you can call the :py:func:`getFields ` method. You can add new fields to a model object with :py:func:`registerField `, which is convenient for dumping. A model can also be used to compute stresses from a strain field:: +To know what fields are available, you can call the :py:class:`list` function on +a model (``list(model)``). You can add new fields to a model object with the +``[]`` operator:``model['new_field'] = new_field``, which is convenient for +dumping fields that are computed outside of Tamaas. + +.. note:: + The fields ``traction`` and ``displacement`` are always registered in models, + and are accessible via :py:attr:`model.traction + ` and :py:attr:`model.displacement + `. + +A model can also be used to compute stresses from a strain field:: import numpy as np - strain = np.zeros(model.getDiscretization() + [6]) # Mandel--Voigt notation + strain = np.zeros(model.shape + [6]) # Mandel--Voigt notation stress = np.zeros_like(strain) model.applyElasticity(stress, strain) +.. tip:: + ``print(model)`` gives a lot of information about the model: the model type, + shape, registered fields, and more! + + +Integral operators +------------------ + +Integral operators are a central part of Tamaas: they are carefully designed for +performance in periodic system. When a :py:class:`Model ` +object is used with contact solvers or with a residual object (for plasticty), +the objects using the model register integral operators with the model, so the +user typically does not have to worry about creating integral operators. + +Integral operators are accessed through the :py:attr:`operators +` property of a model object. The ``[]`` +operator gives access to the operators, and ``list(model.operators)`` gives the +list of registered operators:: + + # Accessing operator + elasticity = model.operators['hooke'] + + # Applying operator + elasticity(strain, stress) + + # Print all registered operators + print(list(model.operators)) + +.. note:: + At model creation, these operators are automatically registered: + + - ``hooke``: Hooke's elasticity law + - ``von_mises``: computes Von Mises stresses + - ``deviatoric``: computes the deviatoric part of a stress tensor + - ``eigenvalues``: computes the eigenvalues of a symetric tensor field + + :cpp:class:`Westergaard ` operators are automatically + registered when :py:meth:`solveNeumann + ` or :py:meth:`solveDirichlet + ` are called. + Model dumpers ------------- The submodule `tamaas.dumpers` contains a number of classes to save model data into different formats: :py:class:`UVWDumper ` Dumps a model to `VTK `_ format. Requires the `UVW `_ python package which you can install with pip:: pip install uvw This dumper is made for visualization with VTK based software like `Paraview `_. :py:class:`NumpyDumper ` Dumps a model to a compressed Numpy file. :py:class:`H5Dumper ` Dumps a model to a compressed `HDF5 `_ file. Requires the `h5py `_ package. The dumpers are initialzed with a basename and the fields that you wish to write to file (optionally you can set ``all_fields`` to ``True`` to dump all fields in the model). By default, each write operation creates a new file in a separate directory (e.g. :py:class:`UVWDumper ` creates a ``paraview`` directory). To write to a specific file you can use the `dump_to_file` method. Here is a usage example:: from tamaas.dumpers import UVWDumper, H5Dumper # Create dumper uvw_dumper = UVWDumper('rough_contact_example', 'stress', 'plastic_strain') # Dump model uvw_dumper << model # Or alternatively model.addDumper(H5Dumper('rough_contact_archive', all_fields=True)) model.addDumper(uvw_dumper) model.dump() The last ``model.dump()`` call will call both dumpers. The resulting files will have the following hierachy:: ./paraview/rough_contact_example_0000.vtr ./paraview/rough_contact_example_0001.vtr ./hdf5/rough_contact_archive_0000.h5 - + .. note:: Currently, only :py:class:`H5Dumper ` supports parallel output with MPI. diff --git a/doc/sphinx/source/rough_surfaces.rst b/doc/sphinx/source/rough_surfaces.rst index 7a4a802..fd82e80 100644 --- a/doc/sphinx/source/rough_surfaces.rst +++ b/doc/sphinx/source/rough_surfaces.rst @@ -1,109 +1,111 @@ Random rough surfaces ===================== The generation of stochatsticly rough surfaces is controlled in Tamaas by two abstract classes: :cpp:class:`tamaas::SurfaceGenerator` and :cpp:class:`tamaas::Filter`. The former provides access lets the user set the surface sizes and random seed, while the latter encodes the information of the spectrum of the surface. Two surface generation methods are provided: - :cpp:class:`tamaas::SurfaceGeneratorFilter` implements a Fourier filtering algorithm (see `Hu & Tonder `_), - :cpp:class:`tamaas::SurfaceGeneratorRandomPhase` implements a random phase filter. Both of these rely on a :cpp:class:`tamaas::Filter` object to provided the filtering information (usually power spectrum density coefficients). Tamaas provides two such classes and allows for Python subclassing: - :cpp:class:`tamaas::Isopowerlaw` provides a roll-off powerlaw, - :cpp:class:`tamaas::RegularizedPowerlaw` provides a powerlaw with a regularized rolloff. Tamaas also provided routines for surface statistics. Generating a surface -------------------- Let us now see how to generate a surface. Frist create a filter object and set the surface sizes:: import tamaas as tm # Create spectrum object spectrum = tm.Isopowerlaw2D() # Set spectrum parameters spectrum.q0 = 4 spectrum.q1 = 4 spectrum.q2 = 32 spectrum.hurst = 0.8 The ``spectrum`` object can be queried for information, such as the root-mean-square of heights, the various statistical moments, the spectrum bandwidth, etc. Then we create a generator object and build the random surface:: - generator = tm.SurfaceGeneratorFilter2D() - generator.setSizes([128, 128]) - generator.setSpectrum(spectrum) + generator = tm.SurfaceGeneratorFilter2D([128, 128]) + generator.spectrum = spectrum generator.random_seed = 0 surface = generator.buildSurface() .. important:: The ``surface`` object is a :py:class:`numpy.ndarray` wrapped around internal memory in the ``generator`` object, so a subsequent call to :py:func:`buildSurface ` may change its content. Note that if ``generator`` goes out of scope its memory will not be freed if there is still a live reference to the surface data. -.. note:: - If ran in an MPI context, the method :py:meth:`setSizes - ` expects the *global* shape of - the surface. +.. important:: + If ran in an MPI context, the constructor of + :py:class:`SurfaceGeneratorFilter2D + ` (and others) expects the *global* + shape of the surface. The shape can also be changed with ``generator.shape = + [64, 64]``. Custom filter ------------- Tamaas provides several classes that can be derived directly with Python classes, and :cpp:class:`tamaas::Filter` is one of them. Since it provides a single pure virtual method :cpp:func:`computeFilter `, it is easy to write a sub-class. Here we implement a class that takes a user-defined auto-correlation function and implements the :cpp:func:`computeFilter ` virtual function:: class AutocorrelationFilter(tm.Filter2D): def __init__(self, autocorrelation): tm.Filter2D.__init__(self) self.autocorrelation = autocorrelation.copy() def computeFilter(self, filter_coefficients): shifted_ac = np.fft.ifftshift(self.autocorrelation) filter_coefficients[...] = np.sqrt(np.fft.rfft2(shifted_ac)) Here ``filter_coefficients`` is also a :py:class:`numpy.ndarray` and is therefore easily manipulated. The creation of the surface then follows the same pattern as previously:: # Create spectrum object autocorrelation = ... # set your desired autocorrelation spectrum = AutocorrelationFilter(autocorrelation) generator = tm.SurfaceGenerator2D() - generator.setSpectrum(spectrum) + generator.shape = autocorrelation.shape + generator.spectrum = spectrum surface = generator.buildSurface() The lifetime of the ``spectrum`` object is associated to the ``generator`` when :py:func:`setSpectrum ` is called. Surface Satistics ----------------- Tamaas provides the C++ class :cpp:class:`tamaas::Statistics` and its wrapper :py:class:`Statistics2D ` to compute statistics on surfaces, including: - power spectrum density - autocorrelation - spectrum moments - root-mean-square of heights :math:`\sqrt{\langle h^2 \rangle}` - root-mean-square of slopes (computed in Fourier domain) :math:`\sqrt{\langle |\nabla h|^2\rangle}` All these quantities are computed in a discretization-independent manner: increasing the number of points in the surface should not drastically change the computed values (for a given spectrum). This allows to refine the discretization as much as possible to approximate a continuum. Note that the autocorrelation and PSD are fft-shifted. Here is a sample code plotting the PSD and autocorrelation:: psd = tm.Statistics2D.computePowerSpectrum(surface) psd = np.fft.fftshift(psd, axes=0) # Shifting only once axis because of R2C transform import matplotlib.pyplot as plt from matplotlib.colors import LogNorm plt.imshow(psd.real, norm=LogNorm()) acf = tm.Statistics2D.computeAutocorrelation(surface) acf = np.fft.fftshift(acf) plt.figure() plt.imshow(acf) plt.show() See ``examples/statistics.py`` for more usage examples of statistics. diff --git a/examples/adhesion.py b/examples/adhesion.py index dfc8ea4..fbc4973 100644 --- a/examples/adhesion.py +++ b/examples/adhesion.py @@ -1,127 +1,122 @@ #!/usr/bin/env python3 # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import tamaas as tm import matplotlib.pyplot as plt import numpy as np import argparse from matplotlib.colors import ListedColormap class AdhesionPython(tm.Functional): """ Functional class that extends a C++ class and implements the virtual methods """ def __init__(self, rho, gamma): tm.Functional.__init__(self) self.rho = rho self.gamma = gamma def computeF(self, gap, pressure): return -self.gamma * np.sum(np.exp(-gap / self.rho)) def computeGradF(self, gap, gradient): gradient += self.gamma * np.exp(-gap / self.rho) / self.rho parser = argparse.ArgumentParser() parser.add_argument('--local-functional', dest="py_adh", action="store_true", help="use the adhesion functional written in python") args = parser.parse_args() # Surface size n = 1024 # Surface generator -sg = tm.SurfaceGeneratorFilter2D() -sg.setSizes([n, n]) -sg.setRandomSeed(0) +sg = tm.SurfaceGeneratorFilter2D([n, n]) +sg.random_seed = 0 # Spectrum -spectrum = tm.Isopowerlaw2D() -sg.setFilter(spectrum) +sg.spectrum = tm.Isopowerlaw2D() # Parameters -spectrum.q0 = 16 -spectrum.q1 = 16 -spectrum.q2 = 64 -spectrum.hurst = 0.8 +sg.spectrum.q0 = 16 +sg.spectrum.q1 = 16 +sg.spectrum.q2 = 64 +sg.spectrum.hurst = 0.8 # Generating surface surface = sg.buildSurface() -# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface) surface /= n -# print(spectrum.rmsSlopes()) -# print(tm.SurfaceStatistics.computeRMSSlope(surface)) plt.imshow(surface) plt.title('Rough surface') # Creating model model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n]) # Solver solver = tm.PolonskyKeerRey(model, surface, 1e-12, tm.PolonskyKeerRey.gap, tm.PolonskyKeerRey.gap) adhesion_params = { "rho": 2e-3, "surface_energy": 2e-5 } # Use the python derived from C++ functional class if args.py_adh: adhesion = AdhesionPython(adhesion_params["rho"], adhesion_params["surface_energy"]) # Use the C++ class else: adhesion = tm.ExponentialAdhesionFunctional(surface) adhesion.setParameters(adhesion_params) solver.addFunctionalTerm(adhesion) # Solve for target pressure g_target = 5e-2 solver.solve(g_target) -tractions = model.getTraction() +tractions = model.traction plt.figure() plt.imshow(tractions) plt.colorbar() plt.title('Contact tractions') plt.figure() zones = np.zeros_like(tractions) tol = 1e-6 zones[tractions > tol] = 1 zones[tractions < -tol] = -1 plt.imshow(zones, cmap=ListedColormap(['white', 'gray', 'black'])) plt.colorbar(ticks=[-2/3, 0, 2/3]).set_ticklabels([ 'Adhesion', 'No Contact', 'Contact' ]) plt.title('Contact and Adhesion Zones') plt.show() diff --git a/examples/legacy/adhesion.py b/examples/legacy/adhesion.py deleted file mode 100644 index 64ed1d5..0000000 --- a/examples/legacy/adhesion.py +++ /dev/null @@ -1,161 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - -import sys -import tamaas as tm -import numpy as np -import matplotlib.pyplot as plt - -def plotSurface(surface): - fig = plt.figure() - ax = fig.add_subplot(111) - img = ax.imshow(surface) - #fig.colorbar(img) - -def constructHertzProfile(size, curvature): - radius = 1. / curvature - x = np.linspace(-0.5, 0.5, size) - y = np.linspace(-0.5, 0.5, size) - x, y = np.meshgrid(x, y) - surface = np.sqrt(radius**2 - x**2 - y**2) - surface -= surface.min() - return surface.copy() - -def computeHertzDisplacement(e_star, contact_size, max_pressure, size): - x = np.linspace(-0.5, 0.5, size) - y = np.linspace(-0.5, 0.5, size) - x, y = np.meshgrid(x, y) - disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - (x**2 + y**2)) - return disp.copy() - -def main(): - grid_size = 128 - curvature = 0.5 - effective_modulus = 1. - load = 0.1 - - surface_energy = 0 - rho = 2.071e-7 - - surface = constructHertzProfile(grid_size, curvature) - # SG = tm.SurfaceGeneratorFilterFFT() - # SG.getGridSize().assign(grid_size) - # SG.getHurst().assign(0.8) - # SG.getRMS().assign(0.002); - # SG.getQ0().assign(8); - # SG.getQ1().assign(8); - # SG.getQ2().assign(16); - # SG.getRandomSeed().assign(156); - # SG.Init() - # surface = SG.buildSurface() - print "Max height {}".format(surface.max()) - print "Min height {}".format(surface.min()) - - bem = tm.BemGigi(surface) - bem.setDumpFreq(1) - functional = tm.ExponentialAdhesionFunctional(bem) - functional.setParameter('rho', rho) - functional.setParameter('surface_energy', surface_energy) - bem.setEffectiveModulus(effective_modulus) - bem.addFunctional(functional) - bem.computeEquilibrium(1e-6, load) - - tractions = bem.getTractions() - print "Average pressure = {}".format(tractions.mean()) - # bem.computeTrueDisplacements() - t_displacements = bem.getTrueDisplacements() - t_gap = bem.getGap() - - plotSurface(tractions) - - plt.figure() - plt.plot(surface[grid_size/2, :]) - plt.title("Surface") - plt.figure() - plt.plot(tractions[grid_size/2, :]) - plt.title("Pressure") - plt.figure() - plt.plot(t_gap[grid_size/2, :]) - plt.title("Gap") - plt.figure() - plt.plot(t_displacements[grid_size/2, :]) - plt.title("Displacement") - plt.figure() - plt.plot(t_displacements[grid_size/2, :] - surface[grid_size/2, :]) - plt.title("Disp-surf") - plotSurface(t_displacements) - - plt.show() - return 0 - - # Testing contact area against Hertz solution for solids of revolution - contact_area = tm.SurfaceStatistics.computeContactArea(tractions) - hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.) - hertz_area = np.pi * hertz_contact_size**2 - area_error = np.abs(hertz_area - contact_area) / hertz_area - print "Area error: {}".format(area_error) - - # Testing maximum pressure - max_pressure = tractions.max() - hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi - pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure - print "Max pressure error: {}".format(pressure_error) - - # Testing displacements - hertz_displacements = computeHertzDisplacement(effective_modulus, - hertz_contact_size, - hertz_max_pressure, - grid_size) - - # Selecing only the points that are in contact - contact_indexes = [(i, j, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)] - contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes)) - - # Displacements of bem are centered around the mean of the whole surface - # and Hertz displacements are not centered, so we need to compute mean - # on the contact zone for both arrays - bem_mean = 0. - hertz_mean = 0. - for index in contact_indexes: - bem_mean += displacements[index] - hertz_mean += hertz_displacements[index] - - bem_mean /= len(contact_indexes) - hertz_mean /= len(contact_indexes) - # Correction applied when computing error - correction = hertz_mean - bem_mean - - # Computation of error of displacement in contact zone - error = 0. - hertz_norm = 0. - for index in contact_indexes: - error += (hertz_displacements[index] - displacements[index] - correction)**2 - hertz_norm += (hertz_displacements[index] - hertz_mean)**2 - - displacement_error = np.sqrt(error / hertz_norm) - print "Displacement error (in contact zone): {}".format(displacement_error) - - if area_error > 1e-2 or pressure_error > 1e-2 or displacement_error > 1e-4: - return 1 - - return 0 - -if __name__ == "__main__": - sys.exit(main()) diff --git a/examples/legacy/cluster_statistics.py b/examples/legacy/cluster_statistics.py deleted file mode 100644 index 8b1bca2..0000000 --- a/examples/legacy/cluster_statistics.py +++ /dev/null @@ -1,104 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - -import numpy as np - - -def plotSurface(surf): - fig = plt.figure() - axe = fig.add_subplot(111) - img = axe.imshow(surf) - cbar = fig.colorbar(img) - - -################################################################ -# surface generation -################################################################ - -from tamaas import * - -n = 128 -SG = SurfaceGeneratorFilterFFT() -SG.getGridSize().assign(n) -SG.getHurst().assign(0.8) -SG.getRMS().assign(1.); -SG.getQ0().assign(4); -SG.getQ1().assign(4); -SG.getQ2().assign(64); -SG.getRandomSeed().assign(20); -SG.Init() -s = SG.buildSurface() -rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(s) -s *= 1./rms_slopes_spectral -s -= s.max() - -################################################################ -# surface plot -################################################################ - -import matplotlib.pyplot as plt - -#fig1 = plotSurface(s) - -################################################################ -# contact contact -################################################################ - -bem = BemPolonski(s) -bem.setEffectiveModulus(1.) - -load = 0.1 -bem.computeEquilibrium(1e-13,load) - -tractions = bem.getTractions() -displacements = bem.getDisplacements() - -#fig2 = plotSurface(tractions) -#fig3 = plotSurface(displacements) - -################################################################ -# cluster statistics -################################################################ - - -area = ContactArea(tractions.shape[0],1.); -area.getSurface()[tractions > 0.] = 1 -# Cluster detection -print "===============================" -print "Detect contact clusters: " -area.detectContactClusters() -print "Found " + str(area.getNumClusters()) -clusters = ContactClusterCollection(area); -cluster_list = area.getContactClusters(); -cluster_areas = [float(c.getA())/float(n)**2 for c in cluster_list] -cluster_perimeters = [c.getP() for c in cluster_list] -print "cluster_areas",cluster_areas -print "cluster_perimeters",cluster_perimeters -print "cluster_total_area",clusters.getTotalArea() -print "cluster_total_perimeter",clusters.getTotalPerimeter() -print "nb_clusters_with_hole",clusters.getNbClustersWithHole() -print "nb_clusters",clusters.getNbClusters() -################################################################ -value,bins = np.histogram(cluster_areas,bins=50) -fig = plt.figure() -axe = fig.add_subplot(111) -axe.loglog(bins[:-1],value,'-o') -################################################################ -plt.show() diff --git a/examples/legacy/contact.py b/examples/legacy/contact.py deleted file mode 100644 index 0adbc84..0000000 --- a/examples/legacy/contact.py +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - - -def plotSurface(surf): - fig = plt.figure() - axe = fig.add_subplot(111) - img = axe.imshow(surf) - cbar = fig.colorbar(img) - - -################################################################ -# surface generation -################################################################ - -from tamaas import * - - -SG = SurfaceGeneratorFilterFFT() -SG.getGridSize().assign(128) -SG.getHurst().assign(0.8) -SG.getRMS().assign(1.); -SG.getQ0().assign(4); -SG.getQ1().assign(4); -SG.getQ2().assign(64); -SG.getRandomSeed().assign(20); -SG.Init() -s = SG.buildSurface() -rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(s) -s *= 1./rms_slopes_spectral -s -= s.max() - -################################################################ -# surface plot -################################################################ - -import matplotlib.pyplot as plt - -fig1 = plotSurface(s) - -################################################################ -# contact contact -################################################################ - -bem = BemPolonski(s) -bem.setEffectiveModulus(1.) - -load = 0.1 -bem.computeEquilibrium(1e-13,load) - -tractions = bem.getTractions() -displacements = bem.getDisplacements() - -fig2 = plotSurface(tractions) -fig3 = plotSurface(displacements) -plt.show() diff --git a/examples/legacy/fractal_surface.py b/examples/legacy/fractal_surface.py deleted file mode 100644 index 0a07a9c..0000000 --- a/examples/legacy/fractal_surface.py +++ /dev/null @@ -1,102 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - -from tamaas import * -import argparse - -parser = argparse.ArgumentParser() -parser.add_argument("--rescale", help="Rescale surface for RMS(slopes) = 1", action='store_true') -parser.add_argument("--N", help="Number of points", type=int, default=512) -parser.add_argument("--k0", help="Roll-off wave number", type=int, default=4) -parser.add_argument("--k1", help="Low cutoff wave number", type=int, default=4) -parser.add_argument("--k2", help="High cutoff wave number", type=int, default=32) -parser.add_argument("--rms", help="RMS(heights)", default=1.) -parser.add_argument("--H", help="Hurst exponent", default=0.8) - -args = parser.parse_args() - -#generate surface -SG = SurfaceGeneratorFilterFFT() -SG.getGridSize().assign(args.N) -SG.getHurst().assign(args.H) -SG.getRMS().assign(args.rms); -SG.getQ0().assign(args.k0); -SG.getQ1().assign(args.k1); -SG.getQ2().assign(args.k2); -SG.getRandomSeed().assign(156); -SG.Init() -a = SG.buildSurface() - -if args.rescale: - rms_slopes = SurfaceStatistics.computeSpectralRMSSlope(a) - a /= rms_slopes - -#compute and print surface statistics -class Stats: - - def __getitem__(self,key): - return self.__dict__[key] - - -stats = Stats() - -stats.size = SG.getGridSize() -stats.hurst = SG.getHurst().value() -stats.rms = SG.getRMS() -stats.k0 = SG.getQ0() -stats.k1 = SG.getQ1().value() -stats.k2 = SG.getQ2().value() -stats.seed = SG.getRandomSeed() -stats.rms_spectral = SurfaceStatistics.computeSpectralStdev(a); -stats.rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(a); -stats.rms_geometric = a.std(ddof=1) -stats.rms_slopes_geometric = SurfaceStatistics.computeRMSSlope(a); -stats.moments = SurfaceStatistics.computeMoments(a); -stats.m0 = stats['rms_spectral']**2 -stats.m2 = stats.moments[0] -stats.m4 = stats.moments[1] -stats.alpha = stats['m0']*stats['m4']/(stats['m2']**2) -stats.L = 1. -stats.m0prime = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,1. , stats.L) -stats.moment_A = stats.m0/stats.m0prime -stats.analytic_m0 = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L); -stats.analytic_m2 = SurfaceStatistics.computeAnalyticFractalMoment(2,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L); -stats.analytic_m4 = SurfaceStatistics.computeAnalyticFractalMoment(4,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L); -stats.analytic_alpha = stats.analytic_m0*stats.analytic_m4/(stats.analytic_m2*stats.analytic_m2); - -print """ -[N] {size} -[rms] {rms} -[rmsSpectral] {rms_spectral} -[rmsSlopeSpectral] {rms_slopes_spectral} -[rmsSlopeGeometric] {rms_slopes_geometric} -[Hurst] {hurst} -[k1] {k1} -[k2] {k2} -[moment A] {moment_A} -[m0] {m0} -[analytic m0] {analytic_m0} -[m2] {m2} -[analytic m2] {analytic_m2} -[m4] {m4} -[analytic m4] {analytic_m4} -[alpha] {alpha} -[analytic_alpha] {analytic_alpha} -[seed] {seed} diff --git a/examples/legacy/hertz.py b/examples/legacy/hertz.py deleted file mode 100644 index 2a30094..0000000 --- a/examples/legacy/hertz.py +++ /dev/null @@ -1,116 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - -import tamaas - - -def dumpHertzPrediction(file, bem, pressure, r, E ): - A = tamaas.SurfaceStatistics.computeContactArea(bem.getTractions()) - Aratio = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions()) - radius = np.sqrt(A/np.pi) - #L = bem.getSurface().getL() - L = 1. - load = pressure * L*L - radius_hertz = pow(0.75*load*r/E,1./3.) - p0_hertz = 1./np.pi*pow(6.*E*E*load/r/r,1./3.) - p0 = tamaas.SurfaceStatistics.computeMaximum(bem.getTractions()) - n = bem.getTractions().shape[0] - computed_load = bem.getTractions().sum()/n/n*L*L - file.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\t{9}\n".format(pressure,load,computed_load,Aratio,A,radius,radius_hertz,radius_hertz/radius,p0,p0_hertz)) - - -################################################################ - -def main(argv): - - parser = argparse.ArgumentParser(description='Hertz test for the Boundary element method of Stanley and Kato') - - parser.add_argument("--N",type=int, help="Surface size", required=True) - parser.add_argument("--r",type=float, help="radius of hertz sphere", required=True) - parser.add_argument("--steps",type=int, help="number of steps within which the pressure is applied.", required=True) - parser.add_argument("--precision",type=float, help="relative precision, convergence if | (f_i - f_{i-1})/f_i | < Precision.", required=True) - parser.add_argument("--e_star",type=float, help="effective elastic modulus" , required=True) - parser.add_argument("--L",type=float, help="size of the surface" , required=True) - parser.add_argument("--max_pressure",type=float, help="maximal load requested" , required=True) - parser.add_argument("--plot_surface",type=bool, help="request output of text files containing the contact pressures on the surface" , default=False) - parser.add_argument("--nthreads",type=int, help="request a number of threads to use via openmp to compute" , default=1) - - args = parser.parse_args() - arguments = vars(args) - - n = arguments["N"] - r = arguments["r"] - max_pressure = arguments["max_pressure"] - Ninc = arguments["steps"] - epsilon = arguments["precision"] - Estar = arguments["e_star"] - L = arguments["L"] - plot_surface = arguments["plot_surface"] - nthreads = arguments["nthreads"] - - pressure = 0. - dp = max_pressure/float(Ninc) - - s = np.zeros((n,n)) - print s.shape - for i in range(0,n): - for j in range(0,n): - x = 1.*i - n/2 - y = 1.*j - n/2 - d = (x*x + y*y)*1./n/n*L*L - - if d < r*r: s[i,j] = - r + np.sqrt(r*r-d) - else: s[i,j] = - r - - - print "\n::: DATA ::::::::::::::::::::::::::::::\n" - print " [N] {0}\n".format(n) - print " [r] {0}\n".format(r) - print " [Pext] {0}\n".format(max_pressure) - print " [Steps] {0}\n".format(Ninc) - print " [Precision] {0}\n".format(epsilon) - - bem = tamaas.BemPolonski(s) - bem.setEffectiveModulus(Estar) - - file = open("hertz-prediction",'w') - - file.write("#pressure\tload\tcomputed-load\tarea_ratio\tarea\tradius\tradius-hertz\tradius/radius-hertz\tp0\tp0-hertz\n") - - for i in range(0,Ninc): - - pressure += dp - - bem.computeEquilibrium(epsilon,pressure) - A = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions()) - - dumpHertzPrediction(file,bem,pressure,r,Estar) - - if A == 1.0: - file.close() - break - - tamaas.dumpTimes() - -################################################################ -import sys -import argparse -import numpy as np -main(sys.argv) diff --git a/examples/legacy/surface_1d.py b/examples/legacy/surface_1d.py deleted file mode 100644 index 6c48dca..0000000 --- a/examples/legacy/surface_1d.py +++ /dev/null @@ -1,83 +0,0 @@ -#!/usr/bin/env python -# @file -# @section LICENSE -# -# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), -# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as published -# by the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. -# -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see . - -from __future__ import division -import tamaas as tm -from numpy.fft import fft, fftshift, fftfreq -import matplotlib.pyplot as plt - -tm.initialize() - - -def draw_triangle(pos, w, ax, style, alpha=1.5, base=10): - h = alpha * w - ax.loglog([base**pos[0], base**(pos[0]+w)], - [base**pos[1], base**pos[1]], style) - ax.loglog([base**pos[0], base**pos[0]], - [base**pos[1], base**(pos[1]+h)], style) - ax.loglog([base**(pos[0]+w), base**pos[0]], - [base**pos[1], base**(pos[1]+h)], style) - # ax.text(base**(pos[0]+w/3), base**(pos[1]+h/3), "$-{}$".format(alpha), - # horizontalalignment='center', - # verticalalignment='center') - ax.text(base**(pos[0]+w/2), base**(pos[1]-0.2), '$1$', - horizontalalignment='center', - verticalalignment='top') - ax.text(base**(pos[0]-0.15), base**(pos[1]+h/2), '${}$'.format(alpha), - horizontalalignment='right', - verticalalignment='center') - - -iso = tm.Isopowerlaw1D() -helper = tm.ParamHelper(iso) - -params = { - "Q0": 16, - "Q1": 16, - "Q2": 512, - "Hurst": 0.8 -} - -helper.set(params) - -gen = tm.SurfaceGeneratorFilter1D() -gen.setFilter(iso) -gen.setSizes([2048]) -# gen.setRandomSeed(0) # uncomment for fixed seed - -surf = gen.buildSurface() - -plt.plot(surf) - -surf_fft = fft(surf) -psd = surf_fft*surf_fft.conj() -psd = fftshift(psd) - -plt.figure() -freqs = fftshift(fftfreq(psd.size, d=1/psd.size)) -plt.plot(freqs[psd.size/2+1:], psd.real[psd.size/2+1:]) -draw_triangle([5, -4], 2, plt.gca(), '-k', 2*(params["Hurst"]+1), 2) -plt.yscale('log') -plt.xscale('log', basex=2) -plt.ylim([10**(-2), 10**8]) - -plt.show() - -tm.finalize() diff --git a/examples/pipe_tools/surface b/examples/pipe_tools/surface index a0d77f7..18bdfe5 100755 --- a/examples/pipe_tools/surface +++ b/examples/pipe_tools/surface @@ -1,81 +1,78 @@ #!/usr/bin/env python3 # -*- mode: python; coding: utf-8 -*- # vim: set ft=python: """ Create a random self-affine rough surface from command line parameters. """ import argparse import sys import time import tamaas as tm import numpy as np __author__ = "Lucas Frérot" __copyright__ = ( "Copyright (©) 2019-2020, EPFL (École Polytechnique Fédérale de Lausanne)," "\nLaboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" ) __license__ = "AGPL" __email__ = "lucas.frerot@gmail.com" parser = argparse.ArgumentParser( description="Generate a self-affine rough surface" ) parser.add_argument("--cutoffs", "-K", nargs=3, type=int, help="Long, rolloff and short wavelength cutoffs", metavar=('k_l', 'k_r', 'k_s'), required=True) parser.add_argument("--sizes", nargs=2, type=int, help="Number of points", metavar=('nx', 'ny'), required=True) parser.add_argument("--hurst", "-H", type=float, help="Hurst exponent", required=True) parser.add_argument("--rms", type=float, help="Root-mean-square of slopes", default=1.) parser.add_argument("--seed", type=int, help="Random seed", default=int(time.time())) parser.add_argument("--generator", help="Generation method", choices=('random_phase', 'filter'), default='random_phase') parser.add_argument("--output", "-o", help="Output file name (compressed if .gz)") args = parser.parse_args() -spectrum = tm.Isopowerlaw2D() -spectrum.q0 = args.cutoffs[0] -spectrum.q1 = args.cutoffs[1] -spectrum.q2 = args.cutoffs[2] -spectrum.hurst = args.hurst - if args.generator == 'random_phase': - generator = tm.SurfaceGeneratorRandomPhase2D() + generator = tm.SurfaceGeneratorRandomPhase2D(args.sizes) elif args.generator == 'filter': - generator = tm.SurfaceGeneratorFilter2D() + generator = tm.SurfaceGeneratorFilter2D(args.sizes) else: raise ValueError('Unknown generator method {}'.format(args.generator)) -generator.setSizes(args.sizes) +generator.spectrum = tm.Isopowerlaw2D() +generator.spectrum.q0 = args.cutoffs[0] +generator.spectrum.q1 = args.cutoffs[1] +generator.spectrum.q2 = args.cutoffs[2] +generator.spectrum.hurst = args.hurst generator.random_seed = args.seed -generator.setSpectrum(spectrum) -surface = generator.buildSurface() / spectrum.rmsSlopes() * args.rms +surface = generator.buildSurface() / generator.spectrum.rmsSlopes() * args.rms output = args.output if args.output is not None else sys.stdout np.savetxt(output, surface) diff --git a/examples/plasticity.py b/examples/plasticity.py index 31f1fd4..139973e 100644 --- a/examples/plasticity.py +++ b/examples/plasticity.py @@ -1,86 +1,85 @@ #!/usr/bin/env python3 # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import numpy as np import tamaas as tm from tamaas.dumpers import H5Dumper as Dumper from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver, ToleranceManager -tm.initialize(1) - # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 51, 51] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model.E = 1. model.nu = 0.3 # Setup for plasticity residual = tm.ModelFactory.createResidual(model, sigma_y=0.1 * model.E, hardening=0.01 * model.E) # Possibly change integration method # residual.setIntegrationMethod(tm.integration_method.cutoff, 1e-12) # Setup non-linear solver with variable tolerance epsolver = ToleranceManager(1e-5, 1e-9, 1/4)(Solver)(residual) # Setup for contact x = np.linspace(0, system_size[1], discretization[1], endpoint=False, dtype=tm.dtype) y = np.linspace(0, system_size[2], discretization[2], endpoint=False, dtype=tm.dtype) xx, yy = np.meshgrid(x, y, indexing='ij') R = 0.2 surface = -((xx - flat_domain[0] / 2)**2 + (yy - flat_domain[1] / 2)**2) / (2 * R) +# Extract local part of the global surface local_shape = tm.mpi.local_shape(surface.shape) offset = tm.mpi.local_offset(surface.shape) local_surface = surface[offset:offset+local_shape[0], :].copy() csolver = tm.PolonskyKeerRey(model, local_surface, 1e-12, tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.pressure) # EPIC setup epic = tm.EPICSolver(csolver, epsolver, 1e-7) # Dumper dumper_helper = Dumper('hertz', 'displacement', 'stress', 'plastic_strain') model.addDumper(dumper_helper) loads = np.linspace(0.001, 0.005, 3) for i, load in enumerate(loads): epic.acceleratedSolve(load) model.dump() tm.Logger().get(tm.LogLevel.info) \ << "---> Solved load step {}/{}".format(i+1, len(loads)) diff --git a/examples/rough_contact.py b/examples/rough_contact.py index bdf967f..bc0ca84 100644 --- a/examples/rough_contact.py +++ b/examples/rough_contact.py @@ -1,71 +1,64 @@ #!/usr/bin/env python3 # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import tamaas as tm import matplotlib.pyplot as plt # Initialize threads and fftw tm.set_log_level(tm.LogLevel.info) # Show progression of solver # Surface size n = 512 # Surface generator -sg = tm.SurfaceGeneratorFilter2D() -sg.setSizes([n, n]) +sg = tm.SurfaceGeneratorFilter2D([n, n]) sg.random_seed = 0 # Spectrum -spectrum = tm.Isopowerlaw2D() -sg.setFilter(spectrum) +sg.spectrum = tm.Isopowerlaw2D() # Parameters -spectrum.q0 = 16 -spectrum.q1 = 16 -spectrum.q2 = 64 -spectrum.hurst = 0.8 +sg.spectrum.q0 = 16 +sg.spectrum.q1 = 16 +sg.spectrum.q2 = 64 +sg.spectrum.hurst = 0.8 # Generating surface surface = sg.buildSurface() -# surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface) -surface /= n -# print(spectrum.rmsSlopes()) -# print(tm.SurfaceStatistics.computeRMSSlope(surface)) +surface /= tm.Statistics2D.computeSpectralRMSSlope(surface) plt.imshow(surface) # Creating model model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n]) # Solver -solver = tm.PolonskyKeerRey(model, surface, 1e-12, - tm.PolonskyKeerRey.pressure, - tm.PolonskyKeerRey.pressure) +solver = tm.PolonskyKeerRey(model, surface, 1e-12) # Solve for target pressure p_target = 0.1 solver.solve(p_target) plt.figure() -plt.imshow(model.getTraction()) +plt.imshow(model.traction) plt.title('Contact tractions') -print(model.getTraction().mean()) +print(model.traction.mean()) plt.show() diff --git a/examples/saturated.py b/examples/saturated.py index 81f4a59..dfad863 100644 --- a/examples/saturated.py +++ b/examples/saturated.py @@ -1,67 +1,66 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import tamaas as tm import numpy as np import matplotlib.pyplot as plt grid_size = 256 load = 1.6 p_sat = 100 iso = tm.Isopowerlaw2D() iso.q0 = 4 iso.q1 = 4 iso.q2 = 16 iso.hurst = 0.8 -sg = tm.SurfaceGeneratorFilter2D() +sg = tm.SurfaceGeneratorFilter2D([grid_size, grid_size]) sg.random_seed = 2 -sg.setSizes([grid_size, grid_size]) -sg.setSpectrum(iso) +sg.spectrum = iso surface = sg.buildSurface() surface -= np.max(surface) model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [1., 1.], [grid_size, grid_size]) model.E = 1. model.nu = 0 esolver = tm.PolonskyKeerRey(model, surface, 1e-12, tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.pressure) esolver.solve(load) -elastic_tractions = model['traction'].copy() +elastic_tractions = model.traction.copy() plt.imshow(elastic_tractions) plt.title('Elastic contact tractions') plt.colorbar() -model['traction'][:] = 0. +model.traction[:] = 0. solver = tm.KatoSaturated(model, surface, 1e-12, p_sat) -solver.setDumpFrequency(1) -solver.setMaxIterations(200) +solver.dump_freq = 1 +solver.max_iter = 200 solver.solve(load) plt.figure() -plt.imshow(model['traction']) +plt.imshow(model.traction) plt.title('Saturated contact tractions') plt.colorbar() plt.show() diff --git a/examples/statistics.py b/examples/statistics.py index 798a591..f47b220 100644 --- a/examples/statistics.py +++ b/examples/statistics.py @@ -1,86 +1,83 @@ #!/usr/bin/env python3 # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import os import numpy as np import matplotlib.pyplot as plt import tamaas as tm from matplotlib.colors import LogNorm -# Initialize threads and fftw tm.set_log_level(tm.LogLevel.info) # Show progression of solver # Surface size n = 512 # Surface generator -sg = tm.SurfaceGeneratorFilter2D() -sg.setSizes([n, n]) +sg = tm.SurfaceGeneratorFilter2D([n, n]) sg.random_seed = 0 # Spectrum -spectrum = tm.Isopowerlaw2D() -sg.setFilter(spectrum) +sg.spectrum = tm.Isopowerlaw2D() # Parameters -spectrum.q0 = 16 -spectrum.q1 = 16 -spectrum.q2 = 64 -spectrum.hurst = 0.8 +sg.spectrum.q0 = 16 +sg.spectrum.q1 = 16 +sg.spectrum.q2 = 64 +sg.spectrum.hurst = 0.8 # Generating surface surface = sg.buildSurface() # Computing PSD and shifting for plot psd = tm.Statistics2D.computePowerSpectrum(surface) psd = np.fft.fftshift(psd, axes=0) plt.imshow(psd.real, norm=LogNorm()) plt.gca().set_title('Power Spectrum Density') plt.gcf().tight_layout() # Computing autocorrelation and shifting for plot acf = tm.Statistics2D.computeAutocorrelation(surface) acf = np.fft.fftshift(acf) plt.figure() plt.imshow(acf) plt.gca().set_title('Autocorrelation') plt.gcf().tight_layout() plt.show() # Write the rough surface in paraview's VTK format try: import uvw try: os.mkdir('paraview') except FileExistsError: pass x = np.linspace(0, 1, n, endpoint=True) y = x.copy() with uvw.RectilinearGrid(os.path.join('paraview', 'surface.vtr'), (x, y)) as grid: grid.addPointData(uvw.DataArray(surface, range(2), 'surface')) except ImportError: print("uvw not installed, will not write VTK file") pass diff --git a/examples/stresses.py b/examples/stresses.py index ca51412..9fbf92e 100644 --- a/examples/stresses.py +++ b/examples/stresses.py @@ -1,119 +1,122 @@ #!/usr/bin/env python3 # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import argparse import os import numpy as np import tamaas as tm from tamaas.dumpers import H5Dumper as Dumper from tamaas.dumpers._helper import hdf5toVTK parser = argparse.ArgumentParser( description="Hertzian tractios applied on elastic half-space") parser.add_argument("radius", type=float, help="Radius of sphere") parser.add_argument("load", type=float, help="Applied normal force") parser.add_argument("name", help="Output file name") parser.add_argument("--plots", help='Show surface normal stress', action="store_true") args = parser.parse_args() # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 127, 127] system_size = [0.25, 1., 1.] # Material contants E = 1. # Young's modulus nu = 0.3 # Poisson's ratio E_star = E/(1-nu**2) # Hertz modulus # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model.E = E model.nu = nu # Setup for integral operators residual = tm.ModelFactory.createResidual(model, 0, 0) # Coordinates x = np.linspace(0, system_size[1], discretization[1], endpoint=False) y = np.linspace(0, system_size[2], discretization[2], endpoint=False) x, y = np.meshgrid(x, y, indexing='ij') center = [0.5, 0.5] r = np.sqrt((x-center[0])**2 + (y-center[1])**2) -local_size = model.getBoundaryDiscretization() +# Span of local data +local_size = model.boundary_shape local_offset = tm.mpi.local_offset(r.shape) local_slice = np.s_[local_offset:local_offset+local_size[0], :] # Sphere R = args.radius P = args.load # Contact area a = (3*P*R/(4*E_star))**(1/3) p_0 = 3 * P / (2 * np.pi * a**2) # Pressure definition -traction = model.getTraction() +traction = model.traction hertz_traction = np.zeros(discretization[1:]) hertz_traction[r < a] = p_0 * np.sqrt(1 - (r[r < a] / a)**2) traction[..., 2] = hertz_traction[local_slice] # Array references -displacement = model.getDisplacement() -stress = residual.getStress() -gradient = residual.getVector() - -# Applying operator -boussinesq = model.getIntegralOperator("boussinesq") -boussinesq_gradient = model.getIntegralOperator("boussinesq_gradient") -boussinesq.apply(traction, displacement) -boussinesq_gradient.apply(traction, gradient) -model.applyElasticity(stress, gradient) +displacement = model.displacement +stress = model['stress'] +gradient = np.zeros_like(stress) + +# Getting integral operators +boussinesq = model.operators["boussinesq"] +boussinesq_gradient = model.operators["boussinesq_gradient"] + +# Applying operators +boussinesq(traction, displacement) +boussinesq_gradient(traction, gradient) +model.operators["hooke"](gradient, stress) # Dumper dumper_helper = Dumper(args.name, 'stress') model.addDumper(dumper_helper) model.dump() # Converting HDF dump to VTK with tm.mpi.sequential(): if tm.mpi.rank() == 0: hdf5toVTK(os.path.join('hdf5', 'stress_0000.h5'), args.name) if args.plots: import matplotlib.pyplot as plt # noqa fig, ax = plt.subplots(1, 2) fig.suptitle("Rank {}".format(tm.mpi.rank())) ax[0].set_title("Traction") ax[1].set_title("Normal Stress") ax[0].imshow(traction[..., 2]) ax[1].imshow(stress[0, ..., 2]) fig.tight_layout() plt.show() diff --git a/python/SConscript b/python/SConscript index 6859b58..088b8b9 100644 --- a/python/SConscript +++ b/python/SConscript @@ -1,162 +1,157 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @file # @section LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . from __future__ import print_function from os.path import abspath, join from SCons.Script import Import, Split, Copy, Dir Import('main_env') # Pybind11 wrapper env_pybind = main_env.Clone(SHLIBPREFIX='') # Remove pedantic warnings cxx_flags = env_pybind['CXXFLAGS'] try: del cxx_flags[cxx_flags.index('-pedantic')] except ValueError: pass env_pybind.Tool(pybind11) pybind_sources = Split(""" tamaas_module.cpp wrap/core.cpp wrap/percolation.cpp wrap/surface.cpp wrap/model.cpp wrap/solvers.cpp wrap/compute.cpp wrap/mpi.cpp wrap/test_features.cpp """) -# Adding legacy wrap code -if env_pybind['legacy_bem']: - env_pybind.AppendUnique(CPPDEFINES=['TAMAAS_LEGACY_BEM']) - pybind_sources += ["wrap/bem.cpp"] - # Setting paths to find libTamaas env_pybind.AppendUnique(LIBPATH=[Dir(join('#${build_dir}', 'src'))]) # Link against a static libTamaas if env_pybind['build_static_lib']: env_pybind.PrependUnique(LIBS=['Tamaas']) # keep other libs for link env_pybind['RPATH'] = "" # no need for rpath w/ static lib # Link against a dynamic libTamaas else: env_pybind.AppendUnique(RPATH=[ - "'$$$$ORIGIN/../../src'", #< path to lib in build_dir - "'$$$$ORIGIN/../../..'", #< path to lib in install prefix + "'$$$$ORIGIN/../../src'", # path to lib in build_dir + "'$$$$ORIGIN/../../..'", # path to lib in install prefix ]) env_pybind['LIBS'] = ['Tamaas'] # discard other libs for link # Building the pybind library tamaas_wrap = env_pybind.Pybind11Module( target='tamaas/_tamaas', source=pybind_sources, ) # For some reason link happens too early Import('libTamaas') env_pybind.Depends(tamaas_wrap, libTamaas) # Copying the __init__.py file with extra python classes copy_env = env_pybind.Clone() # Copying additional python files python_files = """ compute.py dumpers/__init__.py dumpers/_helper.py nonlinear_solvers/__init__.py """.split() targets = [tamaas_wrap] targets += [ copy_env.Command(join('tamaas', f), join(abspath(str(Dir('#python/tamaas'))), f), Copy("$TARGET", "$SOURCE")) for f in python_files ] targets.append(copy_env.Command('MANIFEST.in', '#python/MANIFEST.in', Copy("$TARGET", "$SOURCE"))) subst_env = env_pybind.Clone( SUBST_DICT={ '@version@': '$version', '@authors@': str(copy_env['authors']), '@email@': '$email', # TODO change when issue with unicode fixed # '@copyright@': '$copyright', # '@maintainer@': '$maintainer', } ) subst_env.Tool('textfile') targets.append(subst_env.Substfile('setup.py.in')) targets.append(subst_env.Substfile('tamaas/__init__.py.in')) # Defining alias for python builds main_env.Alias('build-python', targets) # Checking if we can use pip to install (more convenient for end-user) conf = Configure(main_env, custom_tests={'CheckPythonModule': CheckPythonModule}) has_pip = conf.CheckPythonModule('pip') install_env = conf.Finish() # Current build directory install_env['PYDIR'] = Dir('.') # Setting command line for installation if has_pip: install_env['PYINSTALLCOM'] = '${py_exec} -m pip install -U $PYOPTIONS .' install_env['PYDEVELOPCOM'] = \ '${py_exec} -m pip install $PYOPTIONS -e .[all]' else: install_env['PYINSTALLCOM'] = '${py_exec} setup.py install $PYOPTIONS' install_env['PYDEVELOPCOM'] = '${py_exec} setup.py develop $PYOPTIONS' install_env['py_version'] = get_python_version(install_env) install_env.PrependENVPath( 'PYTHONPATH', install_env.subst('${prefix}/lib/python${py_version}/site-packages')) # Specify install target python_install = install_env.Command( join('$prefix', 'dummy_target'), targets, install_env['PYINSTALLCOM'], PYOPTIONS='--prefix ${prefix}', chdir=install_env['PYDIR']) python_install_dev = install_env.Command( join('$prefix', 'dummy_target_local'), targets, install_env['PYDEVELOPCOM'], PYOPTIONS='--user', chdir=install_env['PYDIR']) main_env.Alias('install-python', python_install) main_env.Alias('dev', python_install_dev) diff --git a/python/cast.hh b/python/cast.hh index b97750c..b861add 100644 --- a/python/cast.hh +++ b/python/cast.hh @@ -1,211 +1,182 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ -#ifndef __CAST_HH__ -#define __CAST_HH__ +#ifndef CAST_HH +#define CAST_HH /* -------------------------------------------------------------------------- */ #include "grid_base.hh" +#include "grid.hh" #include "numpy.hh" -#include "surface.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace pybind11 { // Format descriptor necessary for correct wrap of tamaas complex type template struct format_descriptor< tamaas::complex, detail::enable_if_t::value>> { static constexpr const char c = format_descriptor::c; static constexpr const char value[3] = {'Z', c, '\0'}; static std::string format() { return std::string(value); } }; #ifndef PYBIND11_CPP17 template constexpr const char format_descriptor< tamaas::complex, detail::enable_if_t::value>>::value[3]; #endif namespace detail { // declare tamaas complex as a complex type for pybind11 template struct is_complex> : std::true_type {}; template struct is_fmt_numeric, detail::enable_if_t::value>> : std::true_type { static constexpr int index = is_fmt_numeric::index + 3; }; static inline handle policy_switch(return_value_policy policy, handle parent) { switch (policy) { case return_value_policy::copy: case return_value_policy::move: return handle(); case return_value_policy::automatic_reference: // happens in python-derived // classes case return_value_policy::reference: return none(); case return_value_policy::reference_internal: return parent; default: TAMAAS_EXCEPTION("Policy is not handled"); } } template handle grid_to_python(const tamaas::Grid& grid, return_value_policy policy, handle parent) { parent = policy_switch(policy, parent); // reusing variable std::vector sizes(dim); std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin()); if (grid.getNbComponents() != 1) sizes.push_back(grid.getNbComponents()); return array(sizes, grid.getInternalData(), parent).release(); } /** * Type caster for grid classes * inspired by https://tinyurl.com/y8m47qh3 from T. De Geus * and pybind11/eigen.h */ template