diff --git a/SConstruct b/SConstruct index 7d7e21c..a013aa2 100644 --- a/SConstruct +++ b/SConstruct @@ -1,430 +1,430 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 import SCons from os.path import join, abspath # 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.1.1", 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@epfl.ch', - copyright=u"Copyright (©) 2016-20 EPFL " + copyright=u"Copyright (©) 2016-2020 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""" FindFFTW(env, ['omp'], 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', '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 # Compiler detection compiler_default = os.getenv('CXX', 'g++') # 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=('omp', 'tbb'), # allowed_values=('omp', 'cuda'), 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', compiler_default), ('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), # Strip binary of extra info BoolVariable('strip_info', 'Strip binary of added information', 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 """ # 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['PRINT_CMD_LINE_FUNC'] = pretty_cmd_print # Include paths main_env.AppendUnique(CPPPATH=['#/src', '#/src/core', '#/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']) main_env['CXXFLAGS'] = Split(main_env['CXXFLAGS']) main_env.AppendUnique(CXXFLAGS='-std=c++14 -Wall -Wextra -pedantic'.split()) # Append openmp flags regardless of backend for fftw_omp main_env.AppendUnique(CXXFLAGS=omp_flags[cxx]) main_env.AppendUnique(LINKFLAGS=omp_flags[cxx]) # 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 # Flags and options if main_env['build_python']: main_env.AppendUnique(CPPDEFINES=['USE_PYTHON']) if main_env['legacy_bem']: main_env.AppendUnique(CPPDEFINES=['LEGACY_BEM']) if main_env['build_type'] == 'debug': main_env.AppendUnique(CPPDEFINES=['TAMAAS_DEBUG']) # Define the scalar types main_env.AppendUnique(CPPDEFINES={'REAL_TYPE': '${real_type}', 'INT_TYPE': '${integer_type}'}) if main_env['real_type'] == 'long double': main_env.AppendUnique(CPPDEFINES=['LONG_PRECISION']) # Compilation flags cxxflags_dict = { "debug": Split("-g -O0"), "profiling": Split("-g -O3 -fno-omit-frame-pointer"), "release": Split("-O3") } # Link flags for shared libs shlinkflags_dict = { "debug": [], "profiling": Split("-g -O3 -fno-omit-frame-pointer"), "release": [] } 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}') shlinkflags_dict[build_type].append('-fsanitize=${sanitizer}') main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type]) main_env.AppendUnique(SHLINKFLAGS=shlinkflags_dict[build_type]) main_env.AppendUnique(LINKFLAGS=shlinkflags_dict[build_type]) # Adding main_env['LIBPATH'] = [abspath(join(build_dir, 'src')), abspath(join(main_env['prefix'], 'lib'))] main_env['RPATH'] = "$LIBPATH" if main_env['should_configure']: detect_dependencies(main_env) # Writing information file main_env.Tool('textfile') main_env['SUBST_DICT'] = get_git_subst() main_env['SUBST_DICT'].update({ '@build_type@': '$build_type', '@build_dir@': abspath(build_dir), }) if main_env['strip_info']: for k in main_env['SUBST_DICT']: main_env['SUBST_DICT'][k] = "" # 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 = main_env.Command('#.phony_doc', '', 'make -C {}'.format(Dir('#/doc'))) 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') diff --git a/examples/adhesion.py b/examples/adhesion.py index 07cbea9..f2ea9cf 100644 --- a/examples/adhesion.py +++ b/examples/adhesion.py @@ -1,123 +1,123 @@ #!/usr/bin/env python3 # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 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() # Initialize threads and fftw tm.initialize() # 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) # 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.figure() zones = np.zeros_like(tractions) tol = 1e-6 zones[tractions > tol] = 1 zones[tractions < -tol] = -1 plt.imshow(zones, cmap='Greys') # Cleanup threads tm.finalize() plt.show() diff --git a/examples/legacy/adhesion.py b/examples/legacy/adhesion.py index 2f834d0..8bc9036 100644 --- a/examples/legacy/adhesion.py +++ b/examples/legacy/adhesion.py @@ -1,161 +1,161 @@ #!/usr/bin/env python # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 51ec9bb..a03370b 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-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 fe0a200..25abd57 100644 --- a/examples/legacy/contact.py +++ b/examples/legacy/contact.py @@ -1,73 +1,73 @@ #!/usr/bin/env python # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 4cb6444..b68fb4c 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-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 bf25b12..cb5bfaa 100644 --- a/examples/legacy/hertz.py +++ b/examples/legacy/hertz.py @@ -1,116 +1,116 @@ #!/usr/bin/env python # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 1d3ea27..8cdd108 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-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 ec236fe..ff57702 100644 --- a/examples/plasticity.py +++ b/examples/plasticity.py @@ -1,74 +1,74 @@ #!/usr/bin/env python3 # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 UVWDumper as Dumper from tamaas.nonlinear_solvers import DFSANESolver as Solver tm.initialize(2) # 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) epsolver = Solver(residual, model) # Setup for contact x = np.linspace(0, system_size[1], discretization[1], endpoint=False) y = np.linspace(0, system_size[2], discretization[2], endpoint=False) 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) csolver = tm.PolonskyKeerRey(model, 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() print("---> Solved load step {}/{}".format(i+1, len(loads))) diff --git a/examples/rough_contact.py b/examples/rough_contact.py index abfb80e..6afc382 100644 --- a/examples/rough_contact.py +++ b/examples/rough_contact.py @@ -1,74 +1,74 @@ #!/usr/bin/env python3 # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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.initialize() 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()) print(model.getTraction().mean()) # Cleanup threads tm.finalize() plt.show() diff --git a/examples/saturated.py b/examples/saturated.py index fb2416e..c31dcc8 100644 --- a/examples/saturated.py +++ b/examples/saturated.py @@ -1,67 +1,67 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 tm.initialize() 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.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.colorbar() plt.show() tm.finalize() diff --git a/examples/statistics.py b/examples/statistics.py index 528c00a..c6d296d 100644 --- a/examples/statistics.py +++ b/examples/statistics.py @@ -1,87 +1,87 @@ #!/usr/bin/env python3 # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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.initialize() 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 4049bc0..0a92298 100644 --- a/examples/stresses.py +++ b/examples/stresses.py @@ -1,96 +1,96 @@ #!/usr/bin/env python3 # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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 import argparse from tamaas.dumpers import UVWDumper 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") args = parser.parse_args() tm.initialize() # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [127, 127, 127] system_size = [1., 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) # 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[r < a, 2] = p_0 * np.sqrt(1 - (r[r < a] / a)**2) # 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 = UVWDumper(args.name, 'stress') model.addDumper(dumper_helper) model.dump() print("Done") tm.finalize() diff --git a/python/SConscript b/python/SConscript index 57d1410..85436ce 100644 --- a/python/SConscript +++ b/python/SConscript @@ -1,143 +1,143 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @file # @section LICENSE # -# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), +# Copyright (©) 2016-2020 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/test_features.cpp """) # Adding legacy wrap code if main_env['legacy_bem']: pybind_sources += ["wrap/bem.cpp"] # Building the pybind library tamaas_wrap = env_pybind.Pybind11Module( target='tamaas/_tamaas', source=pybind_sources, LIBS=['Tamaas'], ) # 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 $PYOPTIONS .' install_env['PYDEVELOPCOM'] = '${py_exec} -m pip install $PYOPTIONS -e .' 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 5001747..c672217 100644 --- a/python/cast.hh +++ b/python/cast.hh @@ -1,211 +1,211 @@ /** * @file * @section LICENSE * - * Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), + * Copyright (©) 2016-2020 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