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