diff --git a/.gitignore b/.gitignore index 9460480..7f24725 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,6 @@ *~ .sconf_temp .sconsign.dblite *.log build* -site_scons/site_tools/swig.pyc +*.pyc diff --git a/SConstruct b/SConstruct index bf278af..d9ad04c 100644 --- a/SConstruct +++ b/SConstruct @@ -1,107 +1,108 @@ import os,sys #path = '/usr/lib/scons/' #sys.path.append(path) #from SCons import * #from SCons.Environment import * #from SCons.Variables import * #from SCons.Util import * #script_dir = GetOption('file') #if script_dir == []: # script_dir = '.' #else: # script_dir = script_dir[0] colors = {} colors['cyan'] = '' colors['purple'] = '' colors['blue'] = '' colors['green'] = '' colors['yellow'] = '' colors['red'] = '' colors['end'] = '' #colors['cyan'] = '\033[96m' #colors['purple'] = '\033[95m' #colors['blue'] = '\033[94m' #colors['green'] = '\033[92m' #colors['yellow'] = '\033[93m' #colors['red'] = '\033[91m' #colors['end'] = '\033[0m' main_env = Environment(CXX= 'g++', ENV = os.environ) vars = Variables('build-setup.conf') vars.Add(EnumVariable('build_type', 'Build type', 'release', allowed_values=('release', 'profiling', 'debug'), ignorecase=2)) vars.Add('prefix', 'Prefix where to install', '/usr/local') vars.Add(BoolVariable('timer', 'Activate the timer possibilities', False)) vars.Add(BoolVariable('verbose', 'Activate verbosity', False)) vars.Update(main_env) Help(vars.GenerateHelpText(main_env)) build_type = main_env['build_type'] build_dir = 'build-' + main_env['build_type'] print "Building in ", build_dir timer_flag = main_env['timer'] verbose = main_env['verbose'] if not verbose: main_env['SHCXXCOMSTR'] = unicode('{0}[Compiling] {1}$SOURCE{2}'.format(colors['green'],colors['blue'],colors['end'])) main_env['SHLINKCOMSTR'] = unicode('{0}[Linking] {1}$TARGET{2}'.format(colors['purple'],colors['blue'],colors['end'])) main_env['SWIGCOMSTR'] = unicode('{0}[Swig] {1}$SOURCE{2}'.format(colors['yellow'],colors['blue'],colors['end'])) main_env.AppendUnique(CPPPATH=['#/src', '#/python']) main_env.AppendUnique(LIBS='gomp') # Flags and options main_env.AppendUnique(CXXFLAGS=['-std=c++11', '-Wall', '-fopenmp']) main_env.AppendUnique(SHLINKFLAGS=['-fopenmp']) if timer_flag: main_env.AppendUnique(CPPDEFINES=['USING_TIMER']) cxxflags_dict = { "debug":Split("-g -O0"), "profiling":Split("-g -pg -O2"), "release":Split("-O3") } shlinkflags_dict = { "debug":[], "profiling":['-pg'], "release":[] } main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type]) main_env.AppendUnique(SHLINKFLAGS=shlinkflags_dict[build_type]) main_env['LIBPATH'] = [os.path.abspath(os.path.join(build_dir,'src'))] main_env['RPATH'] = "$LIBPATH" # save new values vars.Save('build-setup.conf', main_env) Export('main_env') SConscript('src/SConscript',variant_dir=os.path.join(build_dir,'src'),duplicate=False) SConscript('python/SConscript',variant_dir=os.path.join(build_dir,'python'),duplicate=False) +SConscript('tests/SConscript',variant_dir=os.path.join(build_dir,'tests'),duplicate=False) # saving the env file try: env_file = open(os.path.join(build_dir,'tamaas_environement.sh'),'w') env_file.write(""" export PYTHONPATH=$PYTHONPATH:{0}/python export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:{0}/src """.format(os.path.abspath(build_dir))) env_file.close() except: pass diff --git a/python-examples/hertz.py b/python-examples/hertz.py index b84716b..ef3c264 100644 --- a/python-examples/hertz.py +++ b/python-examples/hertz.py @@ -1,146 +1,123 @@ #!/usr/bin/python # -*- coding: utf-8 -*- ##* # # @author Guillaume Anciaux # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas 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 Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. 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("--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()) -# std::stringstream fname -# fname << n << "-" << r << "-" -# << pressure -# bem.getTractions().dumpToParaview(fname.str()) -# fname << "-disp" -# bem.getTrueDisplacements().dumpToParaview(fname.str()) -# -# std::cout << "* " << i << "/" << Ninc << ", pressure = " << pressure << " area " << A << std::endl dumpHertzPrediction(file,bem,pressure,r,Estar) -# if (plot_surface) { -# std::stringstream prefix -# prefix << "-A-" << A -# << "-pressure-" << pressure -# << "-step-" << i -# -# std::string filename = parser.makeFilename("pressures" + prefix.str(),"plot") -# bem.getTractions().dumpToTextFile(filename) -# std::string filename_disp = parser.makeFilename("displacement"+ prefix.str(),"plot") -# bem.getDisplacements().dumpToTextFile(filename_disp) -# std::string filename_gap = parser.makeFilename("gap"+ prefix.str(),"plot") -# bem.getGap().dumpToTextFile(filename_gap) -# } - - if A == 1.0: file.close() break tamaas.dumpTimes() ################################################################ import sys import argparse import numpy as np main(sys.argv) diff --git a/src/SConscript b/src/SConscript index 004baf7..f8589c7 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,79 +1,80 @@ Import('main_env') import distutils.sysconfig import os fftw_include = "" fftw_library = "" if 'FFTW_INCLUDE' in main_env['ENV']: fftw_include = main_env['ENV']['FFTW_INCLUDE'] if 'FFTW_LIBRARY' in main_env['ENV']: fftw_library = main_env['ENV']['FFTW_LIBRARY'] env = main_env.Clone( tools = [fftw], FFTW_LIBRARY_WISH = ['main', 'omp'], FFTW_INCLUDE_DIR = fftw_include, FFTW_LIBRARY_DIR = fftw_library, ) # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_voss.cpp surface_generator_crenel.cpp surface_generator_ellipsoid.cpp surface_generator_filter.cpp surface_generator_filter_bessel.cpp surface_generator_filter_fft.cpp """.split() #env.SharedLibrary('Generator',generator_list) # Lib SURFACE surface_list = """ map_2d.cpp map_2d_square.cpp surface.cpp surface_timer.cpp tamaas.cpp """.split() #env.SharedLibrary('Surface',surface_list) # Lib PERCOLATION percolation_list = """ cluster_grow.cpp contact_area.cpp contact_cluster.cpp contact_cluster_collection.cpp """.split() #env.SharedLibrary('Percolation',percolation_list) # BEM PERCOLATION bem_list = """ bem_kato.cpp bem_polonski.cpp bem_gigi.cpp bem_gigipol.cpp bem_penalty.cpp bem_uzawa.cpp bem_fft_base.cpp functional.cpp elastic_energy_functional.cpp adhesion_functional.cpp fftransform.cpp +fft_plan_manager.cpp """.split() #env.SharedLibrary('BEM',bem_list) rough_contact_list = generator_list + surface_list + percolation_list + bem_list current_files = set(os.listdir(Dir('.').srcnode().abspath)) unexpected_files = list(current_files - set(rough_contact_list + ['SConscript'])) unexpected_headers = [ f for f in unexpected_files if os.path.splitext(f)[1] == '.hh' ] unexpected_files = [ f for f in unexpected_files if not os.path.splitext(f)[1] == '.hh' ] #print "Unexpected files : {}".format(unexpected_files) #print "Include path : {}".format([str(x) for x in env['CPPPATH']]) #print unexpected_headers env.SharedLibrary('Tamaas',rough_contact_list) diff --git a/src/bem_fft_base.cpp b/src/bem_fft_base.cpp index 3193179..8638fe3 100644 --- a/src/bem_fft_base.cpp +++ b/src/bem_fft_base.cpp @@ -1,305 +1,284 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" #include "surface_statistics.hh" #include "elastic_energy_functional.hh" +#include "fft_plan_manager.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ BemFFTBase::BemFFTBase(Surface & p): BemInterface(p), functional(new ElasticEnergyFunctional(*this)), surface_tractions(p.size(),p.getL()), surface_displacements(p.size(),p.getL()), - surface_spectral_tractions(p.size(),p.getL()), - surface_spectral_displacements(p.size(),p.getL()), - pressure_transform(p.size(), surface_tractions, surface_spectral_tractions), - displacement_transform(p.size(), surface_displacements, surface_spectral_displacements), + surface_spectral_input(p.size(),p.getL()), + surface_spectral_output(p.size(),p.getL()), surface_spectral_influence_disp(p.size(),p.getL()), true_displacements(p.size(),p.getL()), old_displacements(p.size(),p.getL()), gap(p.size(),p.getL()), e_star(std::nan("")), nthreads(1), tmp_coeff(1./M_PI){ // this->setNumberOfThreads(this->nthreads); max_iterations = 100000; dump_freq = 100; std::cerr << this << " is setting up the surfaces"; std::cerr << &p << " has size " << p.size() << std::endl; } /* -------------------------------------------------------------------------- */ BemFFTBase::~BemFFTBase() { delete functional; } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getTractions() const{ return surface_tractions; } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getDisplacements()const{ return surface_displacements; } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getSpectralInfluenceOverDisplacement(){ return surface_spectral_influence_disp; } /* -------------------------------------------------------------------------- */ -const SurfaceComplex & BemFFTBase::getSpectralDisplacements()const{ - return surface_spectral_displacements; -} -/* -------------------------------------------------------------------------- */ -const SurfaceComplex & BemFFTBase::getSpectralTractions()const{ - return surface_spectral_tractions; -} -/* -------------------------------------------------------------------------- */ void BemFFTBase::setNumberOfThreads(UInt nthreads){ this->nthreads = nthreads; omp_set_num_threads(nthreads); } /* -------------------------------------------------------------------------- */ void BemFFTBase::setEffectiveModulus(Real e_star){ this->e_star = e_star; this->computeSpectralInfluenceOverDisplacement(); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeSpectralInfluenceOverDisplacement(){ STARTTIMER("computeSpectralInfluenceOverDisplacement"); if (std::isnan(tmp_coeff)) SURFACE_FATAL("constant tmp_coeff is nan"); if (std::isnan(e_star)) SURFACE_FATAL("constant e_star is nan"); UInt n = surface.size(); surface_spectral_influence_disp(0) = 0.; // everything but the external lines and central lines #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { for (UInt j = 1; j < n/2; ++j) { Real d = sqrt(i*i+j*j); surface_spectral_influence_disp(i,j) = 1./d; surface_spectral_influence_disp(n-i,j) = 1./d; surface_spectral_influence_disp(i,n-j) = 1./d; surface_spectral_influence_disp(n-i,n-j) = 1./d; } } // external line (i or j == 0) #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { Real d = i; surface_spectral_influence_disp(i,0) = 1./d; surface_spectral_influence_disp(n-i,0) = 1./d; surface_spectral_influence_disp(0,i) = 1./d; surface_spectral_influence_disp(0,n-i) = 1./d; } //internal line (i or j == n/2) #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { Real d = sqrt(i*i+n/2*n/2); surface_spectral_influence_disp(i,n/2) = 1./d; surface_spectral_influence_disp(n-i,n/2) = 1./d; surface_spectral_influence_disp(n/2,i) = 1./d; surface_spectral_influence_disp(n/2,n-i) = 1./d; } { //i == n/2 j == n/2 Real d = sqrt(2*n/2*n/2); surface_spectral_influence_disp(n/2,n/2) = 1./d; //i == 0 j == n/2 d = n/2; surface_spectral_influence_disp(0,n/2) = 1./d; //i == n/2 j == 0 surface_spectral_influence_disp(n/2,0) = 1./d; } UInt size = n*n; Real L = surface.getL(); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { surface_spectral_influence_disp(i) *= tmp_coeff/this->e_star*L; } STOPTIMER("computeSpectralInfluenceOverDisplacement"); } /* -------------------------------------------------------------------------- */ -void BemFFTBase::computeSpectralDisplacements(){ - - STARTTIMER("computeSpectralDisplacement"); - - UInt n = surface.size(); - UInt size = n*n; - - // std::cerr << "000000000000 " << surface_spectral_tractions(n*n-1) << std::endl; - // std::cerr << "000000000000 " << surface_spectral_displacements(n*n-1) << std::endl; +void BemFFTBase::computeDisplacementsFromTractions(){ - surface_spectral_displacements = surface_spectral_tractions; - // std::cerr << "AAAAAAAAAAAA " << surface_spectral_tractions(n*n-1) << std::endl; - // std::cerr << "AAAAAAAAAAAA " << surface_spectral_displacements(n*n-1) << std::endl; - // std::cerr << "AAAAAAAAAAAA " << &surface_spectral_displacements(n*n-1) << std::endl; - // std::cerr << "AAAAAAAAAAAA " << &surface_spectral_displacements(0) << std::endl; + STARTTIMER("computeDisplacementsFromTractions"); -#pragma omp parallel for - for (UInt i = 0; i < size; ++i) { - surface_spectral_displacements(i) *= surface_spectral_influence_disp(i); - } - - // std::cerr << "BBBBBBBBBBBB " << surface_spectral_displacements(n*n-1) << std::endl; - // std::cerr << "BBBBBBBBBBBB " << surface_spectral_influence_disp(n*n-1) << std::endl; + this->applyInfluenceFunctions(surface_tractions,surface_displacements); - - STOPTIMER("computeSpectralDisplacement"); + STOPTIMER("computeDisplacementFromTractions"); } /* -------------------------------------------------------------------------- */ - -void BemFFTBase::computeSpectralTractions(){ - - STARTTIMER("computeSpectralTraction"); - - pressure_transform.forwardTransform(); - - STOPTIMER("computeSpectralTraction"); - +void BemFFTBase::computeTractionsFromDisplacements(){ + + STARTTIMER("computeDisplacementsFromTractions"); + + this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions); + + STOPTIMER("computeDisplacementFromTractions"); + } /* -------------------------------------------------------------------------- */ +void BemFFTBase::applyInfluenceFunctions(Surface & input, Surface & output){ -void BemFFTBase::computeDisplacements(){ + STARTTIMER("applyInfluenceFunctions"); - STARTTIMER("computeDisplacement"); + FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); + + plan1.forwardTransform(); - this->computeSpectralTractions(); - this->computeSpectralDisplacements(); - this->computeRealDisplacementsFromSpectral(); + surface_spectral_output *= surface_spectral_influence_disp; - STOPTIMER("computeDisplacement"); + FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); -} + plan2.backwardTransform(); -/* -------------------------------------------------------------------------- */ - -void BemFFTBase::computeRealTractionsFromSpectral(){ + STOPTIMER("applyInfluenceFunctions"); - STARTTIMER("computeRealTractionsFromSpectral"); +} - pressure_transform.backwardTransform(); +/* -------------------------------------------------------------------------- */ - STOPTIMER("computeRealDisplacementFromSpectral"); +void BemFFTBase::applyInverseInfluenceFunctions(Surface & input, Surface & output){ -} + STARTTIMER("applyInfluenceFunctions"); -/* -------------------------------------------------------------------------- */ + FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); + + plan1.forwardTransform(); + UInt n = surface_spectral_input.size(); + UInt size = n*n; -void BemFFTBase::computeRealDisplacementsFromSpectral(){ + surface_spectral_output(0) = 0; + +#pragma omp parallel for + for (UInt i = 1; i < size; ++i) { + surface_spectral_output(i) /= surface_spectral_influence_disp(i); + } - STARTTIMER("computeRealDisplacementFromSpectral"); + FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); - displacement_transform.backwardTransform(); + plan2.backwardTransform(); - STOPTIMER("computeRealDisplacementFromSpectral"); + STOPTIMER("applyInfluenceFunctions"); } /* -------------------------------------------------------------------------- */ - void BemFFTBase::computeGaps() { UInt size = surface.size(); #pragma omp parallel for for (UInt i = 0 ; i < size*size ; i++) { - gap(i) = true_displacements(i) - surface(i); + //gap(i) = true_displacements(i) - surface(i); + gap(i) = surface(i) - true_displacements(i); } } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeTrueDisplacements(){ this->true_displacements = this->surface_displacements; UInt n = surface.size(); UInt size = n*n; Real shift = 1e300; for (UInt i = 0; i < size; ++i) { if (surface_displacements(i) - shift - surface(i) < 0.){ shift = surface_displacements(i) - surface(i); } } for (UInt i = 0; i < size; ++i) { true_displacements(i) = surface_displacements(i) - shift; } } /* -------------------------------------------------------------------------- */ void BemFFTBase::setFunctional(Functional *new_funct) { delete this->functional; this->functional = new_funct; } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_fft_base.hh b/src/bem_fft_base.hh index fe111cf..fa2d381 100644 --- a/src/bem_fft_base.hh +++ b/src/bem_fft_base.hh @@ -1,167 +1,157 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_FFT_BASE_H #define BEM_FFT_BASE_H /* -------------------------------------------------------------------------- */ #include "bem_interface.hh" #include #include "surface_statistics.hh" #include "fftransform.hh" #include "functional.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemFFTBase : public BemInterface { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemFFTBase(Surface & p); virtual ~BemFFTBase(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: + //! compute the displacement from traction input + virtual void computeDisplacementsFromTractions(); + //! compute the tractions from displacement input + virtual void computeTractionsFromDisplacements(); + //! apply influence funtions + virtual void applyInfluenceFunctions(Surface & input, Surface & output); + //! apply inverse influence funtions + virtual void applyInverseInfluenceFunctions(Surface & input, Surface & output); //! compute the influcence coefficients of pressure over displacement virtual void computeSpectralInfluenceOverDisplacement(); - //! compute real tractions from fourier represnetation (IFFT) - virtual void computeRealTractionsFromSpectral(); - //! compute real displacements from fourier represnetation (IFFT) - virtual void computeRealDisplacementsFromSpectral(); - //! compute the displacement in real space - virtual void computeDisplacements(); - //! compute the displacement fourier space - virtual void computeSpectralDisplacements(); - //! compute the tractions in fourier space - virtual void computeSpectralTractions(); - //! compute the gradient the objective function (with respect to traction) + //! compute the displacements virtual void computeTrueDisplacements(); //! compute the gap virtual void computeGaps(); //! get the number of iterations virtual UInt getNbIterations() const = 0; //! get the convergence steps of iterations virtual const std::vector & getConvergenceIterations() const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! return the influcence coefficients virtual const Surface & getSpectralInfluenceOverDisplacement(); //! return the actual tractions virtual const Surface & getTractions() const; //! return the actual displacements virtual const Surface & getDisplacements()const; - //! return the actual displacements - virtual const SurfaceComplex & getSpectralDisplacements()const; - //! return the actual displacements - virtual const SurfaceComplex & getSpectralTractions()const; //! get the computed true displacement const Surface & getTrueDisplacements(){return true_displacements;}; //! get the computed gap const Surface & getGap(){return gap;}; //! set the number of threads to be used to perform the calculations void setNumberOfThreads(UInt nthreads); //! get the number of threads UInt getNumberOfThreads() const { return nthreads; } //! set the elastic effective modulus void setEffectiveModulus(Real e_star); //! set the maximal number of iterations void setMaxIterations(UInt max_iter){this->max_iterations = max_iter;}; //! set dump freq void setDumpFreq(Real dump_freq){this->dump_freq = dump_freq;}; //! set functional void setFunctional(::tamaas::Functional * new_funct); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! Functional for optimization Functional * functional; //! surface tractions Surface surface_tractions; //! surface displacement Surface surface_displacements; - //! surface tractions - SurfaceComplex surface_spectral_tractions; - //! surface displacement - SurfaceComplex surface_spectral_displacements; - //! FF Transform for pressures - FFTransform pressure_transform; - //! FF Transform for displacements - FFTransform displacement_transform; - //! surface influcence coefficient over displacement + //! surface input in spectral + SurfaceComplex surface_spectral_input; + //! surface output in spectral + SurfaceComplex surface_spectral_output; + //! surface influence coefficient over displacement Surface surface_spectral_influence_disp; //! shifted displacement in order to get the true displacement Surface true_displacements; //! old displacements (for stopping criterion) Surface old_displacements; //! computing of the gap Surface gap; //! Effective young modulus Real e_star; //! Number of threads to use during calculations UInt nthreads; //! maximum of iterations UInt max_iterations; //! number of iterations UInt nb_iterations; //! number of iterations std::vector convergence_iterations; //! frequency at which to dump iteration info UInt dump_freq; /* ------------------------------------------------------------------------ */ /* temporary trick about some obscure scaling factor */ /* ------------------------------------------------------------------------ */ /* 0 2 = 2*1 1 4 = 2*2 2 10 = 2*5 3 20 = 2*10 4 34 = 2*12 */ //#define tmp_coeff (2.*n/M_PI) //#define tmp_coeff (n/M_PI) const Real tmp_coeff; }; __END_TAMAAS__ #endif /* BEM_FFT_BASE_H */ diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp index b29fc5a..b62f5d5 100644 --- a/src/bem_gigi.cpp +++ b/src/bem_gigi.cpp @@ -1,356 +1,339 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_gigi.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemGigi::computeEquilibrium(Real epsilon, Real mean_displacement) { this->computeSpectralInfluenceOverDisplacement(); Real Rold = 1.; Real f = 1e300; this->search_direction = 0.; this->true_displacements =0; this->true_displacements = mean_displacement; this->computeGaps(); this->optimizeToMeanDisplacement(mean_displacement); this->computeGaps(); convergence_iterations.clear(); nb_iterations = 0; max_iterations =50; Real R=0; while (f> epsilon && nb_iterations++ < max_iterations) { this->computeGaps(); this->functional->computeGradF(); this->updateSearchDirection(R / Rold); Rold = R; Real alpha = this->computeOptimalStep(); std::cout << "alpha vaut "<< alpha<< std::endl; this->old_displacements=this->true_displacements; this->updateDisplacements(alpha); this->computeGaps(); this->optimizeToMeanDisplacement(mean_displacement); //espace admissible this->computeGaps(); f=computeStoppingCriterion(); } this->computePressures(mean_displacement); return f; } /* -------------------------------------------------------------------------- */ Real BemGigi::computeStoppingCriterion() { Real crit=0.; Real disp_norm = 0.; UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for reduction(+:crit, disp_norm) for (UInt i = 0; i < size; ++i) { crit += (true_displacements(i)-old_displacements(i))*(true_displacements(i)-old_displacements(i)); disp_norm += (true_displacements(i)*true_displacements(i)); } return fabs(crit) / disp_norm; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void BemGigi::optimizeToMeanDisplacement(Real imposed_mean) { Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0); UInt n = surface.size(); UInt size = n * n; // Initial guesses for upper and lower bound Real step_min = -10; Real step_max = 10; // Gaps for upper and lower bound Real gap_min = positiveGapAverage(step_min); Real gap_max = positiveGapAverage(step_max); UInt max_expansion = 8; // Expand bounds if necessary for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10) gap_max = positiveGapAverage(step_max); for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10) gap_min = positiveGapAverage(step_min); Real g = 0.; Real epsilon = 1e-12; Real step = 0; while (fabs(g - target_value) > epsilon) { step = (step_min + step_max) / 2.; g = positiveGapAverage(step); if (g > target_value) step_max = step; else if (g < target_value) step_min = step; else { step_max = step; step_min = step; } } step = (step_min + step_max) / 2.; #pragma omp parallel for for (UInt i = 0 ; i < size ; i++) { gap(i) += step; if (gap(i) < 0) gap(i) = 0; true_displacements(i) = gap(i) + surface(i); } } /* -------------------------------------------------------------------------- */ Real BemGigi::positiveGapAverage(Real shift) { UInt n = surface.size(); Real res = 0; #pragma omp parallel for reduction(+: res) for (UInt i = 0 ; i < n * n ; i++) { Real shifted_gap = gap(i) + shift; res += shifted_gap * (shifted_gap > 0); } return res / (n * n); } /* -------------------------------------------------------------------------- */ void BemGigi::updateSearchDirection(Real factor) { STARTTIMER("updateSearchDirection"); UInt n = surface.size(); UInt size = n*n; Surface & gradF = this->functional->getGradF(); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->search_direction(i)= gradF(i); } STOPTIMER("updateSearchDirection"); } /* -------------------------------------------------------------------------- */ Real BemGigi::computeOptimalStep() { STARTTIMER("computeOptimalStep"); search_direction.FFTTransform(spectral_search_direction, nthreads); UInt n = surface.size(); UInt size = n*n; spectral_search_direction(0) = 0; #pragma omp parallel for for (UInt i = 1; i < size; ++i) { spectral_search_direction(i) /= this->surface_spectral_influence_disp(i); } // /!\ does not contain the spectral search direction anymore spectral_search_direction.FFTITransform(nthreads); // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0); // spectral_search_direction -= average; Surface & gradF = this->functional->getGradF(); Real numerator = 0., denominator = 0.; #pragma omp parallel for reduction(+: numerator, denominator) for (UInt i = 0; i < size; ++i) { numerator += gradF(i) * search_direction(i); denominator += spectral_search_direction(i).real() * search_direction(i); } Real alpha = numerator / denominator; STOPTIMER("computeOptimalStep"); return alpha; } /* -------------------------------------------------------------------------- */ Real BemGigi::computeTau() { STARTTIMER("computeOptimalStep"); search_direction.FFTTransform(spectral_search_direction, nthreads); UInt n = surface.size(); UInt size = n*n; spectral_search_direction(0) = 0; #pragma omp parallel for for (UInt i = 1; i < size; ++i) { spectral_search_direction(i) /= this->surface_spectral_influence_disp(i); } // /!\ does not contain the spectral search direction anymore spectral_search_direction.FFTITransform(nthreads); // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0); // spectral_search_direction -= average; Surface & gradF = this->functional->getGradF(); Real numerator = 0., denominator = 0.; #pragma omp parallel for reduction(+: numerator, denominator) for (UInt i = 0; i < size; ++i) { if (this->gap(i) > 0) { numerator += gradF(i) * search_direction(i); denominator += spectral_search_direction(i).real() * search_direction(i);} } Real tau = numerator / denominator; STOPTIMER("computeOptimalStep"); return tau; } /* -------------------------------------------------------------------------- */ void BemGigi::updateDisplacements(Real alpha) { STARTTIMER("updateDisplacements"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->true_displacements(i) -= alpha*this->search_direction(i); } STOPTIMER("updateDisplacements"); } /* -------------------------------------------------------------------------- */ void BemGigi::emptyOverlap() { STARTTIMER("emptyoverlap"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (gap(i) < 0) this->true_displacements(i) = this->surface(i); } STOPTIMER("emptyoverlap"); } /* -------------------------------------------------------------------------- */ void BemGigi::enforceMeanDisplacement(Real mean_displacement) { STARTTIMER("enforceMeanDisplacement"); UInt n = surface.size(); UInt size = n*n; Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0); Real factor = mean_displacement / average; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->true_displacements(i) *= factor; } STOPTIMER("enforceMeanDisplacement"); } /* -------------------------------------------------------------------------- */ void BemGigi::computePressures(Real mean_displacement) { - UInt n = surface.size(); - UInt size = n*n; - - - this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads); - this->surface_spectral_tractions(0) = 0; - -#pragma omp parallel for - for (UInt i = 1 ; i < size ; ++i) { - this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i); - } - - this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads); - - - - - + this->computeTractionsFromDisplacements(); } __END_TAMAAS__ diff --git a/src/bem_gigipol.cpp b/src/bem_gigipol.cpp index 2ac91cf..3456d7a 100644 --- a/src/bem_gigipol.cpp +++ b/src/bem_gigipol.cpp @@ -1,415 +1,408 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_gigipol.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemGigipol::computeEquilibrium(Real epsilon, Real mean_displacement) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; this->surface_spectral_r = 0.; this->search_direction = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; UInt n = surface.size(); UInt size = n*n; this->true_displacements =0; this->true_displacements = mean_displacement; this->computeGaps(); this->optimizeToMeanDisplacement(mean_displacement); convergence_iterations.clear(); nb_iterations = 0;max_iterations=500; while(f > epsilon && nb_iterations < max_iterations) { fPrevious = f; this->functional->computeGradF(); this->search_direction= this->functional->getGradF(); Real gbar = this->computeMeanPressuresInNonContact(); this->search_direction -= gbar; Real G = this->computeG(); std::cout << "G vaut "<< G<< std::endl; this->updateT(G,Gold,delta); Real tau = this->computeTau(); this->old_displacements=this->true_displacements; Gold = G; #pragma omp parallel for for (unsigned int i = 0; i < size; ++i) { this->true_displacements(i)-=tau*this->surface_t(i); } //Projection on admissible space this->computeGaps(); #pragma omp parallel for for (unsigned int i = 0; i < size; ++i) { if (this->gap(i)<0){ this->true_displacements(i)=this->surface(i);} } this->computeGaps(); delta = this->updateDisplacements(tau, mean_displacement); // std::cout << "delta vaut "<< delta << std::endl; //Projection on admissible space this->computeGaps(); this->enforceMeanDisplacement(mean_displacement); this->computeGaps(); this->computePressures(mean_displacement); f=this->computeStoppingCriterion(); std::cout << "f vaut "<< f << std::endl; if (nb_iterations % dump_freq ==1000){ Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions); std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fPrevious << " " << A << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computePressures(mean_displacement); return f; } /* -------------------------------------------------------------------------- */ Real BemGigipol::computeStoppingCriterion() { Real crit=0.; UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for reduction(+:crit) for (UInt i = 0; i < size; ++i) { crit += std::abs(this->gap(i)*(this->surface_tractions(i))); ; } return crit ; } /* -------------------------------------------------------------------------- */ Real BemGigipol::computeG(){ STARTTIMER("computeG"); unsigned int n = surface.size(); unsigned int size = n*n; Real res = 0.; #pragma omp parallel for reduction(+: res) for (unsigned int i = 0; i < size; ++i) { if (this->gap(i) > 0.) { Real val = this->search_direction(i); res += val*val;} } STOPTIMER("computeG"); return res; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ Real BemGigipol::computeMeanPressuresInNonContact(){ STARTTIMER("computeMeanPressuresInNonContact"); unsigned int n = surface.size(); unsigned int size = n*n; Real res = 0.; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,res) for (unsigned int i = 0; i < size; ++i) { if (this->gap(i) == 0. ) continue; ++nb_contact; res += this->search_direction(i); } res /= nb_contact; STOPTIMER("computeMeanPressuresInNonContact"); return res; } /* -------------------------------------------------------------------------- */ void BemGigipol::optimizeToMeanDisplacement(Real imposed_mean) { Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0); UInt n = surface.size(); UInt size = n * n; // Initial guesses for upper and lower bound Real step_min = -10; Real step_max = 10; // Gaps for upper and lower bound Real gap_min = positiveGapAverage(step_min); Real gap_max = positiveGapAverage(step_max); UInt max_expansion = 8; // Expand bounds if necessary for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10) gap_max = positiveGapAverage(step_max); for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10) gap_min = positiveGapAverage(step_min); Real g = 0.; Real epsilon = 1e-12; Real step = 0; while (fabs(g - target_value) > epsilon) { step = (step_min + step_max) / 2.; g = positiveGapAverage(step); if (g > target_value) step_max = step; else if (g < target_value) step_min = step; else { step_max = step; step_min = step; } } step = (step_min + step_max) / 2.; #pragma omp parallel for for (UInt i = 0 ; i < size ; i++) { gap(i) += step; if (gap(i) < 0) gap(i) = 0; true_displacements(i) = gap(i) + surface(i); } } /* -------------------------------------------------------------------------- */ Real BemGigipol::positiveGapAverage(Real shift) { UInt n = surface.size(); Real res = 0; #pragma omp parallel for reduction(+: res) for (UInt i = 0 ; i < n * n ; i++) { Real shifted_gap = gap(i) + shift; res += shifted_gap * (shifted_gap > 0); } return res / (n * n); } /* -------------------------------------------------------------------------- */ void BemGigipol::updateT(Real G, Real Gold, Real delta){ STARTTIMER("updateT"); unsigned int n = surface.size(); unsigned int size = n*n; Real factor = delta*G/Gold; #pragma omp parallel for for (unsigned int i = 0; i < size; ++i) { if (this->gap(i) == 0.) this->surface_t(i) = 0.; else{ this->surface_t(i) *= factor; this->surface_t(i) += this->search_direction(i); } } STOPTIMER("updateT"); } /* -------------------------------------------------------------------------- */ Real BemGigipol::computeTau() { STARTTIMER("computeOptimalStep"); surface_spectral_r = surface_t; surface_spectral_r.FFTTransform(nthreads); unsigned int n = surface.size(); unsigned int size = n*n; surface_spectral_r(0)=0; #pragma omp parallel for for (unsigned int i = 1; i < size; ++i) { surface_spectral_r(i) /= this->surface_spectral_influence_disp(i); } surface_spectral_r.FFTITransform(surface_r,nthreads); Real rbar = 0; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,rbar) for (unsigned int i = 0; i < size; ++i) { if (this->gap(i) == 0) continue; ++nb_contact; rbar += surface_r(i); } rbar /= nb_contact; surface_r -= rbar; Real tau_sum1 = 0.,tau_sum2 = 0.; #pragma omp parallel for reduction(+: tau_sum1, tau_sum2) for (unsigned int i = 0; i < size; ++i) { if (this->gap(i) == 0) continue; tau_sum1 += this->search_direction(i)*surface_t(i); tau_sum2 += surface_r(i)*surface_t(i); } Real tau = tau_sum1/tau_sum2; STOPTIMER("computeTau"); return tau; } /* -------------------------------------------------------------------------- */ Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) { STARTTIMER("updateDisplacements"); unsigned int n = surface.size(); unsigned int size = n*n; //compute number of interpenetration without contact UInt nb_iol = 0; #pragma omp parallel for reduction(+: nb_iol) for (unsigned int i = 0; i < size; ++i) { if ( this->search_direction(i) <0. && this->gap(i) == 0.){ this->true_displacements(i) -= tau* this->search_direction(i); ++nb_iol; } } Real delta = 0; if (nb_iol > 0) delta = 0.; else delta = 1.; return delta; STOPTIMER("updateDisplacements"); } /* -------------------------------------------------------------------------- */ void BemGigipol::enforceMeanDisplacement(Real mean_displacement) { STARTTIMER("enforceMeanDisplacement"); UInt n = surface.size(); UInt size = n*n; Real moyenne_surface=SurfaceStatistics::computeAverage(this->surface, 0); Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0); Real factor = (mean_displacement-moyenne_surface) / (average-moyenne_surface); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->true_displacements(i) = factor*(this->true_displacements(i)-this->surface(i)) +this->surface(i); } STOPTIMER("enforceMeanDisplacement"); } /* -------------------------------------------------------------------------- */ void BemGigipol::computePressures(Real mean_displacement) { UInt n = surface.size(); UInt size = n*n; Real rho = this->functional->getParameter("rho"); Real surface_energy = this->functional->getParameter("surface_energy"); - this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads); - this->surface_spectral_tractions(0) = 0; - -#pragma omp parallel for - for (UInt i = 1 ; i < size ; ++i) { - this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i); - } - - this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads); - + this->applyInverseInfluenceFunctions(this->true_displacements, this->surface_tractions); Real mini=3000; for (UInt i = 0; i < size; ++i) { - if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho); + if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho); } -this->surface_tractions-=mini; + this->surface_tractions-=mini; } +/* -------------------------------------------------------------------------- */ + __END_TAMAAS__ diff --git a/src/bem_gigipol.hh b/src/bem_gigipol.hh index 890a46b..fec1b37 100644 --- a/src/bem_gigipol.hh +++ b/src/bem_gigipol.hh @@ -1,110 +1,110 @@ /** * * @author Valentine Rey * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_GIGIPOL_H #define BEM_GIGIPOL_H /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemGigipol : public BemFFTBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemGigipol(Surface & p) : BemFFTBase(p), search_direction(p.size(), p.getL()), + spectral_search_direction(p.size(), p.getL()), + mean_disp_search_direction(p.size(), p.getL()), surface_t(p.size(),p.getL()), surface_r(p.size(),p.getL()), surface_spectral_r(p.size(),p.getL()), - spectral_search_direction(p.size(), p.getL()), - mean_disp_search_direction(p.size(), p.getL()), - pold(p.size(),p.getL()){ + pold(p.size(),p.getL()){ e_star = 1.; max_iterations = 10000; } virtual ~BemGigipol(){} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: Real computeEquilibrium(Real epsilon, Real mean_displacement); Real computeStoppingCriterion(); Real computeG(); Real computeMeanPressuresInNonContact(); void optimizeToMeanDisplacement(Real imposed_mean); Real positiveGapAverage(Real shift); void updateT(Real G, Real Gold, Real delta); Real updateDisplacements(Real tau, Real mean_displacement); void enforceMeanDisplacement(Real mean_displacement); void computePressures(Real mean_displacement); Real computeTau() ; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! get the number of iterations UInt getNbIterations() const{return this->nb_iterations;}; //! get the convergence steps of iterations const std::vector & getConvergenceIterations() const{return this->convergence_iterations;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! exploration direction for the Gigi algo Surface search_direction; //! projected exploration direction for the Gigi algo (fourier version) SurfaceComplex spectral_search_direction; //! exploration direction for optimization of mean displacement Surface mean_disp_search_direction; //! exploration direction for the polonski algo Surface surface_t; //! projected exploration direction for the polonski algo Surface surface_r; //! projected exploration direction for the polonski algo (fourier version) SurfaceComplex surface_spectral_r; //! previous pressure Surface pold; }; __END_TAMAAS__ #endif /* BEM_GIGIPOL_H */ diff --git a/src/bem_interface.hh b/src/bem_interface.hh index 20c5b26..104eb57 100644 --- a/src/bem_interface.hh +++ b/src/bem_interface.hh @@ -1,100 +1,102 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_INTERFACE_H #define BEM_INTERFACE_H #include "surface.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemInterface { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemInterface(Surface & surface): surface(surface), spectral_surface(surface.size(),surface.getL()) { this->computeSpectralSurface(); }; virtual ~BemInterface(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: - //! compute the displacement from the pressure field - virtual void computeDisplacements() = 0; + //! compute the displacements from the tractions + virtual void computeDisplacementsFromTractions(){SURFACE_FATAL("TO BE IMPLEMENTED")}; + //! compute the tractions from the displacements + virtual void computeTractionsFromDisplacements(){SURFACE_FATAL("TO BE IMPLEMENTED")}; public: //! compute the spectral surface virtual void computeSpectralSurface(){ surface.FFTTransform(spectral_surface); }; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! return the actual tractions virtual const Surface & getTractions() const = 0; //! return the actual displacements virtual const Surface & getDisplacements() const = 0; //! return the manipulated surface const Surface & getSurface() const{return this->surface;}; //! return the manipulated surface in fourier space const SurfaceComplex & getSpectralSurface() const{return this->spectral_surface;}; //! return the manipulated surface void setSurface(Surface & surface){this->surface = surface;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! the surface profile we work with Surface & surface; //! the surface profile in fourier space SurfaceComplex spectral_surface; }; __END_TAMAAS__ #endif /* BEM_INTERFACE_H */ diff --git a/src/bem_kato.cpp b/src/bem_kato.cpp index f45d90d..5133229 100644 --- a/src/bem_kato.cpp +++ b/src/bem_kato.cpp @@ -1,253 +1,253 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_kato.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemKato::linescan(Real scale,Real pressure){ updatePressure(scale); shiftPressure(pressure); truncatePressure(); this->functional->computeF(); Real res = this->functional->getF(); return res; } /* -------------------------------------------------------------------------- */ Real BemKato::linesearch(Real & hmax, Real fold, Real pressure, int search_flag){ if (search_flag){ Real h = hmax; // Real fold = bem.computeF(); //if (fold == 0) fold = 1e300; Real f = linescan(h,pressure); if (f < fold) return f; while (f > fold){ h *= 0.5; if (h < 1e-3) throw 1; f = linescan(h,pressure); } f = linescan(h,pressure); // if (hmax / h > 10) hmax /=2; return f; } return linescan(hmax,pressure); } /* -------------------------------------------------------------------------- */ Real BemKato::computeEquilibrium(Real epsilon, Real pressure){ // UInt n = surface.size(); // UInt size = n*n; - this->computeDisplacements(); + this->computeDisplacementsFromTractions(); this->functional->computeGradF(); this->backupTractions(); Real f = 1.; Real fPrevious = 1e300; this->nb_iterations = 0; this->convergence_iterations.clear(); while(f > epsilon && this->nb_iterations < this->max_iterations) { fPrevious = f; - this->computeDisplacements(); + this->computeDisplacementsFromTractions(); this->functional->computeGradF(); this->backupTractions(); try { f = linescan(1.,pressure); } catch (int e){ std::cout << " line search failed " << std::endl; f = linescan(1,pressure); nb_iterations = 0; break; } if (nb_iterations % dump_freq == 0) { // std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fold << " " << ((f-fold)/forigin) << std::endl; std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fPrevious << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); this->computeGaps(); return f; } /* -------------------------------------------------------------------------- */ void BemKato::updatePressure(Real scale){ STARTTIMER("updatePressure"); UInt n = surface.size(); UInt size = n*n; const Surface & gradF = this->functional->getGradF(); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions(i) = this->surface_tractions_backup(i) - gradF(i) * scale; } STOPTIMER("updatePressure"); } /* -------------------------------------------------------------------------- */ void BemKato::backupTractions(){ STARTTIMER("switchPressure"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions_backup(i) = this->surface_tractions(i); } STOPTIMER("switchPressure"); } /* -------------------------------------------------------------------------- */ Real BemKato::positivePressure(Real step){ STARTTIMER("positivePressure"); UInt n = surface.size(); UInt size = n*n; Real p_tot = 0.0; #pragma omp parallel for reduction(+:p_tot) for (UInt i = 0; i < size; ++i) { Real sh_press = this->surface_tractions(i) + step; if (sh_press > max_pressure) p_tot += max_pressure; else p_tot += sh_press*(sh_press > 0); } STOPTIMER("positivePressure"); return p_tot/n/n; } /* -------------------------------------------------------------------------- */ void BemKato::shiftPressure(Real required_pressure){ Real step_min = -10; Real step_max = 10; STARTTIMER("shiftPressureInitialGuess"); Real p_max = positivePressure(step_max); Real p_min = positivePressure(step_min); for(UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max*=10) { p_max = positivePressure(step_max); } for(UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min*=10) { p_min = positivePressure(step_min); } Real p = positivePressure(0.0); Real epsilon = 1e-12; STOPTIMER("shiftPressureInitialGuess"); STARTTIMER("shiftPressureDichotomy"); while (fabs(step_min - step_max) > epsilon){ Real step = (step_min + step_max)/2; p = positivePressure(step); if (p > required_pressure) step_max = step; else if (p < required_pressure) step_min = step; else { step_max = step; step_min = step; } } STOPTIMER("shiftPressureDichotomy"); STARTTIMER("shiftPressureTrueShift"); UInt n = surface.size(); UInt size = n*n; //shift the pressure so that satisfies the constraint #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions(i) += (step_max+step_min)/2; } STOPTIMER("shiftPressureTrueShift"); } /* -------------------------------------------------------------------------- */ void BemKato::truncatePressure(){ STARTTIMER("truncatePressure"); UInt n = surface.size(); UInt size = n*n; //shift the pressure so that satisfies the constraint #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) > max_pressure) this->surface_tractions(i) = max_pressure; else this->surface_tractions(i) *= (this->surface_tractions(i) > 0); } STOPTIMER("truncatePressure"); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_penalty.cpp b/src/bem_penalty.cpp index 44e5833..b46bae6 100644 --- a/src/bem_penalty.cpp +++ b/src/bem_penalty.cpp @@ -1,241 +1,218 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_penalty.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemPenalty::computeEquilibrium(Real epsilon, Real mean_displacement, Real penalization_parameter) { this->computeSpectralInfluenceOverDisplacement(); this->search_direction = 0.; this->true_displacements =0; this->true_displacements = SurfaceStatistics::computeMaximum(this->surface); this->computeGaps(); Real f = 1e300; std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl; convergence_iterations.clear(); nb_iterations = 0; max_iterations =5000; while (f> epsilon && nb_iterations++ < max_iterations) { this->functional->computeGradF(); this->computeSearchDirection(mean_displacement,penalization_parameter); Real alpha=0.0001; this->old_displacements=this->true_displacements; this->updateUnknown( alpha, mean_displacement); this->computeGaps(); f = computeStoppingCriterion(); if (nb_iterations % dump_freq == 0) { std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << std::fixed << std::endl; } convergence_iterations.push_back(f); } this->computePressures(mean_displacement); return f; } /* -------------------------------------------------------------------------- */ Real BemPenalty::computeStoppingCriterion() { Real crit=0.; Real disp_norm = 0.; UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for reduction(+:crit, disp_norm) for (UInt i = 0; i < size; ++i) { crit += this->search_direction(i)*this->search_direction(i); disp_norm += (true_displacements(i)*true_displacements(i)); } return crit / disp_norm; } /* -------------------------------------------------------------------------- */ void BemPenalty::computeSearchDirection(Real mean_displacement,Real penalization_parameter) { STARTTIMER("computeOptimalStep"); UInt n = surface.size(); UInt size = n*n; Surface & gradF = this->functional->getGradF(); #pragma omp parallel for for (UInt i = 1; i < size; ++i) { this->search_direction(i)=gradF(i); if (gap(i)<0) { this->search_direction(i)+=penalization_parameter*gap(i); } } } /* -------------------------------------------------------------------------- */ Real BemPenalty::computeOptimalStep() { STARTTIMER("computeOptimalStep"); search_direction.FFTTransform(spectral_search_direction, nthreads); UInt n = surface.size(); UInt size = n*n; spectral_search_direction(0) = 0; #pragma omp parallel for for (UInt i = 1; i < size; ++i) { spectral_search_direction(i) /= this->surface_spectral_influence_disp(i); } // /!\ does not contain the spectral search direction anymore spectral_search_direction.FFTITransform(nthreads); // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0); //spectral_search_direction -= average; Real numerator = 0., denominator = 0.; #pragma omp parallel for reduction(+: numerator, denominator) for (UInt i = 0; i < size; ++i) { numerator += search_direction(i) * search_direction(i); denominator += spectral_search_direction(i).real() * search_direction(i); } Real alpha = numerator / denominator; STOPTIMER("computeOptimalStep"); return alpha; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void BemPenalty::updateUnknown(Real alpha, Real mean_displacement) { STARTTIMER("updateDisplacements"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->true_displacements(i) -= alpha*this->search_direction(i); } Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0); for (UInt i = 0; i < size; ++i) { this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement; } STOPTIMER("updateDisplacements"); } /* -------------------------------------------------------------------------- */ void BemPenalty::computePressures(Real mean_displacement) { - UInt n = surface.size(); - UInt size = n*n; - - Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0); - - this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads); - this->surface_spectral_tractions(0) = 0; - -#pragma omp parallel for - for (UInt i = 1 ; i < size ; ++i) { - this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i); - } - - this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads); - - - - - - - - - - + this->computeTractionsFromDisplacements(); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp index a5fd79d..26ec1d3 100644 --- a/src/bem_polonski.cpp +++ b/src/bem_polonski.cpp @@ -1,314 +1,314 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_polonski.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ BemPolonski::BemPolonski(Surface & p) : BemFFTBase(p), surface_t(p.size(),p.getL()), surface_r(p.size(),p.getL()), surface_spectral_r(p.size(),p.getL()), surface_r_transform(p.size(), surface_r, surface_spectral_r), pold(p.size(),p.getL()){ e_star = 1.; max_iterations = 10000; } BemPolonski::~BemPolonski() {} Real BemPolonski::computeEquilibrium(Real epsilon, Real pressure){ this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; this->surface_spectral_r = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions); if (current_pressure <= 0.) surface_tractions = pressure; this->enforcePressureBalance(pressure); convergence_iterations.clear(); nb_iterations = 0; while(f > epsilon && nb_iterations < max_iterations) { fPrevious = f; - this->computeDisplacements(); + this->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); //std::cout << "G vaut "<< G<< std::endl; this->updateT(G,Gold,delta); Real tau = this->computeTau(); pold = this->surface_tractions; Gold = G; delta = this->updateTractions(tau); //Projection on admissible space this->enforcePressureBalance(pressure); f = this->computeF(); if (nb_iterations % dump_freq == 0){ Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions); std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fPrevious << " " << A << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } return f; } /* -------------------------------------------------------------------------- */ void BemPolonski::computeGaps() { UInt size = surface.size(); #pragma omp parallel for for (UInt i=0; i < size*size; i++) { gap(i) = surface_displacements(i) - surface(i); } } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeMeanGapsInContact(){ STARTTIMER("computeMeanGapsInContact"); UInt n = surface.size(); UInt size = n*n; Real res = 0.; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; ++nb_contact; res += this->gap(i); } res /= nb_contact; STOPTIMER("computeMeanGapsInContact"); return res; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeG(){ STARTTIMER("computeG"); UInt n = surface.size(); UInt size = n*n; Real res = 0.; #pragma omp parallel for reduction(+: res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; Real val = this->gap(i); res += val*val; } STOPTIMER("computeG"); return res; } /* -------------------------------------------------------------------------- */ void BemPolonski::updateT(Real G, Real Gold, Real delta){ STARTTIMER("updateT"); UInt n = surface.size(); UInt size = n*n; Real factor = delta*G/Gold; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.; else{ this->surface_t(i) *= factor; this->surface_t(i) += this->gap(i); } } STOPTIMER("updateT"); } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeTau(){ STARTTIMER("computeTau"); surface_r = surface_t; surface_r_transform.forwardTransform(); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { surface_spectral_r(i) *= this->surface_spectral_influence_disp(i); } surface_r_transform.backwardTransform(); Real rbar = 0; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,rbar) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; ++nb_contact; rbar += surface_r(i); } rbar /= nb_contact; surface_r -= rbar; Real tau_sum1 = 0.,tau_sum2 = 0.; #pragma omp parallel for reduction(+: tau_sum1, tau_sum2) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; tau_sum1 += this->gap(i)*surface_t(i); tau_sum2 += surface_r(i)*surface_t(i); } Real tau = tau_sum1/tau_sum2; STOPTIMER("computeTau"); return tau; } /* -------------------------------------------------------------------------- */ Real BemPolonski::updateTractions(Real tau){ STARTTIMER("updateTractions"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; this->surface_tractions(i) -= tau*this->surface_t(i); if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0; } //compute number of interpenetration without contact UInt nb_iol = 0; #pragma omp parallel for reduction(+: nb_iol) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.){ this->surface_tractions(i) -= tau*this->gap(i); ++nb_iol; } } Real delta = 0; if (nb_iol > 0) delta = 0.; else delta = 1.; STOPTIMER("updateTractions"); return delta; } /* -------------------------------------------------------------------------- */ void BemPolonski::enforcePressureBalance(Real applied_pressure){ STARTTIMER("enforcePressureBalance"); UInt n = surface.size(); UInt size = n*n; Real pressure = 0.; #pragma omp parallel for reduction(+: pressure) for (UInt i = 0; i < size; ++i) { pressure += this->surface_tractions(i); } pressure *= 1./size; for (UInt i = 0; i < size; ++i) { this->surface_tractions(i) *= applied_pressure/pressure; } STOPTIMER("enforcePressureBalance"); } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeF() { UInt n = surface.size(); UInt size = n*n; Real res = 0; Real d_max = SurfaceStatistics::computeMaximum(surface_displacements); Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions)); Real s_max = SurfaceStatistics::computeMaximum(surface); #pragma omp parallel for reduction(+:res) for (UInt i = 0; i < size; ++i) { res += std::abs(surface_tractions(i) * (surface_displacements(i)-d_max-surface(i)+s_max)); // res += // this->surface_tractions[i].real() // *(surface_displacements[i].real() - surface[i].real()); } return res/std::abs(t_sum); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_uzawa.cpp b/src/bem_uzawa.cpp index 093f9ce..97130b2 100644 --- a/src/bem_uzawa.cpp +++ b/src/bem_uzawa.cpp @@ -1,320 +1,303 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_uzawa.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemUzawa::computeEquilibrium(Real epsilon, - Real mean_displacement, Real penalization_parameter) { + Real mean_displacement, Real penalization_parameter) { this->computeSpectralInfluenceOverDisplacement(); - this->search_direction = 0.; - this->true_displacements =0; - this->true_displacements =mean_displacement; + this->search_direction = 0.; + this->true_displacements = 0; + this->true_displacements = mean_displacement; this->computeGaps(); - this->multipliers =0; + this->multipliers = 0; - Real f_2 = 1e300; - Real new_pen = 1e300; - Real old_pen = 1e300; + Real f_2 = 1e300; + Real new_pen = 1e300; + Real old_pen = 1e300; - convergence_iterations.clear(); - UInt nb_iterations_1 = 0; - UInt nb_iterations_2 = 0; + convergence_iterations.clear(); + UInt nb_iterations_1 = 0; + UInt nb_iterations_2 = 0; - Real max_iterations_1 =500; - Real max_iterations_2 =1000; + Real max_iterations_1 = 500; + Real max_iterations_2 = 1000; while (new_pen> epsilon && nb_iterations_1++ < max_iterations_1) { - std::cout << "iter ext "<< nb_iterations_1 << std::endl; + std::cout << "iter ext "<< nb_iterations_1 << std::endl; while (f_2> 1e-12 && nb_iterations_2++ < max_iterations_2) { - std::cout << "iter int "<< nb_iterations_2 << std::endl; - this->functional->computeGradF(); - this->computeSearchDirection(mean_displacement,penalization_parameter); + std::cout << "iter int "<< nb_iterations_2 << std::endl; + this->functional->computeGradF(); + this->computeSearchDirection(mean_displacement,penalization_parameter); - Real alpha=1/(10*penalization_parameter); - this->old_displacements=this->true_displacements; - this->updateUnknown( alpha, mean_displacement); - this->computeGaps(); + Real alpha=1/(10*penalization_parameter); + this->old_displacements=this->true_displacements; + this->updateUnknown( alpha, mean_displacement); + this->computeGaps(); - f_2 = computeStoppingCriterion(); + f_2 = computeStoppingCriterion(); - if (nb_iterations_2 % dump_freq == 0) { - std::cout << std::scientific << std::setprecision(10) - << nb_iterations << " " - << f_2 << std::fixed << std::endl; + if (nb_iterations_2 % dump_freq == 0) { + std::cout << std::scientific << std::setprecision(10) + << nb_iterations << " " + << f_2 << std::fixed << std::endl; + } } - } this->updateMultipliers(penalization_parameter); old_pen=computeInterpenetration(this->old_displacements); new_pen=computeInterpenetration(this->true_displacements); std::cout << "old penetration is "<< old_pen << std::endl; std::cout << "new penetration is "<< new_pen << std::endl; penalization_parameter=this->updatePenalization(penalization_parameter, old_pen, new_pen); //to avoid ill conditioned system if (penalization_parameter>1.0e8) new_pen=0.; nb_iterations_2 = 0; - f_2 = 1e300; + f_2 = 1e300; std::cout << "penalization is "<< penalization_parameter << std::endl; - } + } this->computeGaps(); - this->computePressures(mean_displacement); + this->computePressures(); return new_pen; } /* -------------------------------------------------------------------------- */ Real BemUzawa::computeStoppingCriterion() { - Real crit=0.; - Real disp_norm = 0.; + Real crit=0.; + Real disp_norm = 0.; - UInt n = surface.size(); - UInt size = n*n; + UInt n = surface.size(); + UInt size = n*n; - #pragma omp parallel for reduction(+:crit, disp_norm) - for (UInt i = 0; i < size; ++i) { +#pragma omp parallel for reduction(+:crit, disp_norm) + for (UInt i = 0; i < size; ++i) { - crit += (this->search_direction(i))*(this->search_direction(i)); - disp_norm += (true_displacements(i)*true_displacements(i)); - } + crit += (this->search_direction(i))*(this->search_direction(i)); + disp_norm += (true_displacements(i)*true_displacements(i)); + } - return crit / disp_norm; + return crit / disp_norm; } /* -------------------------------------------------------------------------- */ void BemUzawa::computeSearchDirection(Real mean_displacement,Real penalization_parameter) { STARTTIMER("computeOptimalStep"); UInt n = surface.size(); UInt size = n*n; - Surface & gradF = this->functional->getGradF(); + Surface & gradF = this->functional->getGradF(); #pragma omp parallel for for (UInt i = 1; i < size; ++i) { - this->search_direction(i)=gradF(i); + this->search_direction(i)=gradF(i); - if (this->multipliers(i)+penalization_parameter*gap(i)<=0){ + if (this->multipliers(i)+penalization_parameter*gap(i)<=0){ - this->search_direction(i)=this->search_direction(i)+ this->multipliers(i)+penalization_parameter*gap(i); + this->search_direction(i)=this->search_direction(i)+ this->multipliers(i)+penalization_parameter*gap(i); } -} + } } /* -------------------------------------------------------------------------- */ Real BemUzawa::computeOptimalStep() { STARTTIMER("computeOptimalStep"); search_direction.FFTTransform(spectral_search_direction, nthreads); UInt n = surface.size(); UInt size = n*n; spectral_search_direction(0) = 0; #pragma omp parallel for - for (UInt i = 1; i < size; ++i) { + for (UInt i = 1; i < size; ++i) { spectral_search_direction(i) /= this->surface_spectral_influence_disp(i); } - // /!\ does not contain the spectral search direction anymore + // /!\ does not contain the spectral search direction anymore spectral_search_direction.FFTITransform(nthreads); - // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0); -//spectral_search_direction -= average; + // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0); + //spectral_search_direction -= average; - Real numerator = 0., denominator = 0.; + Real numerator = 0., denominator = 0.; #pragma omp parallel for reduction(+: numerator, denominator) - for (UInt i = 0; i < size; ++i) { - numerator += search_direction(i) * search_direction(i); + for (UInt i = 0; i < size; ++i) { + numerator += search_direction(i) * search_direction(i); denominator += spectral_search_direction(i).real() * search_direction(i); - } + } - Real alpha = numerator / denominator; + Real alpha = numerator / denominator; STOPTIMER("computeOptimalStep"); - return alpha; + return alpha; } /* -------------------------------------------------------------------------- */ void BemUzawa::updateMultipliers(Real penalization_parameter) { STARTTIMER("updateMultipliers"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { - if (this->multipliers(i)+penalization_parameter*gap(i)>0.){ - this->multipliers(i) = 0; - } - else{ - this->multipliers(i)+=penalization_parameter*gap(i); - } + if (this->multipliers(i)+penalization_parameter*gap(i)>0.){ + this->multipliers(i) = 0; + } + else{ + this->multipliers(i)+=penalization_parameter*gap(i); + } } STOPTIMER("updateMultipliers"); } /* -------------------------------------------------------------------------- */ Real BemUzawa::updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen) { STARTTIMER("updatePenalization"); - if (new_pen>0.25*old_pen){ - penalization_parameter*=5; + if (new_pen>0.25*old_pen){ + penalization_parameter*=5; } - return penalization_parameter; + return penalization_parameter; STOPTIMER("updatePenalization"); } /* -------------------------------------------------------------------------- */ Real BemUzawa::computeInterpenetration( Surface & displacements ) { STARTTIMER("updateMultipliers"); UInt n = surface.size(); UInt size = n*n; Real res=0.; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { - if (0>this->true_displacements(i)-surface(i)){ - res+=(this->true_displacements(i)-surface(i))*(this->true_displacements(i)-surface(i)); - } + if (0>this->true_displacements(i)-surface(i)){ + res+=(this->true_displacements(i)-surface(i))*(this->true_displacements(i)-surface(i)); + } -} + } - return sqrt(res); + return sqrt(res); STOPTIMER("updateMultipliers"); } /* -------------------------------------------------------------------------- */ void BemUzawa::updateUnknown(Real alpha, Real mean_displacement) { STARTTIMER("updateDisplacements"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->true_displacements(i) -= alpha*this->search_direction(i); } -Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0); + Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0); -for (UInt i = 0; i < size; ++i) { - this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement; + for (UInt i = 0; i < size; ++i) { + this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement; -} + } STOPTIMER("updateDisplacements"); } /* -------------------------------------------------------------------------- */ -void BemUzawa::computePressures(Real mean_displacement) { - UInt n = surface.size(); - UInt size = n*n; - - Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0); - - this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads); - this->surface_spectral_tractions(0) = 0; - -#pragma omp parallel for - for (UInt i = 1 ; i < size ; ++i) { - this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i); - } - - this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads); - - - - +void BemUzawa::computePressures() { + this->applyInverseInfluenceFunctions(this->true_displacements, this->surface_tractions); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_uzawa.hh b/src/bem_uzawa.hh index 474a720..2ec22d2 100644 --- a/src/bem_uzawa.hh +++ b/src/bem_uzawa.hh @@ -1,102 +1,100 @@ /** * * @author Valentine Rey * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_UZAWA_H #define BEM_UZAWA_H /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemUzawa : public BemFFTBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemUzawa(Surface & p) : BemFFTBase(p), search_direction(p.size(), p.getL()), spectral_search_direction(p.size(), p.getL()), multipliers(p.size(),p.getL()), surface_r(p.size(),p.getL()), surface_spectral_r(p.size(),p.getL()), mean_disp_search_direction(p.size(), p.getL()) { e_star = 1.; max_iterations = 10000; }; virtual ~BemUzawa(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: //! compute the equilibrium situation Real computeEquilibrium(Real epsilon,Real pressure,Real penalization_parameter); void computeSearchDirection(Real mean_displacement,Real penalization_parameter); - void updateMultipliers(Real penalization_parameter); - Real updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen); - Real computeInterpenetration( Surface & displacements ) ; + void updateMultipliers(Real penalization_parameter); + Real updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen); + Real computeInterpenetration( Surface & displacements ) ; void updateUnknown(Real tau, Real mean_displacement); Real computeStoppingCriterion(); - void computePressures(Real mean_displacement); - Real computeOptimalStep(); - - + void computePressures(); + Real computeOptimalStep(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! get the number of iterations UInt getNbIterations() const{return this->nb_iterations;}; //! get the convergence steps of iterations const std::vector & getConvergenceIterations() const{return this->convergence_iterations;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! exploration direction for the Gigi algo Surface search_direction; //! projected exploration direction for the Gigi algo (fourier version) SurfaceComplex spectral_search_direction; //! exploration direction for the polonski algo Surface multipliers; //! projected exploration direction for the polonski algo Surface surface_r; //! projected exploration direction for the polonski algo (fourier version) SurfaceComplex surface_spectral_r; //! exploration direction for optimization of mean displacement Surface mean_disp_search_direction; }; __END_TAMAAS__ #endif /* BEM_Uzawa_H */ diff --git a/src/fft_plan_manager.cpp b/src/fft_plan_manager.cpp new file mode 100644 index 0000000..2e15d15 --- /dev/null +++ b/src/fft_plan_manager.cpp @@ -0,0 +1,89 @@ +/** + * + * @author Guillaume Anciaux + * + * @section LICENSE + * + * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de + * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des + * Solides) + * + * Tamaas is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Tamaas 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 Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Tamaas. If not, see . + * + */ +/* -------------------------------------------------------------------------- */ +#include "fft_plan_manager.hh" +/* -------------------------------------------------------------------------- */ + + +__BEGIN_TAMAAS__ +/* -------------------------------------------------------------------------- */ + +FFTPlanManager::FFTPlanManager(){ +} +/* -------------------------------------------------------------------------- */ + +FFTPlanManager::~FFTPlanManager(){ +} +/* -------------------------------------------------------------------------- */ + + +FFTPlanManager & FFTPlanManager::get(){ + if (FFTPlanManager::singleton == NULL) + FFTPlanManager::singleton = new FFTPlanManager(); + + return *FFTPlanManager::singleton; +} +/* -------------------------------------------------------------------------- */ + +FFTransform & FFTPlanManager::createPlan(Surface & input, + SurfaceComplex & output){ + + auto index = std::make_pair(const_cast(input.getInternalData()), + const_cast(output.getInternalData())); + + auto it = plans.find(index); + auto end = plans.end(); + + if (it == end){ + plans[index] = new FFTransform(input.size(),input,output); + } + + return *plans[index]; +} + +/* -------------------------------------------------------------------------- */ + +void FFTPlanManager::destroyPlan(Surface & input, + SurfaceComplex & output){ + + auto index = std::make_pair(const_cast(input.getInternalData()), + const_cast(output.getInternalData())); + + auto it = plans.find(index); + auto end = plans.end(); + + if (it != end){ + delete (plans[index]); + plans.erase(index); + } +} + +/* -------------------------------------------------------------------------- */ +FFTPlanManager * FFTPlanManager::singleton = NULL; +/* -------------------------------------------------------------------------- */ + + + +__END_TAMAAS__ diff --git a/src/bem_interface.hh b/src/fft_plan_manager.hh similarity index 56% copy from src/bem_interface.hh copy to src/fft_plan_manager.hh index 20c5b26..90d5bde 100644 --- a/src/bem_interface.hh +++ b/src/fft_plan_manager.hh @@ -1,100 +1,77 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ - /* -------------------------------------------------------------------------- */ -#ifndef BEM_INTERFACE_H -#define BEM_INTERFACE_H - -#include "surface.hh" +#ifndef FFT_PLAN_MANAGER_H +#define FFT_PLAN_MANAGER_H +/* -------------------------------------------------------------------------- */ +#include +#include "fftransform.hh" /* -------------------------------------------------------------------------- */ - __BEGIN_TAMAAS__ +/* -------------------------------------------------------------------------- */ +class FFTPlanManager { - -class BemInterface { - + /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ -public: - - BemInterface(Surface & surface): - surface(surface), - spectral_surface(surface.size(),surface.getL()) { - this->computeSpectralSurface(); - }; - virtual ~BemInterface(){}; - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - protected: - //! compute the displacement from the pressure field - virtual void computeDisplacements() = 0; -public: + FFTPlanManager(); - //! compute the spectral surface - virtual void computeSpectralSurface(){ - surface.FFTTransform(spectral_surface); - }; +public: + virtual ~FFTPlanManager(); +public: + /* ------------------------------------------------------------------------ */ - /* Accessors */ + /* Methods */ /* ------------------------------------------------------------------------ */ -public: - - //! return the actual tractions - virtual const Surface & getTractions() const = 0; - //! return the actual displacements - virtual const Surface & getDisplacements() const = 0; - //! return the manipulated surface - const Surface & getSurface() const{return this->surface;}; - //! return the manipulated surface in fourier space - const SurfaceComplex & getSpectralSurface() const{return this->spectral_surface;}; - //! return the manipulated surface - void setSurface(Surface & surface){this->surface = surface;}; - + static FFTPlanManager & get(); + FFTransform & createPlan(Surface & input, SurfaceComplex & output); + void destroyPlan(Surface & input, SurfaceComplex & output); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ -protected: +private: + + std::map,FFTransform *> plans; - //! the surface profile we work with - Surface & surface; - //! the surface profile in fourier space - SurfaceComplex spectral_surface; + static FFTPlanManager * singleton; }; +/* -------------------------------------------------------------------------- */ __END_TAMAAS__ +/* -------------------------------------------------------------------------- */ + -#endif /* BEM_INTERFACE_H */ +#endif /* FFT_PLAN_MANAGER_H */ diff --git a/src/map_2d.cpp b/src/map_2d.cpp index 9ff44ae..2650acb 100644 --- a/src/map_2d.cpp +++ b/src/map_2d.cpp @@ -1,246 +1,246 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "map_2d.hh" #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template void Map2d::rescaleHeights(Real C){ /* compute actual RMS */ for (UInt i = 0 ; i < n[0] ; ++i) for (UInt j = 0 ; j < n[1] ; ++j){ this->at(i,j) *= C; } } /* -------------------------------------------------------------------------- */ template void Map2d::setGridSize(UInt n0, UInt n1){ n[0] = n0; n[1] = n1; data.resize(n[0]*n[1]); data.assign(n[0]*n[1],T(0.)); } /* -------------------------------------------------------------------------- */ template void Map2d::setGridSize(const UInt n[2]){ setGridSize(n[0],n[1]); } /* -------------------------------------------------------------------------- */ template void Map2d::operator=(const T & e){ for (UInt i = 0 ; i < n[0]*n[1] ; ++i){ this->at(i) = e; } } /* -------------------------------------------------------------------------- */ template Map2d::Map2d(UInt n0, UInt n1, Real L0, Real L1) { setGridSize(n0,n1); this->L[0] = L0; this->L[1] = L1; this->plan = NULL; generateFFTPlan(); } /* -------------------------------------------------------------------------- */ template Map2d::Map2d(const UInt * n, const Real *L) { setGridSize(n[0],n[1]); this->L[0] = L[0]; this->L[1] = L[1]; this->plan = NULL; generateFFTPlan(); } /* -------------------------------------------------------------------------- */ template Map2d::~Map2d() { destroyFFTPlan(); } /* -------------------------------------------------------------------------- */ template void Map2d::generateFFTPlan() { // this should not do anything } /* -------------------------------------------------------------------------- */ template <> void Map2d::generateFFTPlan() { fftw_complex * in = (fftw_complex*) &data[0]; this->plan = fftw_plan_dft_2d(n[0], n[1], in, in, FFTW_FORWARD, FFTW_ESTIMATE); } /* -------------------------------------------------------------------------- */ template void Map2d::destroyFFTPlan() { if (this->plan) { fftw_destroy_plan(this->plan); this->plan = NULL; } } /* -------------------------------------------------------------------------- */ template void Map2d::FFTTransform(int nthreads){ SURFACE_FATAL("impossible situation"); throw; } /* -------------------------------------------------------------------------- */ template -void Map2d::FFTTransform(Map2d > & output, int nthreads){ +void Map2d::FFTTransform(Map2d > & output, int nthreads)const{ output = *this; output.FFTTransform(nthreads); } /* -------------------------------------------------------------------------- */ template -void Map2d::FFTITransform(Map2d > & output, int nthreads){ +void Map2d::FFTITransform(Map2d > & output, int nthreads)const{ output = *this; output.FFTITransform(nthreads); } /* -------------------------------------------------------------------------- */ template -void Map2d::FFTITransform(Map2d & output, int nthreads){ +void Map2d::FFTITransform(Map2d & output, int nthreads) const{ const UInt * n = this->sizes(); const Real * Ls = this->getLs(); Map2d > _output(n,Ls); _output = *this; _output.FFTITransform(nthreads); for (UInt i = 0; i < n[0]*n[1]; i++) { output(i) = _output(i).real(); } } /* -------------------------------------------------------------------------- */ template <> void Map2d::FFTTransform(int nthreads){ STARTTIMER("FFTW"); fftw_execute(this->plan); STOPTIMER("FFTW"); } /* -------------------------------------------------------------------------- */ template void Map2d::FFTITransform(int nthreads){ SURFACE_FATAL("impossible situation"); throw; } /* -------------------------------------------------------------------------- */ template <> void Map2d::FFTITransform(int nthreads){ #pragma omp parallel for for (UInt i = 0; i < n[0]*n[1]; ++i) { data[i] = std::conj(data[i]); } FFTTransform(nthreads); #pragma omp parallel for for (UInt i = 0; i < n[0]*n[1]; ++i) { data[i] = std::conj(data[i]); data[i] *= 1./(n[0]*n[1]); } } /* -------------------------------------------------------------------------- */ template void Map2d::dumpToTextFile(std::string filename) { UInt n1 = this->size(0); UInt n2 = this->size(1); std::cout << "Writing surface file " << filename << std::endl; std::ofstream file(filename.c_str()); if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename); file.precision(15); for (UInt i = 0; i < n1; i++) { for (UInt j=0; j< n2; j++) { file << i << " " << j << " "; this->writeValue(file,this->at(i,j)); file << std::endl; } file << std::endl; } file.close(); } /* -------------------------------------------------------------------------- */ template void Map2d::writeValue(std::ostream & os, const T & val) { os << val; } /* -------------------------------------------------------------------------- */ template <> void Map2d::writeValue(std::ostream & os, const complex & val) { os << std::scientific << val.real(); os << " " << val.imag(); } /* -------------------------------------------------------------------------- */ template class Map2d; template class Map2d; template class Map2d; template class Map2d; template class Map2d; template class Map2d; __END_TAMAAS__ diff --git a/src/map_2d.hh b/src/map_2d.hh index a7487ec..6dd4270 100644 --- a/src/map_2d.hh +++ b/src/map_2d.hh @@ -1,482 +1,591 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __MAP_2D_HH__ #define __MAP_2D_HH__ /* -------------------------------------------------------------------------- */ #include #include #include #include "tamaas.hh" #include #include "my_vector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template class Map2d { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ // class iterator; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.); Map2d(const UInt * n, const Real * L); + virtual ~Map2d(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template std::vector > getNeighborIndexes(UInt i, UInt j); template std::vector > getNeighborIndexes(std::pair i); // //! return an iterator over the neighbors // iterator beginNeighborhood(); // //! return the end iterator over the neighbors // iterator endNeighborhood(); //! generate fft plan void generateFFTPlan(); //! destroy fft plan void destroyFFTPlan(); //! perform fft transform - void FFTTransform(Map2d > & output,int nthreads = 1); + void FFTTransform(Map2d > & output,int nthreads = 1) const; //! perform inverse fft transform - void FFTITransform(Map2d > & output, int nthreads = 1); + void FFTITransform(Map2d > & output, int nthreads = 1) const; //! perform inverse fft transform - void FFTITransform(Map2d & output, int nthreads = 1); + void FFTITransform(Map2d & output, int nthreads = 1) const; //! perform fft transform void FFTTransform(int nthreads = 1); // //! inverse fft transform void FFTITransform(int nthreads = 1); //compute compute component by component product and put result in object template void fullProduct(Map2d & s1,Map2d &s2); //! set the grid size void setGridSize(const UInt n[2]); //! set the grid size void setGridSize(UInt n0,UInt n1); //! rescale height profile void rescaleHeights(Real factor); //! Write heights array to text file usable with plot programs void dumpToTextFile(std::string filename); + //! copy method + template class surf, typename T2> + void copy(const surf & s); // set all the members of the map to a given value void operator=(const T & e); - + // multiply the map by scalar + template + void operator*=(const T2 & e); + // divide the map by scalar + template + void operator/=(const T2 & e); + // multiply the map by scalar + void operator*=(const T & e); + // divide the map by scalar + void operator/=(const T & e); + // multiply the map with another map, component by component + template class surf, typename T2> + void operator*=(const surf & e); + // divide the map with another map, component by component + template class surf, typename T2> + void operator/=(const surf & e); + // copy the map passed as parameter template class surf, typename T2> void operator=(const surf & s); + // temporary hack + void operator=(const Map2d & s); + // copy the map passed as const parameter template class surf, typename T2> void operator=(surf & s); /* -------------------------------------------------------------------------- */ private: //! write an entry to a stringstream void writeValue(std::ostream & os, const T & val); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: // access to the i,j-th value in the surface inline void operator+=(const Map2d & s); // access to the i,j-th value in the surface inline void operator+=(const T & val); // access to the i,j-th value in the surface inline void operator-=(const T & val); // access to the i,j-th value in the surface inline T & operator()(UInt i,UInt j); // access to the i,j-th value in the surface inline T & operator()(std::pair i); // access to the i,j-th value in the surface inline const T & operator()(UInt i,UInt j)const; // access to the i,j-th value in the surface inline const T & operator()(std::pair i)const; // access to the i,j-th value in the surface inline T & at(UInt i,UInt j); // access to the i-th value in the array inline T & at(UInt i); // access to the i,j-th value in the surface inline T & at(std::pair i); // access to the i,j-th value in the surface inline const T & at(UInt i,UInt j) const; // access to the i,j-th value in the surface inline const T & at(std::pair i) const; // access to the i,j-th value in the surface inline const T & at(UInt i) const; // access to the i-th value in the array inline T & operator()(int i); // access to the i-th value in the array // inline T & operator[](std::pair i); // access to the i-th value in the array inline const T & operator()(int i) const; //! get lateral number of discretization points inline UInt size(UInt dir)const {return n[dir];}; //! get both lateral number of discretization points inline const UInt* sizes()const {return n;}; //! get the lateral size of the sample inline const Real * getLs() const {return L;}; //! get the lateral size of the sample inline Real getL(UInt dir) const {return L[dir];}; //! get the lateral size of the sample inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];}; //! get the lateral size of the sample inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;}; const T * getInternalData() const{return data.getPtr();}; //! Write heights array to text file usable with plot programs template void dumpToTextFile(std::string filename) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! number of points in each direction UInt n[2]; //! size in each direction Real L[2]; //! internal storage of data MyVector data; //! fftw plan fftw_plan plan; }; /* -------------------------------------------------------------------------- */ //#define _at(I,J,N) I+J*N[1] #define _at(I,J,N) I*N[1]+J /* -------------------------------------------------------------------------- */ template inline T & Map2d::at(UInt i,UInt j){ if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access " << i << " " << j << " " << n[0] << " " << n[1] << " " << _at(i,j,n)); return data[_at(i,j,n)]; } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::at(UInt i,UInt j)const { if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access " << i << " " << j << " " << n[0] << " " << n[1] << " " << _at(i,j,n)); return data[_at(i,j,n)]; } /* -------------------------------------------------------------------------- */ #undef _at /* -------------------------------------------------------------------------- */ template inline T & Map2d::at(UInt i){ return data[i]; } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::at(UInt i) const{ return data[i]; } /* -------------------------------------------------------------------------- */ template inline T & Map2d::operator()(UInt i,UInt j){ return this->at(i,j); } /* -------------------------------------------------------------------------- */ template inline T & Map2d::at(std::pair i){ return this->at(i.first,i.second); } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::at(std::pair i)const { return this->at(i.first,i.second); } /* -------------------------------------------------------------------------- */ template inline T & Map2d::operator()(std::pair i){ return this->at(i.first,i.second); } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::operator()(UInt i,UInt j)const{ return this->at(i,j); } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::operator()(std::pair i)const { return this->at(i.first,i.second); } /* -------------------------------------------------------------------------- */ template inline T & Map2d::operator()(int i){ return data[i]; } /* -------------------------------------------------------------------------- */ template inline const T & Map2d::operator()(int i) const{ return data[i]; } /* -------------------------------------------------------------------------- */ // template // inline T & Map2d::operator[](std::pair c){ // return this->at(c.first,c.second); // } // /* -------------------------------------------------------------------------- */ template void Map2d::operator-=(const T & e){ for (UInt i = 0 ; i < n[0]*n[1] ; ++i){ this->at(i) -= e; } } /* -------------------------------------------------------------------------- */ template void Map2d::operator+=(const T & e){ for (UInt i = 0 ; i < n[0]*n[1] ; ++i){ this->at(i) += e; } } /* -------------------------------------------------------------------------- */ template void Map2d::operator+=(const Map2d & s){ for (UInt i = 0 ; i < n[0]*n[1] ; ++i){ this->at(i) += s.at(i); } } /* -------------------------------------------------------------------------- */ template template std::vector > Map2d::getNeighborIndexes(UInt i, UInt j){ std::vector > vres; std::pair res; for (int ii = -1; ii < 2; ++ii) { res.first = ii+static_cast(i); if (res.first < 0 && !periodic) continue; if (res.first >= (int)n[0] && !periodic) continue; if (res.first < 0 && periodic) res.first += n[0]; if (res.first >= (int)n[0] && periodic) res.first -= n[0]; for (int jj = -1; jj < 2; ++jj) { res.second = jj+static_cast(j); if (res.second < 0 && !periodic) continue; if (res.second >= (int)n[1] && !periodic) continue; if (res.second < 0 && periodic) res.second += n[1]; if (res.second >= (int)n[1] && periodic) res.second -= n[1]; if (ii == 0 && jj == 0) continue; { UInt i = res.first; UInt j = res.second; if (i+j*n[0] >= n[0]*n[1]) SURFACE_FATAL("out of bound access " << i << " " << j << " " << n[0] << " " << n[1] << " " << i+j*n[0]); } vres.push_back(res); } } return vres; } /* -------------------------------------------------------------------------- */ template template std::vector > Map2d ::getNeighborIndexes(std::pair i){ return this->getNeighborIndexes(i.first,i.second); } /* -------------------------------------------------------------------------- */ template template void Map2d::fullProduct(Map2d & s1,Map2d & s2){ auto n1 = s1.sizes(); auto n2 = s2.sizes(); if (n1[0] != n2[0] || n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices"); if (n[0] != n1[0] || n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices") for (UInt i = 0 ; i < n[0] ; ++i) for (UInt j = 0 ; j < n[1] ; ++j){ this->at(i,j) = s1(i,j)* s2(i,j); } } /* -------------------------------------------------------------------------- */ template template void Map2d::dumpToTextFile(std::string filename) const { UInt m = this->n[0]; UInt n = this->n[1]; std::cout << "Writing surface file " << filename << std::endl; std::ofstream file(filename.c_str()); if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename); file.precision(15); for (UInt i = 0; i < m; i++) { for (UInt j=0; j< n; j++) { if (!consider_zeros && std::norm(this->at(i,j)) < zero_threshold*zero_threshold) continue; file << i << " " << j << " "; file << std::scientific << this->at(i,j); file << std::endl; } file << std::endl; } file.close(); } // template // class Map2d::iterator { // public: // iterator(std::vector & data, // std::queue > & indices): // data(data),indices(indices){ // } // inline void operator++(){ // indices.pop(); // } // inline T & operator *(){ // return data[indices.front]; // } // bool operator != (Map2d::iterator & it){ // return (indices.front == it.indices.front); // } // private: // std::vector & data; // std::queue > indices; // bool finished; // }; /* -------------------------------------------------------------------------- */ - template template class surf, typename T2> -void Map2d::operator=(const surf & s){ +void Map2d::copy(const surf & s){ this->setGridSize(s.sizes()); UInt m = this->sizes()[0]; UInt n = this->sizes()[1]; #pragma omp parallel for for (UInt i = 0; i < m*n; i++) { this->at(i) = s(i); } } /* -------------------------------------------------------------------------- */ template template class surf, typename T2> void Map2d::operator=(surf & s){ - this->setGridSize(s.sizes()); - UInt m = this->sizes()[0]; - UInt n = this->sizes()[1]; + this->copy(s); +} +/* -------------------------------------------------------------------------- */ +template +template class surf, typename T2> +void Map2d::operator=(const surf & s){ + this->copy(s); +} + +/* -------------------------------------------------------------------------- */ + +template +void Map2d::operator=(const Map2d & s){ + SURFACE_FATAL("Cannot copy a complex surface into incompatible type"); +} +template <> +inline void Map2d::operator=(const Map2d & s){ + this->copy(s); +} + +/* -------------------------------------------------------------------------- */ +template +template +void Map2d::operator*=(const T2 & e){ + + UInt sz = this->n[0]*this->n[1]; #pragma omp parallel for - for (UInt i = 0; i < m*n; i++) { - this->at(i) = s(i); + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) *= e; + } + +} +/* -------------------------------------------------------------------------- */ + +template +template +void Map2d::operator/=(const T2 & e){ + + UInt sz = this->n[0]*this->n[1]; +#pragma omp parallel for + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) /= e; + } + +} +/* -------------------------------------------------------------------------- */ + +template +void Map2d::operator*=(const T & e){ + + UInt sz = this->n[0]*this->n[1]; +#pragma omp parallel for + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) *= e; + } +} + +/* -------------------------------------------------------------------------- */ + +template +void Map2d::operator/=(const T & e){ + + UInt sz = this->n[0]*this->n[1]; +#pragma omp parallel for + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) /= e; + } + +} + +/* -------------------------------------------------------------------------- */ + +template +template class surf, typename T2> +void Map2d::operator*=(const surf & e){ + + UInt sz = this->n[0]*this->n[1]; +#pragma omp parallel for + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) *= e.at(i); + } + +} +/* -------------------------------------------------------------------------- */ + +template +template class surf, typename T2> +void Map2d::operator/=(const surf & e){ + + UInt sz = this->n[0]*this->n[1]; +#pragma omp parallel for + for (unsigned int i = 0; i < sz; ++i) { + this->at(i) /= e.at(i); } } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #endif /* __MAP_2D_HH__ */ diff --git a/tests/SConscript b/tests/SConscript new file mode 100644 index 0000000..644f79d --- /dev/null +++ b/tests/SConscript @@ -0,0 +1,18 @@ +from os.path import join + +Import('main_env') + +test_files = Split(""" +run_tests.sh +test_hertzian.py +test_westergaard.py +test_fractal.py +test_fractal.verified +""") + +src_dir = "#/tests" +build_dir = 'build-' + main_env['build_type'] + '/tests' + +for file in test_files: + source = join(src_dir, file) + Command(file, source, Copy("$TARGET", "$SOURCE")) diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 89bad37..2aebe22 100755 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -1,32 +1,33 @@ #!/bin/bash n_tests=0 passed_tests=0 for test_file in test_*.py; do ((n_tests++)) echo "Testing $test_file" echo -n "------------ " python $test_file >/dev/null 2>&1 result=$? if [ $result -ne 0 ]; then echo "failure" else echo "success" ((passed_tests++)) fi done python test_fractal.py 2>/dev/null | tail -n 19 > test_fractal.out 2>&1 diff test_fractal.out test_fractal.verified >/dev/null 2>&1 result=$? if [ $result -ne 0 ]; then + echo "" echo "test_fractal.py diff failed" ((passed_tests--)) fi rm test_fractal.out failed_tests=$((n_tests - passed_tests)) echo "" echo "Test results :" echo " - passed $passed_tests/$n_tests" echo " - failed $failed_tests/$n_tests" diff --git a/tests/test_fftransform.py_no b/tests/test_fftransform.py_no index bd579d0..883d8e1 100644 --- a/tests/test_fftransform.py_no +++ b/tests/test_fftransform.py_no @@ -1,55 +1,58 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas 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 Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. If not, see . # ----------------------------------------------------------------------------- import sys -import tamaas as tm import numpy as np import matplotlib.pyplot as plt +from util_test import import_tamaas + def main(): + tm = import_tamaas() + size = 512 # Generate random rough surface SG = tm.SurfaceGeneratorFilterFFT() SG.getGridSize().assign(size) SG.getHurst().assign(0.8) SG.getRMS().assign(1.); SG.getQ0().assign(4); SG.getQ1().assign(4); SG.getQ2().assign(32); SG.getRandomSeed().assign(156); SG.Init() surface = SG.buildSurface() # Compute fourier transform through numpy numpy_spectral = np.fft.fft2(surface) # Compute fourier tranform through tamaas tamaas_spectral = np.zeros((size, size)) #transform = tm.FFTransformReal(size, surface, transform) return 0 if __name__ == '__main__': sys.exit(main()) diff --git a/tests/test_fractal.py b/tests/test_fractal.py index 750b7e1..993ae2f 100644 --- a/tests/test_fractal.py +++ b/tests/test_fractal.py @@ -1,100 +1,101 @@ #!/usr/bin/python # -*- coding: utf-8 -*- ##* # # @author Guillaume Anciaux # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas 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 Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. If not, see . # # ################################################################ -import sys -from tamaas import * +import sys, imp + +import tamaas as tm def main(): #generate surface - SG = SurfaceGeneratorFilterFFT() + SG = tm.SurfaceGeneratorFilterFFT() SG.getGridSize().assign(512) SG.getHurst().assign(0.8) SG.getRMS().assign(1.); SG.getQ0().assign(4); SG.getQ1().assign(4); SG.getQ2().assign(32); SG.getRandomSeed().assign(156); SG.Init() a = SG.buildSurface() #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 = SurfaceStatistics.computeStdev(SurfaceForPythonReal(a)); - stats.rms_slopes_geometric = SurfaceStatistics.computeRMSSlope(a); - stats.moments = SurfaceStatistics.computeMoments(a); + stats.rms_spectral = tm.SurfaceStatistics.computeSpectralStdev(a); + stats.rms_slopes_spectral = tm.SurfaceStatistics.computeSpectralRMSSlope(a); + stats.rms_geometric = tm.SurfaceStatistics.computeStdev(tm.SurfaceForPythonReal(a)); + stats.rms_slopes_geometric = tm.SurfaceStatistics.computeRMSSlope(a); + stats.moments = tm.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.m0prime = tm.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_m0 = tm.SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L); + stats.analytic_m2 = tm.SurfaceStatistics.computeAnalyticFractalMoment(2,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L); + stats.analytic_m4 = tm.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} """.format(**stats.__dict__) return 0 if __name__ == "__main__": sys.exit(main()) diff --git a/tests/test_hertzian.py b/tests/test_hertzian.py index 7164ca1..d2bec8e 100644 --- a/tests/test_hertzian.py +++ b/tests/test_hertzian.py @@ -1,121 +1,122 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas 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 Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. If not, see . # ----------------------------------------------------------------------------- import sys -import tamaas as tm import numpy as np import matplotlib.pyplot as plt +import tamaas as tm + 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.mean() 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 = 1024 curvature = 0.1 effective_modulus = 1. load = 0.0001 surface = constructHertzProfile(grid_size, curvature) bem = tm.BemPolonski(surface) bem.setEffectiveModulus(effective_modulus) bem.computeEquilibrium(1e-12, load) tractions = bem.getTractions() displacements = bem.getDisplacements() # 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-3 or pressure_error > 3e-3 or displacement_error > 1e-4: return 1 return 0 if __name__ == "__main__": sys.exit(main()) diff --git a/tests/test_westergaard.py b/tests/test_westergaard.py index e0a20fd..88d7d3e 100644 --- a/tests/test_westergaard.py +++ b/tests/test_westergaard.py @@ -1,76 +1,77 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas 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 Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. If not, see . # ----------------------------------------------------------------------------- import sys -import tamaas as tm import numpy as np import matplotlib.pyplot as plt +import tamaas as tm + def constructSinProfile(size, mode, amplitude): x = np.linspace(0, 1, size) y = np.linspace(0, 1, size) x, y = np.meshgrid(x, y) surface = amplitude * np.sin(2*np.pi*x*mode) return surface.copy() def main(): grid_size = 256 mode = 1 delta = 0.1 effective_modulus = 1. p_star = np.pi * effective_modulus * delta * mode # Full contact load load = 0.4 * p_star surface = constructSinProfile(grid_size, mode, delta) bem = tm.BemPolonski(surface) bem.setEffectiveModulus(effective_modulus) bem.computeEquilibrium(1e-12, load) tractions = bem.getTractions() displacements = bem.getDisplacements() # Testing contact area against Westergaard solution contact_area = tm.SurfaceStatistics.computeContactArea(tractions) westergaard_contact_size = 2. / (mode * np.pi) * np.arcsin(np.sqrt(load / p_star)) area_error = np.abs(contact_area - mode * westergaard_contact_size) / (mode * westergaard_contact_size) print "Area error: {}".format(area_error) # Testing displacements at full contact load = 1.1 * p_star bem.computeEquilibrium(1e-12, load) contact_area = tm.SurfaceStatistics.computeContactArea(tractions) print "Area at load > full contact : {}".format(contact_area) # Testing pressure distribution at full contact westergaard_pressure = constructSinProfile(grid_size, mode, p_star) tractions_amplitude = tractions.max() - tractions.mean() return 0 if __name__ == "__main__": sys.exit(main())