diff --git a/.gitignore b/.gitignore index 7f24725..03731e9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *~ .sconf_temp .sconsign.dblite *.log build* *.pyc +src/tamaas_info.cpp diff --git a/SConstruct b/SConstruct index fd4e12d..72bc43b 100644 --- a/SConstruct +++ b/SConstruct @@ -1,130 +1,139 @@ import os,sys +from version import write_info_file #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(ENV = os.environ) compiler_default = 'g++' if 'CXX' in os.environ: compiler_default = os.environ['CXX'] 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('CXX', 'Compiler', compiler_default) 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)) +# Save all options, not just those that differ from default +with open('build-setup.conf', 'w') as setup: + for key in vars.keys(): + setup.write("{} = '{}'\n".format(key, main_env[key])) +# instead of +#vars.Save('build-setup.conf', 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']) if 'CXXFLAGS' in os.environ: main_env.AppendUnique(CXXFLAGS = Split(os.environ['CXXFLAGS'])) main_env.AppendUnique(SHLINKFLAGS=['-fopenmp']) if timer_flag: main_env.AppendUnique(CPPDEFINES=['USING_TIMER']) cxxflags_dict = { "debug":Split("-g -O0 -DTAMAAS_DEBUG"), "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" 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'] main_env = main_env.Clone( tools = [fftw], FFTW_LIBRARY_WISH = ['main', 'omp'], FFTW_INCLUDE_DIR = fftw_include, FFTW_LIBRARY_DIR = fftw_library, ) -# save new values -vars.Save('build-setup.conf', main_env) +# writing information file +write_info_file("src/tamaas_info.cpp") + 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/examples/SConstruct b/examples/SConstruct index 0dd887e..c871980 100644 --- a/examples/SConstruct +++ b/examples/SConstruct @@ -1,12 +1,12 @@ import os env = Environment(CC = 'g++', CXXFLAGS = ['-std=c++11', '-g', '-fopenmp', '-O3'], TOOLS = ['default', tamaas], - TAMAAS_PATH = ['..'], + TAMAAS_PATH = '..', TAMAAS_BUILD_TYPE = 'debug', LIBS = ['gomp', 'Tamaas'], ENV = os.environ) env.Program('rough_surface.cc') diff --git a/examples/site_scons/site_init.py b/examples/site_scons/site_init.py index a95334d..61d748e 100644 --- a/examples/site_scons/site_init.py +++ b/examples/site_scons/site_init.py @@ -1,32 +1,34 @@ from SCons.Script import * from os.path import * +import os + def tamaas(env): """Sets correct environement variables for Tamaas. - 'TAMAAS_PATH' is list of potential tamaas repositories - 'TAMAAS_BUILD_TYPE' is the build type of tamaas (release, debug, profiling)""" if env.GetOption("clean"): return - build_type = "release" - if 'TAMAAS_BUILD_TYPE' in env: + try: build_type = env['TAMAAS_BUILD_TYPE'] - path = findPath(env['TAMAAS_PATH'], build_type) + except: + build_type = "release" - if not path: - print "Tamaas ({}) not found in {}".format(build_type, env['TAMAAS_PATH']) + if 'TAMAAS_PATH' not in env: + print "Please set the 'TAMAAS_PATH' variable in your scons environement" Exit(1) - env.AppendUnique(CPPPATH = [abspath(join(path, 'src'))]) - env.AppendUnique(LIBPATH = [abspath(join(path, 'build-{}/src'.format(build_type)))]) - env.AppendUnique(RPATH = [abspath(join(path, 'build-{}/src'.format(build_type)))]) + tm_path = env['TAMAAS_PATH'] + include_dir = tm_path + '/src' + lib_dir = join(tm_path, 'build-{}/src'.format(build_type)) -def findPath(search, build_type): - trial = Environment() - conf = Configure(trial) - for path in search: - trial['CPPPATH'] = join(path, 'src') - trial['LIBPATH'] = join(path, 'build-{}/src'.format(build_type)) - trial['RPATH'] = join(path, 'build-{}/src'.format(build_type)) + env.AppendUnique(CPPPATH = [abspath(include_dir)]) + env.AppendUnique(LIBPATH = [abspath(lib_dir)]) + env.AppendUnique(RPATH = [abspath(lib_dir)]) - if conf.CheckLib("Tamaas", language="CXX"): - return path + conf = Configure(env) + if not conf.CheckLibWithHeader("Tamaas", "tamaas.hh", language="CXX"): + raise SCons.Errors.StopError( + "Tamaas ({}) not found in {}".format(build_type, abspath(tm_path))) + env['TAMAAS_PATH'] = abspath(tm_path) + env = conf.Finish() diff --git a/python/SConscript b/python/SConscript index beb48c9..234716a 100644 --- a/python/SConscript +++ b/python/SConscript @@ -1,36 +1,38 @@ Import('main_env') import distutils.sysconfig +from os.path import abspath from numpy.distutils.misc_util import get_numpy_include_dirs as numpy_dirs from distutils.version import StrictVersion env_swig = main_env.Clone( tools = ['swig'], SWIG = 'swig', SHLIBPREFIX = '', ) if env_swig['SWIGVERSION'] is None: print "Could not find swig version" print env_swig.Dictionary() print env_swig.Environment Exit(1) if StrictVersion(env_swig['SWIGVERSION']) < StrictVersion("3.0"): print "Swig version {} is not supported by tamaas".format(env_swig['SWIGVERSION']) Exit(1) env_swig.AppendUnique(CPPPATH=[distutils.sysconfig.get_python_inc(), numpy_dirs()]) env_swig.AppendUnique(SWIGPATH=['$CPPPATH']) env_swig.AppendUnique(SWIGFLAGS=['-python', '-c++']) verbose = main_env['verbose'] if not verbose: env_swig.AppendUnique(SWIGFLAGS=['-w325', '-w315']) env_swig.AppendUnique(CXXFLAGS=['-Wno-maybe-uninitialized']) # for some warning in wrapper code env_swig.SharedLibrary( target = '_tamaas.so', source = ['tamaas.i'], - LIBS = ['Tamaas'] + LIBS = ['Tamaas'], + RPATH = [abspath('../src')] ) diff --git a/python/bem.i b/python/bem.i index eb21902..7b6e0da 100644 --- a/python/bem.i +++ b/python/bem.i @@ -1,58 +1,59 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Wrapping bem capabilities /* -------------------------------------------------------------------------- */ %include "bem_interface.hh" %include "bem_fft_base.hh" %include "bem_kato.hh" %include "bem_polonski.hh" %include "bem_gigi.hh" %include "bem_penalty.hh" %include "bem_uzawa.hh" %include "bem_gigipol.hh" %include "bem_grid.hh" %include "bem_grid_kato.hh" %include "bem_grid_polonski.hh" %include "bem_grid_teboulle.hh" +%include "bem_grid_condat.hh" %include "functional.hh" %include "meta_functional.hh" %include "elastic_energy_functional.hh" %include "complimentary_term_functional.hh" %include "exponential_adhesion_functional.hh" %include "squared_exponential_adhesion_functional.hh" %include "maugis_adhesion_functional.hh" %include "fftransform.hh" %include "fftransform_fftw.hh" namespace tamaas { %template(FFTransformReal) FFTransform; %template(FFTransformFFTWReal) FFTransformFFTW; } diff --git a/python/surface.i b/python/surface.i index 683cefb..d8d8b2f 100644 --- a/python/surface.i +++ b/python/surface.i @@ -1,380 +1,381 @@ /** * * @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 "map_2d.hh" #include "surface.hh" #include "types.hh" #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include "surface_generator_filter_fft.hh" #include "smart_pointer.hh" #include "surface_statistics.hh" #include "bem_kato.hh" #include "bem_polonski.hh" #include "bem_penalty.hh" #include "bem_gigi.hh" #include "bem_gigipol.hh" #include "bem_uzawa.hh" #include "bem_fft_base.hh" #include "bem_grid.hh" #include "bem_grid_polonski.hh" #include "bem_grid_kato.hh" #include "bem_grid_teboulle.hh" +#include "bem_grid_condat.hh" #include "contact_area.hh" #include "contact_cluster_collection.hh" #include "functional.hh" #include "meta_functional.hh" #include "elastic_energy_functional.hh" #include "complimentary_term_functional.hh" #include "exponential_adhesion_functional.hh" #include "squared_exponential_adhesion_functional.hh" #include "maugis_adhesion_functional.hh" #include "surface_timer.hh" #include "fftransform_fftw.hh" %} %include "cpointer.i" %include "typemaps.i" %include "numpy.i" %init %{ import_array(); %} /* -------------------------------------------------------------------------- */ /* Wrapping returned pointers /* -------------------------------------------------------------------------- */ //%include "typemaps.i" %include "std_string.i" %include "std_vector.i" %include "smart_pointer.hh" %template(VectorReal) std::vector; //%pointer_functions(int, intp); //%pointer_functions(long, longp); //%pointer_functions(unsigned int, uintp); //%pointer_functions(double, doublep); //%pointer_functions(float, floatp); //%pointer_functions(std::complex, comlpexp); %define %my_pointer_class(NAME,TYPE) %template(SmartPointer##NAME) SmartPointer< TYPE >; %typemap(out) TYPE & { #define TOREPLACE SWIGTYPE_p_SmartPointerT_##NAME##_t SmartPointer * ptr = new SmartPointer($1); $result = SWIG_NewPointerObj(SWIG_as_voidptr(ptr), TOREPLACE, 0 | 0 ); #undef TOREPLACE } %enddef %my_pointer_class(int,int) %my_pointer_class(double,double) %my_pointer_class(long,long) %my_pointer_class(vectorReal,std::vector) %apply (unsigned int IN_ARRAY1[ANY]) {(tamaas::UInt n[ANY])}; /* -------------------------------------------------------------------------- */ /* Wrapping Surface class */ /* -------------------------------------------------------------------------- */ %include "tamaas.hh" %ignore tamaas::Grid::iterator::operator++; %ignore *::operator=; %include "grid.hh" %include "grid_hermitian.hh" namespace tamaas { %template (Grid1dReal) Grid; %template (Grid2dReal) Grid; %template (Grid2dComplex) Grid; %template (Grid2dInt) Grid; %template (Grid2dUInt) Grid; %template (GridHermitian2dReal) GridHermitian; } %include "types.hh" %include "grid_hermitian.hh" %include "map_2d.hh" namespace tamaas { %template(Map2dComplex) Map2d; %template(Map2dInt) Map2d; %template(Map2dReal) Map2d; } %include "map_2d_square.hh" namespace tamaas { %template(Map2dSquareComplex) Map2dSquare; %template(Map2dSquareInt) Map2dSquare; %template(Map2dSquareReal) Map2dSquare; } %include "surface.hh" %include "surface_complex.hh" namespace tamaas { %template(SurfaceReal) Surface; %template(SurfaceInt) Surface; %template(SurfaceRealComplex) SurfaceComplex; } %inline %{ __BEGIN_TAMAAS__ template void sanityCheck(PyObject * input) { if (PyArray_TYPE((PyArrayObject*)input) != required_type) throw Exception("Sanity check: python object stores wrong type"); } __END_TAMAAS__ %} %inline %{ __BEGIN_TAMAAS__ template class GridForPython : public Grid { public: GridForPython(T * wrapped_memory, UInt n[d], Real L[d], UInt nb_components): Grid() { - this->__init__(wrapped_memory, n, L); + this->__init__(wrapped_memory, n, nb_components); } void __init__(T * w, UInt n[d], UInt nb_components) { memcpy(this->n, n, d * sizeof(UInt)); memset(this->L, 1, d * sizeof(Real)); this->nb_components = nb_components; this->data.wrapMemory(w, this->computeSize()); } GridForPython(PyObject * input): Grid() { if (!PyArray_Check(input)) { PyObject* obj_type = PyObject_Type(input); std::string type = PyString_AsString(obj_type); throw Exception("GridForPython: incompatible input which is not a numpy: " + type); } UInt _n = PyArray_NDIM((PyArrayObject*)input); // numpy number of dimensions npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input); // numpy dimensions UInt grid_dims[d] = {0}; // tamaas dims if (_n >= d) { for (UInt i = 0 ; i < d ; ++i) grid_dims[i] = ndims[i]; } else { throw Exception("GridForPython: incompatible input size for numpy array: " + _n); } if (_n > d) this->nb_components = ndims[d]; else this->nb_components = 1; PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input); if (iter == NULL) throw; this->__init__(( T *)(iter->dataptr), grid_dims, this->nb_components); Py_DECREF(iter); } virtual ~GridForPython() {} }; __END_TAMAAS__ %} %inline %{ __BEGIN_TAMAAS__ template class SurfaceForPython : public Surface{ public: SurfaceForPython(T * wrapped_memory, UInt n, Real L) : Surface(0,L){ this->__init__(wrapped_memory,n,L); }; void __init__(T * wrapped_memory, UInt n, Real L){ this->n[0] = n; this->n[1] = n; this->data.wrapMemory(wrapped_memory,n*n); } SurfaceForPython(PyObject * input, Real L = 1.) : Surface(0,L){ if (!PyArray_Check(input)){ PyObject* obj_type = PyObject_Type(input); std::string type = PyString_AsString(obj_type); throw Exception("SurfaceForPython: incompatible input which is not a numpy: " + type); } UInt _n = PyArray_NDIM((PyArrayObject*)input); if (_n != 2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimension " << _n); npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input); UInt sz = ndims[0]; UInt sz2 = ndims[1]; if (sz != sz2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimensions " << sz << " " << sz2); PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input); if (iter == NULL) throw; this->__init__(( T *)(iter->dataptr),sz,L); Py_DECREF(iter); }; ~SurfaceForPython(){ }; void resize(UInt new_size){ if (this->size() != new_size) throw Exception("SurfaceForPython: cannot resize a temporary vector"); } }; __END_TAMAAS__ %} %define %grid_typemaps(GRID_CLASS, DATA_TYPE, DATA_TYPECODE, DIM) %typemap(out, fragment="NumPy_Fragments") (GRID_CLASS &) { using namespace tamaas; UInt additional_dim = ($1->getNbComponents() == 1)? 0:1; UInt dim = DIM+additional_dim; npy_intp dims[dim]; for (UInt i = 0 ; i < DIM ; i++) dims[i] = (npy_intp)$1->sizes()[i]; if (additional_dim) dims[dim-1] = $1->getNbComponents(); DATA_TYPE * data = const_cast($1->getInternalData()); PyObject* obj = PyArray_SimpleNewFromData(dim, dims, DATA_TYPECODE, &data[0]); PyArrayObject* array = (PyArrayObject*) obj; if (!array) SWIG_fail; $result = SWIG_Python_AppendOutput($result, obj); } %typemap(in) tamaas::Grid & { using namespace tamaas; if (!PyArray_Check($input)) { PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input)); const char* s = PyString_AsString(objectsRepresentation); std::string __type = s; if ("" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else if ("" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")"); } else { $1 = new GridForPython< DATA_TYPE, DIM >($input); sanityCheck< DATA_TYPECODE >($input); } } %inline %{ namespace tamaas { GRID_CLASS & convertGrid(GRID_CLASS & s){ return s; } } %} %enddef %define %surface_typemaps(SURFACE_CLASS, DATA_TYPE, DATA_TYPECODE) %typemap(out, fragment="NumPy_Fragments") (SURFACE_CLASS &) { using namespace tamaas; UInt nz = $1->size(); npy_intp dims[2] = {nz,nz}; DATA_TYPE * data = const_cast($1->getInternalData()); PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, &data[0]); PyArrayObject* array = (PyArrayObject*) obj; if (!array) SWIG_fail; $result = SWIG_Python_AppendOutput($result, obj); } %typemap(in) tamaas::Surface< DATA_TYPE > & { using namespace tamaas; if (!PyArray_Check($input)){ PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input)); const char* s = PyString_AsString(objectsRepresentation); std::string __type = s; if ("" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else if ("" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else if ("" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")"); } else{ $1 = new SurfaceForPython< DATA_TYPE >($input); sanityCheck< DATA_TYPECODE >($input); } } %inline %{ namespace tamaas { SURFACE_CLASS & convertSurface(SURFACE_CLASS & s){ return s; } } %} %enddef %surface_typemaps(tamaas::Surface, Real, NPY_DOUBLE) %surface_typemaps(tamaas::SurfaceComplex, std::complex, NPY_COMPLEX128 ) %surface_typemaps(tamaas::Surface, unsigned long, NPY_UINT ) %surface_typemaps(tamaas::Surface, UInt, NPY_UINT ) %surface_typemaps(tamaas::Surface, int, NPY_INT ) %grid_typemaps(tamaas::Grid1dReal, Real, NPY_DOUBLE, 1) %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2) %grid_typemaps(tamaas::GridHermitian2dReal, tamaas::complex, NPY_COMPLEX128, 2) %grid_typemaps(tamaas::Grid2dUInt, UInt, NPY_UINT, 2) %grid_typemaps(tamaas::SurfaceComplex, tamaas::complex, NPY_COMPLEX128, 2) namespace tamaas { %template(SurfaceForPythonReal) SurfaceForPython< Real >; } %include "surface_timer.hh" diff --git a/src/SConscript b/src/SConscript index 6d197cd..50842c9 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,78 +1,80 @@ Import('main_env') import os env = main_env # 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 = """ +tamaas_info.cpp bem_kato.cpp bem_polonski.cpp bem_gigi.cpp bem_gigipol.cpp bem_penalty.cpp bem_uzawa.cpp bem_fft_base.cpp functional.cpp meta_functional.cpp elastic_energy_functional.cpp exponential_adhesion_functional.cpp squared_exponential_adhesion_functional.cpp maugis_adhesion_functional.cpp complimentary_term_functional.cpp fftransform.cpp fftransform_fftw.cpp fft_plan_manager.cpp grid.cpp grid_hermitian.cpp types.cpp bem_grid.cpp bem_grid_polonski.cpp bem_grid_kato.cpp bem_grid_teboulle.cpp +bem_grid_condat.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 f659097..88e3e96 100644 --- a/src/bem_fft_base.cpp +++ b/src/bem_fft_base.cpp @@ -1,270 +1,273 @@ /** * * @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 "meta_functional.hh" #include "elastic_energy_functional.hh" #include "fft_plan_manager.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ BemFFTBase::BemFFTBase(Surface & p): BemInterface(p), functional(new MetaFunctional(*this)), surface_tractions(p.size(),p.getL()), surface_displacements(p.size(),p.getL()), 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("")), tmp_coeff(1./M_PI){ // this->setNumberOfThreads(this->nthreads); max_iterations = 5000; dump_freq = 100; this->functional->addFunctionalTerm(new ElasticEnergyFunctional(*this)); std::cerr << this << " is setting up the surfaces"; std::cerr << &p << " has size " << p.size() << std::endl; } /* -------------------------------------------------------------------------- */ BemFFTBase::~BemFFTBase() { FFTPlanManager::get().clean(); } /* -------------------------------------------------------------------------- */ 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; } /* -------------------------------------------------------------------------- */ void BemFFTBase::setEffectiveModulus(Real e_star){ this->e_star = e_star; this->computeSpectralInfluenceOverDisplacement(); } /* -------------------------------------------------------------------------- */ - +Real BemFFTBase::getEffectiveModulus(){ + return this->e_star; +} +/* -------------------------------------------------------------------------- */ 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::computeDisplacementsFromTractions(){ STARTTIMER("computeDisplacementsFromTractions"); this->applyInfluenceFunctions(surface_tractions,surface_displacements); STOPTIMER("computeDisplacementFromTractions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeTractionsFromDisplacements(){ STARTTIMER("computeDisplacementsFromTractions"); this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions); STOPTIMER("computeDisplacementFromTractions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::applyInfluenceFunctions(Surface & input, Surface & output){ STARTTIMER("applyInfluenceFunctions"); FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); plan1.forwardTransform(); surface_spectral_output *= surface_spectral_influence_disp; FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); plan2.backwardTransform(); STOPTIMER("applyInfluenceFunctions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::applyInverseInfluenceFunctions(Surface & input, Surface & output){ STARTTIMER("applyInfluenceFunctions"); FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); plan1.forwardTransform(); surface_spectral_output /= surface_spectral_influence_disp; surface_spectral_output(0) = 0; FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); plan2.backwardTransform(); 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); } } /* -------------------------------------------------------------------------- */ 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::addFunctional(Functional *new_funct) { this->functional->addFunctionalTerm(new_funct); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_fft_base.hh b/src/bem_fft_base.hh index 7d9c5a3..ed20926 100644 --- a/src/bem_fft_base.hh +++ b/src/bem_fft_base.hh @@ -1,152 +1,158 @@ /** * * @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 "meta_functional.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 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; //! get the computed true displacement const Surface & getTrueDisplacements(){return true_displacements;}; //! get the computed gap const Surface & getGap(){return gap;}; //! set the elastic effective modulus void setEffectiveModulus(Real e_star); + //! set the elastic effective modulus + Real getEffectiveModulus(); //! set the maximal number of iterations void setMaxIterations(UInt max_iter){this->max_iterations = max_iter;}; + //! get the maximal number of iterations + UInt getMaxIterations(){return this->max_iterations;}; //! set dump freq - void setDumpFreq(Real dump_freq){this->dump_freq = dump_freq;}; + void setDumpFreq(UInt dump_freq){this->dump_freq = dump_freq;}; + //! get dump freq + UInt getDumpFreq(){return this->dump_freq;}; //! set functional void addFunctional(::tamaas::Functional * new_funct); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! Functional for optimization MetaFunctional * functional; //! surface tractions Surface surface_tractions; //! surface displacement Surface surface_displacements; //! 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; //! 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_grid.cpp b/src/bem_grid.cpp index a5581d8..2c25c9b 100644 --- a/src/bem_grid.cpp +++ b/src/bem_grid.cpp @@ -1,423 +1,432 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "bem_grid.hh" #include "grid.hh" #include "types.hh" #include "fft_plan_manager.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ BemGrid::BemGrid(Surface & surface): E(1.), nu(0.3), surface(surface), surface_normals(surface.sizes(), surface.lengths(), 3), tractions(surface.sizes(), surface.lengths(), 3), displacements(surface.sizes(), surface.lengths(), 3), gaps(surface.sizes(), surface.lengths(), 3), true_tractions(NULL), true_displacements(NULL), dump_freq(100) { UInt hermitian_sizes[2] = {surface.sizes()[0], surface.sizes()[1]/2+1}; influence.setNbComponents(9); spectral_tractions.setNbComponents(3); spectral_displacements.setNbComponents(3); influence.resize(hermitian_sizes); spectral_tractions.resize(hermitian_sizes); spectral_displacements.resize(hermitian_sizes); //computeSurfaceNormals(); } /* -------------------------------------------------------------------------- */ BemGrid::~BemGrid() { FFTPlanManager::get().clean(); delete true_tractions; delete true_displacements; } /* -------------------------------------------------------------------------- */ const Grid2dReal & BemGrid::getTrueTractions() const { if (true_tractions == NULL) TAMAAS_EXCEPTION("True tractions were not allocated"); return *true_tractions; } const Grid2dReal & BemGrid::getTrueDisplacements() const { if (true_displacements == NULL) TAMAAS_EXCEPTION("True displacements were not allocated"); return *true_displacements; } /* -------------------------------------------------------------------------- */ /// This function computes surface slopes in fourier space and fills the normals grid void BemGrid::computeSurfaceNormals() { SurfaceComplex spectral_surface(surface.size(), surface.getL()); GridHermitian spectral_slopes(spectral_surface.sizes(), spectral_surface.lengths(), 3); FFTransform & plan1 = FFTPlanManager::get().createPlan(surface, spectral_surface); FFTransform & plan2 = FFTPlanManager::get().createPlan(surface_normals, spectral_slopes); plan1.forwardTransform(); // Loops should be with signed integers const Int * n = reinterpret_cast(spectral_surface.sizes()); constexpr complex I(0, 1); #pragma omp parallel { #pragma omp for collapse(2) for (Int i = 0 ; i <= n[0] / 2 ; i++) { for (Int j = 0 ; j < n[1] ; j++) { Real q1 = 2*M_PI*j; Real q2 = 2*M_PI*i; spectral_slopes(i, j, 0) = -q1*I*spectral_surface(i, j); spectral_slopes(i, j, 1) = -q2*I*spectral_surface(i, j); spectral_slopes(i, j, 2) = 0; } } #pragma omp for collapse(2) for (Int i = n[0]/2 + 1 ; i < n[0] ; i++) { for (Int j = 0 ; j < n[1] ; j++) { Real q1 = 2*M_PI*j; Real q2 = 2*M_PI*(i-n[0]); spectral_slopes(i, j, 0) = -q1*I*spectral_surface(i, j); spectral_slopes(i, j, 1) = -q2*I*spectral_surface(i, j); spectral_slopes(i, j, 2) = 0; } } } // end #pragma omp parallel // Dirac spectral_slopes(0, 0, 2) = n[0]*n[0]; plan2.backwardTransform(); FFTPlanManager::get().destroyPlan(surface, spectral_surface); FFTPlanManager::get().destroyPlan(surface_normals, spectral_slopes); // Make surface normal VectorProxy sn(nullptr, 3); const UInt N = surface_normals.getNbPoints(); #pragma omp parallel for firstprivate(sn) for (UInt i = 0 ; i < N ; i++) { sn.setPointer(surface_normals.getInternalData() + i * sn.dataSize()); sn /= sn.l2norm(); } } /* -------------------------------------------------------------------------- */ void BemGrid::computeDisplacementsFromTractions() { applyInfluenceFunctions(tractions, displacements); } void BemGrid::computeTractionsFromDisplacements() { applyInverseInfluenceFunctions(tractions, displacements); } /* -------------------------------------------------------------------------- */ void BemGrid::computeGaps(Grid &disp) { const UInt * n = surface.sizes(); #pragma omp parallel for collapse(2) for (UInt i = 0 ; i < n[0] ; i++) { for (UInt j = 0 ; j < n[1] ; j++) { gaps(i, j, 0) = disp(i, j, 0); gaps(i, j, 1) = disp(i, j, 1); gaps(i, j, 2) = disp(i, j, 2) - surface(i, j); } } } /* -------------------------------------------------------------------------- */ void BemGrid::applyInfluenceFunctions(Grid & input, Grid & output) { TAMAAS_ASSERT(input.dataSize() == output.dataSize(), "input and output have different size"); FFTransform & plan_1 = FFTPlanManager::get().createPlan(input, spectral_tractions); FFTransform & plan_2 = FFTPlanManager::get().createPlan(output, spectral_displacements); plan_1.forwardTransform(); complex *F = influence.getInternalData(), *p = spectral_tractions.getInternalData(), *u = spectral_displacements.getInternalData(); VectorProxy spectral_p(nullptr, 3), spectral_u(nullptr, 3); MatrixProxy spectral_F(nullptr, 3, 3); spectral_displacements = 0; const UInt N = influence.getNbPoints(); // OpenMP not very cool with this: best to have a POD iteration variable #pragma omp parallel for firstprivate(spectral_p, spectral_u, spectral_F) for (UInt i = 0 ; i < N ; ++i) { spectral_p.setPointer(p + i * spectral_p.dataSize()); spectral_u.setPointer(u + i * spectral_u.dataSize()); spectral_F.setPointer(F + i * spectral_F.dataSize()); spectral_u.mul(spectral_F, spectral_p); } plan_2.backwardTransform(); } /* -------------------------------------------------------------------------- */ void BemGrid::applyInverseInfluenceFunctions(Grid & __attribute__((unused)) input, Grid & __attribute__((unused)) output) { TAMAAS_EXCEPTION("applyInverseInfluenceFunctions not yet implemented"); } /* -------------------------------------------------------------------------- */ // warning: this function will suck your soul out // loops must be written with signed integers void BemGrid::computeInfluence() { // Iterator over points of surface MatrixProxy it(influence.getInternalData(), 3, 3); // Imaginary unit constexpr complex I(0, 1); const Int size = this->surface.size(); this->influence = 0; // Influence functions matrix, ask Lucas #define INFLUENCE(q_1, q_2, q, E, nu) { \ 2./(E*q*q*q)*(-nu*nu*q_1*q_1 - nu*q_2*q_2 + q*q), \ -q_1*q_2/(E*q*q*q)*(1+nu)*2*nu, \ I*q_1/(E*q*q)*(1+nu)*(1-2*nu), \ -q_1*q_2/(E*q*q*q)*(1+nu)*2*nu, \ 2./(E*q*q*q)*(-nu*nu*q_2*q_2 - nu*q_1*q_1 + q*q), \ I*q_2/(E*q*q)*(1+nu)*(1-2*nu), \ -I*q_1/(E*q*q)*(1+nu)*(1-2*nu), \ -I*q_2/(E*q*q)*(1+nu)*(1-2*nu), \ 2./(E*q)*(1-nu*nu)} // Fill influence for single point #define FILL_INFLUENCE(k, l) do { \ Real q_1 = 2*M_PI*(l); \ Real q_2 = 2*M_PI*(k); \ Real q = std::sqrt(q_1*q_1+q_2*q_2); \ complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \ MatrixProxy inf_mat(influence_1, 3, 3); \ it.setPointer(&influence(i, j, 0)); \ it = inf_mat; \ } while(0) // Fill influence points for "symetric" areas #define FILL_SYM_INFLUENCE(k, l) do { \ Real q_1 = 2*M_PI*(l); \ Real q_2 = 2*M_PI*(k); \ Real q = std::sqrt(q_1*q_1+q_2*q_2); \ complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \ complex influence_2[] = INFLUENCE(q_2, q_1, q, E, nu); \ MatrixProxy inf_mat(influence_1, 3, 3); \ it.setPointer(&influence(i, j, 0)); \ it = inf_mat; \ inf_mat.setPointer(influence_2); \ it.setPointer(&influence(j, i, 0)); \ it = inf_mat; \ } while(0) #pragma omp parallel firstprivate(it) // all the following loops are independent { // Loops over parts of surface /* [ ][ ][ ][ ] [ ][x][ ][ ] [ ][ ][ ][ ] [ ][ ][ ][o] */ #pragma omp for schedule(dynamic) // can't collapse triangluar loops for (Int i = 1 ; i < size/2 ; i++) { for (Int j = i ; j < size/2 ; j++) { FILL_SYM_INFLUENCE(i, j); } } /* [ ][ ][ ][ ] [ ][ ][ ][o] [ ][ ][ ][ ] [ ][x][ ][ ] */ #pragma omp for collapse(2) for (Int i = size/2+1 ; i < size ; i++) { for (Int j = 1 ; j < size/2 ; j++) { FILL_INFLUENCE(i-size, j); } } /* [ ][x][x][o] [x][ ][x][ ] [x][x][x][ ] [ ][ ][ ][ ] */ #pragma omp for for (Int i = 1 ; i <= size/2 ; i++) { { Int j = size/2; FILL_SYM_INFLUENCE(i, j); } { Int j = 0; FILL_SYM_INFLUENCE(i, j); } } /* [ ][ ][ ][ ] [ ][ ][ ][ ] [ ][ ][ ][o] [x][ ][x][ ] */ #pragma omp for for (Int i = size/2+1 ; i < size ; i++) { { Int j = size/2; FILL_INFLUENCE(i-size, j); } { Int j = 0; FILL_INFLUENCE(i-size, j); } } } // End #pragma omp parallel /* State of surface now: [0][x][x][o] [x][x][x][o] [x][x][x][o] [x][x][x][o] Legend: x => explicitely filled o => hermitian symetry 0 => actually 0 */ computeInfluenceLipschitz(); #undef INFLUENCE #undef FILL_INFLUENCE #undef FILL_SYM_INFLUENCE } /* -------------------------------------------------------------------------- */ /// Computing L_2 norm of influence functions void BemGrid::computeInfluenceLipschitz() { const UInt N = influence.getNbPoints(); VectorProxy m(nullptr, 9); Real norm = 0.; #pragma omp parallel for firstprivate(m) reduction(+:norm) for (UInt i = 0 ; i < N ; i++) { m.setPointer(influence.getInternalData() + i * m.dataSize()); norm += std::pow(std::norm(m.l2norm()), 2); } lipschitz = std::sqrt(norm); } /* -------------------------------------------------------------------------- */ /* Friction related computations */ /* -------------------------------------------------------------------------- */ void BemGrid::projectOnFrictionCone(Grid & field, Real mu) { VectorProxy p_proxy(nullptr, 3); const UInt N = field.getNbPoints(); #pragma omp parallel for firstprivate(p_proxy) for (UInt i = 0 ; i < N ; i++) { p_proxy.setPointer(field.getInternalData() + i * p_proxy.dataSize()); projectOnFrictionCone(p_proxy, mu); } } /* -------------------------------------------------------------------------- */ void BemGrid::enforceMeanValue(Grid & field, const Grid & target) { Real red_0 = 0, red_1 = 0, red_2 = 0; const UInt * n = field.sizes(); #pragma omp parallel for collapse(2) reduction(+:red_0, red_1, red_2) for (UInt i = 0 ; i < n[0] ; i++) { for (UInt j = 0 ; j < n[1] ; j++) { red_0 += field(i, j, 0); red_1 += field(i, j, 1); red_2 += field(i, j, 2); } } red_0 /= n[0]*n[1]; red_1 /= n[0]*n[1]; red_2 /= n[0]*n[1]; Real p1_data[3] = {red_0, red_1, red_2}; // Loop for pressure shifting VectorProxy p1(p1_data, 3); VectorProxy p0((Real*)target.getInternalData(), 3); VectorProxy p(nullptr, 3); const UInt N = field.getNbPoints(); #pragma omp parallel for firstprivate(p0, p1) for (UInt i = 0 ; i < N ; i++) { p.setPointer(field.getInternalData() + i * p.dataSize()); p -= p1; p += p0; } } /* -------------------------------------------------------------------------- */ -Real BemGrid::computeMeanValueError(Grid & field, const Grid & target) { +void BemGrid::computeMeanValueError(Grid & field, + const Grid & target, + Grid & error_vec) { Real p0 = 0, p1 = 0, p2 = 0; const UInt * n = field.sizes(); #pragma omp parallel for collapse(2) reduction(+:p0, p1, p2) for (UInt i = 0 ; i < n[0] ; i++) { for (UInt j = 0 ; j < n[1] ; j++) { p0 += field(i, j, 0); p1 += field(i, j, 1); p2 += field(i, j, 2); } } p0 /= n[0]*n[1]; p1 /= n[0]*n[1]; p2 /= n[0]*n[1]; + Real error[3] = {p0 - target(0), p1 - target(1), p2 - target(2)}; + VectorProxy e(error, 3); + error_vec = e; +} + +/* -------------------------------------------------------------------------- */ + +Real BemGrid::computeMeanValueError(Grid & field, const Grid & target) { + Real e[3] = {0.}; + VectorProxy error(e, 3); + computeMeanValueError(field, target, error); Real norm = target(0)*target(0) + target(1)*target(1) + target(2)*target(2); - Real error[3] = {p0 - target(0), p1 - target(1), p2 - target(2)}; - VectorProxy e(error, 3); - Real er = e.l2norm() / std::sqrt(norm); - //std::cout << std::scientific << "Eq error: " << er << std::fixed << std::endl; - return er; + return error.l2norm() / std::sqrt(norm); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_grid.hh b/src/bem_grid.hh index 56ebff1..58f9818 100644 --- a/src/bem_grid.hh +++ b/src/bem_grid.hh @@ -1,180 +1,184 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __BEM_GRID_HH__ #define __BEM_GRID_HH__ /* -------------------------------------------------------------------------- */ #include "surface.hh" #include "grid.hh" #include "types.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ // Necessary typedefs for swig typedef Grid Grid2dReal; typedef Grid Grid2dComplex; typedef Grid Grid2dUInt; typedef GridHermitian GridHermitian2dReal; /** * @brief abstract class for BEMs that use multi component fields */ class BemGrid { public: /// Constructor BemGrid(Surface & surface); /// Destructor virtual ~BemGrid(); public: /* -------------------------------------------------------------------------- */ /* Actual computations */ /* -------------------------------------------------------------------------- */ /// Compute normal vectors to surface void computeSurfaceNormals(); /// Compute influence coefficients void computeInfluence(); /// Compute a Lipschitz constant of gradient void computeInfluenceLipschitz(); /// Apply influence fonctions (p -> u) void applyInfluenceFunctions(Grid & input, Grid & output); /// Apply inverse influence fonctions (u -> p) void applyInverseInfluenceFunctions(Grid & input, Grid & output); /// Compute true tractions virtual void computeTrueTractions() = 0; /// Compute true displacements virtual void computeTrueDisplacements() = 0; /// Compute equilibrium virtual void computeEquilibrium(Real tol, const Grid & p_target) = 0; /// Compute displacements void computeDisplacementsFromTractions(); /// Compute tractions void computeTractionsFromDisplacements(); /// Compute gaps void computeGaps(Grid & disp); /* -------------------------------------------------------------------------- */ /* Field related computations */ /* -------------------------------------------------------------------------- */ /// Project point friction cone static inline void projectOnFrictionCone(VectorProxy & p, Real mu); /// Project field friction cone static void projectOnFrictionCone(Grid & field, Real mu); /// Project on mean value space static void enforceMeanValue(Grid & field, const Grid & target); /// Compute error relative to target mean value static Real computeMeanValueError(Grid & field, const Grid & target); + /// Compute error vector + static void computeMeanValueError(Grid & field, + const Grid & target, + Grid & error_vec); /* -------------------------------------------------------------------------- */ /* Getters */ /* -------------------------------------------------------------------------- */ /// Get influence const GridHermitian2dReal & getInfluence() const {return influence;} /// Get normals const Grid2dReal & getSurfaceNormals() const {return surface_normals;} /// Get tractions const Grid2dReal & getTractions() const {return tractions;} /// Get displacements const Grid2dReal & getDisplacements() const {return displacements;} /// Get true tractions const Grid2dReal & getTrueTractions() const; /// Get true displacements const Grid2dReal & getTrueDisplacements() const; /// Get gaps const Grid2dReal & getGaps() const {return gaps;} /* -------------------------------------------------------------------------- */ /* Setters */ /* -------------------------------------------------------------------------- */ /// Elasticity parameters void setElasticity(Real E_, Real nu_) {E = E_; nu = nu_;} /// Friction coefficient void setMu(Real mu_) {mu = mu_;} /// Max iterations void setMaxIterations(UInt n) {max_iter = n;} /// Dump frequency void setDumpFrequency(UInt n) {dump_freq = n;} protected: /// Elasticity constants Real E, nu; /// Friction coefficient Real mu; /// Lipschitz constant of functional gradient Real lipschitz; /// Max iterations UInt max_iter; /// Surface Surface & surface; /// Normal vectors Grid surface_normals; /// Influcence matrices GridHermitian influence; /// Tractions Grid tractions; /// Displacements Grid displacements; /// Gap Grid gaps; /// True tractions Grid * true_tractions; /// True displacements Grid * true_displacements; /// Spectral tractions GridHermitian spectral_tractions; /// Spectral displacements GridHermitian spectral_displacements; /// Dump frequency UInt dump_freq; }; /* -------------------------------------------------------------------------- */ /* Inline definitions */ /* -------------------------------------------------------------------------- */ inline void BemGrid::projectOnFrictionCone(VectorProxy & p, Real mu) { Real pt = std::sqrt(p(0)*p(0)+p(1)*p(1)); Real pn = p(2); Real x = (pn + mu*pt)/(1+mu*mu); if (pt <= mu * pn && pn >= 0) // inside friction cone return; else if (mu * pt <= std::abs(pn) && pn < 0) // inside "project to origin" cone (undefined normal) //else if (pn < 0) // project to origin all points with pn < 0 p = 0; else { // projectable on cone p(0) *= mu / pt * x; p(1) *= mu / pt * x; p(2) = x; } } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #endif // __BEM_GRID_HH__ diff --git a/src/bem_grid_condat.cpp b/src/bem_grid_condat.cpp new file mode 100644 index 0000000..cc54549 --- /dev/null +++ b/src/bem_grid_condat.cpp @@ -0,0 +1,104 @@ +/** + * + * @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 . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "bem_grid_condat.hh" +/* -------------------------------------------------------------------------- */ + +__BEGIN_TAMAAS__ + +BemGridCondat::BemGridCondat(Surface & surface): + BemGridKato(surface), + grad_step(1.8), + q_step(0.1) +{} + +/* -------------------------------------------------------------------------- */ + +BemGridCondat::~BemGridCondat() {} + +/* -------------------------------------------------------------------------- */ + +void BemGridCondat::updateGradient(Grid &q_vec) { + Real tmp[3] = {0}; + VectorProxy + q(q_vec.getInternalData(), 3), + g(nullptr, 3), + q_N(tmp, 3); + const UInt N = gaps.getNbPoints(); + + q_N = q; + q_N *= q_step; +#pragma omp parallel for firstprivate(g) + for (UInt i = 0 ; i < N ; i++) { + g.setPointer(gaps.getInternalData() + i * g.dataSize()); + g *= grad_step / lipschitz; + g += q; + } +} + +/* -------------------------------------------------------------------------- */ + +void BemGridCondat::computeEquilibrium(Real tol, const Grid & p_target) { + this->computeInfluence(); + Real f = 0; + UInt n_iter = 0; + + TAMAAS_ASSERT(p_target.dataSize() == 3, "Applied pressure must have only 3 components"); + TAMAAS_ASSERT(q_step > 0 && q_step < 1, "q step should be in ]0, 1["); + TAMAAS_ASSERT(grad_step < 2 * q_step, "grad step should be less than 2 q steps"); + this->tractions.uniformSetComponents(p_target); + + Real q_mem[3] = {0}; + Real e_mem[3] = {0}; + Real pn_mem[3] = {0}; + Real pn1_mem[3] = {0}; + Real zero_mem[3] = {0}; + VectorProxy q(q_mem, 3), e(e_mem, 3), pn(pn_mem, 3), pn1(pn1_mem, 3), zero(zero_mem, 3); + Real norm = + p_target(0)*p_target(0) + + p_target(1)*p_target(1) + + p_target(2)*p_target(2); + + do { + computeDisplacementsFromTractions(); + computeGaps(displacements); + updateGradient(q); + beta = computeLambda(); + BemGrid::computeMeanValueError(tractions, zero, pn); + tractions -= gaps; + BemGrid::projectOnFrictionCone(tractions, mu); + BemGrid::computeMeanValueError(tractions, zero, pn1); + pn1 *= 2; + q -= p_target; + q += pn1; + q -= pn; + lambda -= beta; + f = computeStoppingCriterion(); + printFriendlyMessage(n_iter, f); + } while ((f > tol || e.l2norm() / std::sqrt(norm) > tol) && n_iter++ < this->max_iter); +} + +__END_TAMAAS__ diff --git a/src/bem_grid_kato.hh b/src/bem_grid_condat.hh similarity index 60% copy from src/bem_grid_kato.hh copy to src/bem_grid_condat.hh index 806cf15..98a2b52 100644 --- a/src/bem_grid_kato.hh +++ b/src/bem_grid_condat.hh @@ -1,81 +1,60 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ -#ifndef __BEM_GRID_KATO_HH__ -#define __BEM_GRID_KATO_HH__ +#ifndef __BEM_GRID_CONDAT_HH__ +#define __BEM_GRID_CONDAT_HH__ /* -------------------------------------------------------------------------- */ -#include "surface.hh" -#include "grid.hh" -#include "types.hh" -#include "bem_grid.hh" +#include "bem_grid_kato.hh" /* -------------------------------------------------------------------------- */ - __BEGIN_TAMAAS__ /** - * @brief multi component implementation of Kato algorithm + * @brief multi component implementation of Condat algorithm * This method's main unknown is the pressure. + * Boundary condition is applied mean pressure. */ -class BemGridKato : public BemGrid { +class BemGridCondat : public BemGridKato { public: /// Constructor - BemGridKato(Surface & surface); + BemGridCondat(Surface & surface); /// Destructor - virtual ~BemGridKato(); + virtual ~BemGridCondat(); public: /* -------------------------------------------------------------------------- */ /* Actual computations */ /* -------------------------------------------------------------------------- */ - /// Compute true tractions - virtual void computeTrueTractions(); - /// Compute true displacements - virtual void computeTrueDisplacements(); - /// Compute stopping criterion - Real computeStoppingCriterion(); - /// Compute lambda (Lagrange mult. for p_N >= 0) - Real computeLambda(); - /// Compute potential energy - Real computeFunctional(); + /// Update gradient with a vector + virtual void updateGradient(Grid & q); /// Compute equilibrium virtual void computeEquilibrium(Real tol, const Grid & p_target); -/* -------------------------------------------------------------------------- */ -/* Others */ -/* -------------------------------------------------------------------------- */ - /// Print - void printFriendlyMessage(UInt niter, Real conv_crit); - /// Get Lambda - const Surface & getLambda() const {return lambda;} - protected: - Surface lambda; - Real beta; + Real grad_step, q_step; }; __END_TAMAAS__ - -#endif // __BEM_GRID_KATO_HH__ +#endif // __BEM_GRID_CONDAT_HH__ diff --git a/src/bem_grid_kato.cpp b/src/bem_grid_kato.cpp index 29ab3ab..505ca88 100644 --- a/src/bem_grid_kato.cpp +++ b/src/bem_grid_kato.cpp @@ -1,197 +1,211 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "bem_grid_kato.hh" #include "surface_statistics.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ BemGridKato::BemGridKato(Surface & surface): BemGrid(surface), lambda(surface.size(), surface.getL()), beta(0) { this->max_iter = 10000; } /* -------------------------------------------------------------------------- */ BemGridKato::~BemGridKato() {} /* -------------------------------------------------------------------------- */ void BemGridKato::computeTrueTractions() { // nothing to do here, tractions are our main unknown } void BemGridKato::computeTrueDisplacements() { // nothing to do here, tractions are our main unknown if (!this->true_displacements) this->true_displacements = new Grid(this->displacements.sizes(), this->displacements.lengths(), this->displacements.getNbComponents()); *this->true_displacements = this->displacements; const UInt N = this->true_displacements->getNbPoints(); Real *d_start = this->true_displacements->getInternalData(), *g_start = this->gaps.getInternalData(); VectorProxy d(nullptr, 3), g(nullptr, 3); #pragma omp parallel for firstprivate(d, g) for (UInt i = 0 ; i < N ; i++) { d.setPointer(d_start + i * d.dataSize()); g.setPointer(g_start + i * g.dataSize()); d(0) = g(0); d(1) = g(1); d(2) = g(2) + this->surface(i) + beta; } } /* -------------------------------------------------------------------------- */ Real BemGridKato::computeFunctional() { Real *p_start = this->tractions.getInternalData(), *u_start = this->displacements.getInternalData(); VectorProxy p(nullptr, 3), u(nullptr, 3); Real E = 0; const UInt N = this->tractions.getNbPoints(); #pragma omp parallel for reduction(+: E) firstprivate(p, u) for (UInt i = 0 ; i < N ; i++) { p.setPointer(p_start + i * p.dataSize()); u.setPointer(u_start + i * u.dataSize()); E += 0.5 * p.dot(u) - p(2) * this->surface(i); } return E; } /* -------------------------------------------------------------------------- */ Real BemGridKato::computeLambda() { Real *g_start = this->gaps.getInternalData(); VectorProxy g(nullptr, 3), g_T(nullptr, 2); const UInt N = this->gaps.getNbPoints(); const UInt nb_components = this->gaps.getNbComponents(); Real beta = std::numeric_limits::max(); // Computing lambda and beta at the same time #pragma omp parallel for reduction(min: beta) firstprivate(g, g_T) for (UInt i = 0 ; i < N ; i++) { g.setPointer(g_start + i * nb_components); g_T.setPointer(g_start + i * nb_components); lambda(i) = g(2) - mu*g_T.l2norm(); if (lambda(i) < beta) beta = lambda(i); } return beta; } /* -------------------------------------------------------------------------- */ +void BemGridKato::updateGradient(Grid &q_vec) { + VectorProxy + q(q_vec.getInternalData(), 3), + g(nullptr, 3); + const UInt N = gaps.getNbPoints(); + +#pragma omp parallel for firstprivate(g) + for (UInt i = 0 ; i < N ; i++) { + g.setPointer(gaps.getInternalData() + i * g.dataSize()); + g += q; + } +} + +/* -------------------------------------------------------------------------- */ + void BemGridKato::computeEquilibrium(Real tol, const Grid & g_target) { this->computeInfluence(); Real f = 0; UInt n_iter = 0; - TAMAAS_ASSERT(g_target.dataSize() == 3, "Applied pressure must have only 3 components"); - Real surface_mean = SurfaceStatistics::computeAverage(this->surface, 0); - this->surface -= (surface_mean + g_target(2)); + TAMAAS_ASSERT(g_target.dataSize() == 3, "Imposed gap must have only 3 components"); do { this->computeDisplacementsFromTractions(); this->computeGaps(this->displacements); + updateGradient(const_cast&>(g_target)); beta = computeLambda(); this->tractions -= this->gaps; BemGrid::projectOnFrictionCone(this->tractions, mu); lambda -= beta; f = computeStoppingCriterion(); printFriendlyMessage(n_iter, f); } while (f > tol && n_iter++ < this->max_iter); } /* -------------------------------------------------------------------------- */ /// Stopping criterion is othogonality between lambda and pressure Real BemGridKato::computeStoppingCriterion() { Real max = std::numeric_limits::min(); Real ortho = 0; Real *p_start = this->tractions.getInternalData(), *g_start = this->gaps.getInternalData(); VectorProxy p(nullptr, 3), p_T(nullptr, 2), g_T(nullptr, 2); const UInt N = this->tractions.dataSize() / this->tractions.getNbComponents(); const UInt nb_components = this->tractions.getNbComponents(); #pragma omp parallel firstprivate(p, g_T, p_T) { #pragma omp for reduction(max: max) for (UInt i = 0 ; i < this->tractions.dataSize() ; i++) { if (max < this->tractions(i)) max = this->tractions(i); } #pragma omp for reduction(+: ortho) for (UInt i = 0 ; i < N ; i++) { p.setPointer(p_start + i * nb_components); p_T.setPointer(p_start + i * nb_components); g_T.setPointer(g_start + i * nb_components); ortho += std::abs(p(2)*lambda(i)) + std::abs((mu*p(2) - p_T.l2norm())*g_T.l2norm()); } } // end #pragma omp parallel return ortho / max; } /* -------------------------------------------------------------------------- */ void BemGridKato::printFriendlyMessage(UInt niter, Real conv_crit) { if (niter % this->dump_freq == 0) { UInt prec = std::cout.precision(); std::cout << niter << " " << std::scientific << conv_crit << " " << computeFunctional() << std::endl; std::cout << std::fixed; std::cout.precision(prec); } if (std::isnan(conv_crit)) std::cerr << "Something went wrong in the solve phase !" << std::endl; } __END_TAMAAS__ diff --git a/src/bem_grid_kato.hh b/src/bem_grid_kato.hh index 806cf15..1ca47cc 100644 --- a/src/bem_grid_kato.hh +++ b/src/bem_grid_kato.hh @@ -1,81 +1,84 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __BEM_GRID_KATO_HH__ #define __BEM_GRID_KATO_HH__ /* -------------------------------------------------------------------------- */ #include "surface.hh" #include "grid.hh" #include "types.hh" #include "bem_grid.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief multi component implementation of Kato algorithm * This method's main unknown is the pressure. + * Boundary condition is applied mean gap. */ class BemGridKato : public BemGrid { public: /// Constructor BemGridKato(Surface & surface); /// Destructor virtual ~BemGridKato(); public: /* -------------------------------------------------------------------------- */ /* Actual computations */ /* -------------------------------------------------------------------------- */ /// Compute true tractions virtual void computeTrueTractions(); /// Compute true displacements virtual void computeTrueDisplacements(); /// Compute stopping criterion Real computeStoppingCriterion(); /// Compute lambda (Lagrange mult. for p_N >= 0) Real computeLambda(); /// Compute potential energy Real computeFunctional(); + /// Update gradient with a vector + virtual void updateGradient(Grid & q); /// Compute equilibrium - virtual void computeEquilibrium(Real tol, const Grid & p_target); + virtual void computeEquilibrium(Real tol, const Grid & g_target); /* -------------------------------------------------------------------------- */ /* Others */ /* -------------------------------------------------------------------------- */ /// Print void printFriendlyMessage(UInt niter, Real conv_crit); /// Get Lambda const Surface & getLambda() const {return lambda;} protected: Surface lambda; Real beta; }; __END_TAMAAS__ #endif // __BEM_GRID_KATO_HH__ diff --git a/src/bem_grid_teboulle.cpp b/src/bem_grid_teboulle.cpp index 136d08b..7b2adbf 100644 --- a/src/bem_grid_teboulle.cpp +++ b/src/bem_grid_teboulle.cpp @@ -1,83 +1,85 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "bem_grid_teboulle.hh" #include "surface_statistics.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ BemGridTeboulle::BemGridTeboulle(Surface & s): BemGridKato(s), z(s.sizes(), s.lengths(), 3), tractions_old(s.sizes(), s.lengths(), 3) {} /* -------------------------------------------------------------------------- */ BemGridTeboulle::~BemGridTeboulle() {} /* -------------------------------------------------------------------------- */ void BemGridTeboulle::computeEquilibrium(Real tol, const Grid & g_target) { this->computeInfluence(); Real f = 0, l = 0; Real t[2] = {1, 1}; UInt n_iter = 0; TAMAAS_ASSERT(g_target.dataSize() == 3, "Applied pressure must have only 3 components"); - Real surface_mean = SurfaceStatistics::computeAverage(this->surface, 0); - this->surface -= (surface_mean + g_target(2)); + // Real surface_mean = SurfaceStatistics::computeAverage(this->surface, 0); + // this->surface -= (surface_mean + g_target(2)); do { applyInfluenceFunctions(z, displacements); computeGaps(displacements); + updateGradient(const_cast&>(g_target)); beta = computeLambda(); gaps /= lipschitz; z -= gaps; BemGrid::projectOnFrictionCone(z, mu); tractions = z; t[1] = (1. + std::sqrt(4*t[0]*t[0] + 1)) / 2.; l = 1 + (t[0]-1)/t[1]; z -= tractions_old; z *= l; z += tractions_old; t[0] = t[1]; tractions_old = tractions; lambda -= beta; f = computeStoppingCriterion(); printFriendlyMessage(n_iter, f); } while (f > tol && n_iter++ < this->max_iter); - computeDisplacementsFromTractions(); - computeGaps(displacements); + gaps *= lipschitz; + //computeDisplacementsFromTractions(); + //computeGaps(displacements); } __END_TAMAAS__ diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp index 6eada49..7e85655 100644 --- a/src/bem_polonski.cpp +++ b/src/bem_polonski.cpp @@ -1,509 +1,507 @@ /** * * @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 "fft_plan_manager.hh" #include "fftransform.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()), pold(p.size(),p.getL()), p_sat(-1.){ e_star = 1.; max_iterations = 2000; + surface_rms_heights = SurfaceStatistics::computeStdev(this->surface); } /* -------------------------------------------------------------------------- */ BemPolonski::~BemPolonski() {} /* -------------------------------------------------------------------------- */ Real BemPolonski::computeEquilibrium(Real epsilon, Real pressure) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_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->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); 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); // std::cout << std::scientific << std::setprecision(10) // << SurfaceStatistics::computeAverage(surface_tractions) << std::fixed << std::endl; // 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::fixed << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); return f; } /* -------------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ - Real BemPolonski::computeEquilibriuminit(Real epsilon, Real pressure, Surface & init) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; this->surface_tractions=init; Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions); if (current_pressure <= 0.) surface_tractions = pressure; this->enforcePressureBalance(pressure); convergence_iterations.clear(); nb_iterations = 0; this->computeDisplacementsFromTractions(); f=this->computeF(); std::cout << std::scientific << std::setprecision(10) <<" " << f << std::endl; while(f > epsilon && nb_iterations < max_iterations) { fPrevious = f; this->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); 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); // std::cout << std::scientific << std::setprecision(10) // << SurfaceStatistics::computeAverage(surface_tractions) << std::fixed << std::endl; // 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::fixed << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); 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 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) 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 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) 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"); this->applyInfluenceFunctions(surface_t, surface_r); const UInt size = surface_t.size() * surface_t.size(); 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; } } //compute number of interpenetration without contact UInt nb_sl = 0; #pragma omp parallel for reduction(+: nb_sl) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) ==p_sat && 0.< this->gap(i) ) { this->surface_tractions(i) -= tau*this->gap(i); ++nb_sl; } } Real delta = 0; if (nb_iol > 0 || nb_sl> 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; // If we have no saturation, scale to match applied_pressure if (this->p_sat <= 0. || SurfaceStatistics::computeMaximum(this->surface_tractions) < this->p_sat) { Real pressure = 0.; #pragma omp parallel for reduction(+: pressure) for (UInt i = 0; i < size; ++i) { pressure += this->surface_tractions(i); } pressure *= 1./size; this->surface_tractions *= applied_pressure/pressure; } // If we have saturation, use secant method to find zero of saturationFunction else { // // Checking if we can find a zero // Real limit = // this->p_sat * SurfaceStatistics::computeContactArea(this->surface_tractions) // - applied_pressure; // if (limit < 0) SURFACE_FATAL("Saturation pressure is too small for applied pressure"); // // Initial points // Real x_n_2 = 0., x_n_1 = 1., x_n = 0.; // Real f_n_2 = 0., f_n_1 = 0., f_n = 0.; // // Secant loop // do { // f_n_2 = saturationFunction(x_n_2, applied_pressure); // f_n_1 = saturationFunction(x_n_1, applied_pressure); // x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2); // f_n = saturationFunction(x_n, applied_pressure); // x_n_2 = x_n_1; // x_n_1 = x_n; // } while (std::abs(f_n) > 1e-16); // this->surface_tractions *= x_n; // // Truncating pressures //#pragma omp parallel for // for (UInt i = 0 ; i < size ; i++) { // if (this->surface_tractions(i) > this->p_sat) // this->surface_tractions(i) = this->p_sat; // } Real mu_I_sat=0.; #pragma omp parallel for reduction(+: mu_I_sat) for (UInt i = 0; i < size; ++i) { if ( this->p_sat< this->surface_tractions(i) ){ mu_I_sat+=1.; } } Real pc=0.; #pragma omp parallel for reduction(+: pc) for (UInt i = 0; i < size; ++i) { if ( this->surface_tractions(i) < this->p_sat ){ pc+=this->surface_tractions(i);} } #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if ( this->surface_tractions(i) < this->p_sat ){ this->surface_tractions(i)= this->surface_tractions(i)*((size*applied_pressure-this->p_sat*mu_I_sat)/pc);} } #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if ( this->p_sat< this->surface_tractions(i) ){ this->surface_tractions(i)=this->p_sat; } } } STOPTIMER("enforcePressureBalance"); } /* -------------------------------------------------------------------------- */ Real BemPolonski::saturationFunction(Real alpha, Real applied_pressure) { const UInt n = surface.size(); const UInt size = n*n; Real sum = 0.; #pragma omp parallel for reduction(+: sum) for (UInt i = 0 ; i < size ; i++) { Real alpha_p = alpha * this->surface_tractions(i); Real alpha_p_I = (alpha_p > this->p_sat)? alpha_p - this->p_sat : 0.; sum += alpha_p - alpha_p_I; } sum /= size; return sum - applied_pressure; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeF() { UInt n = surface.size(); UInt size = n*n; Real res = 0; - Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions)); + Real t_sum = SurfaceStatistics::computeSum(surface_tractions); computeTrueDisplacements(); #pragma omp parallel for reduction(+:res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) == this->p_sat) continue; - res += std::abs(surface_tractions(i) * ( this->true_displacements(i)-surface(i))); + res += surface_tractions(i) * ( this->true_displacements(i)-surface(i)); // res += // this->surface_tractions[i].real() // *(surface_displacements[i].real() - surface[i].real()); } - return res/std::abs(t_sum); + return std::abs(res/t_sum/surface_rms_heights/size); } /* -------------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ void BemPolonski::computeTrueDisplacements(){ this->applyInfluenceFunctions(this->surface_tractions, this->surface_displacements); this->true_displacements = this->surface_displacements; UInt n = surface.size(); UInt size = n*n; Real shift = 1e300; #pragma omp parallel for reduction(min:shift) for (UInt i = 0; i < size; ++i) { if (surface_displacements(i) - shift - surface(i) < 0. &&( this->p_sat < 0 || this->surface_tractions(i) < this->p_sat )){ shift = surface_displacements(i) - surface(i); } } #pragma omp parallel for for (UInt i = 0; i < size; ++i) { true_displacements(i) = surface_displacements(i) - shift; } } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_polonski.hh b/src/bem_polonski.hh index 38b9e21..64aee47 100644 --- a/src/bem_polonski.hh +++ b/src/bem_polonski.hh @@ -1,106 +1,108 @@ /** * * @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_POLONSKI_H #define BEM_POLONSKI_H /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemPolonski : public BemFFTBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemPolonski(Surface & p); virtual ~BemPolonski(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: //! compute the equilibrium situation Real computeEquilibrium(Real epsilon,Real pressure); Real computeEquilibriuminit(Real epsilon,Real pressure, Surface & init); //! compute the mean of the gaps in the contact region Real computeMeanGapsInContact(); //! compute the non zero gap penalty functional Real computeG(); //! update the search direction void updateT(Real G, Real Gold, Real delta); //! compute the conjugate step size Real computeTau(); //! update the tractions Real updateTractions(Real tau); //! enforce the applied and contact pressures balance void enforcePressureBalance(Real applied_pressure); //! compute the gaps void computeGaps(); //! compute the gaps void computeTrueDisplacements(); //! compute F Real computeF(); //! saturation function Real saturationFunction(Real alpha, Real applied_pressure); //! set saturation pressure void setMaxPressure(Real p) {this->p_sat = p;} /* ------------------------------------------------------------------------ */ /* 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;}; Surface & getSurfaceT(){return surface_t;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! exploration direction for the polonski algo Surface surface_t; //! projected exploration direction for the polonski algo Surface surface_r; //! previous pressure Surface pold; //! saturation pressure Real p_sat; + //! the rms of heights of the surface + Real surface_rms_heights; }; __END_TAMAAS__ #endif /* BEM_POLONSKI_H */ diff --git a/src/grid.cpp b/src/grid.cpp index 4389a1e..d12b5e4 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -1,114 +1,127 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template Grid::Grid(): data() { memset(this->n, 0, dim * sizeof(UInt)); memset(this->L, 0, dim * sizeof(Real)); this->nb_components = 1; } /* -------------------------------------------------------------------------- */ template Grid::Grid(const UInt n[dim], const Real L[dim], UInt nb_components): data() { memcpy(this->n, n, dim * sizeof(UInt)); memcpy(this->L, L, dim * sizeof(Real)); this->nb_components = nb_components; this->resize(this->n); } /* -------------------------------------------------------------------------- */ template Grid::~Grid() {} /* -------------------------------------------------------------------------- */ template void Grid::resize(const UInt *n) { if (n != this->n) memcpy(this->n, n, dim * sizeof(UInt)); UInt size = this->computeSize(); data.resize(size); data.assign(size, T(0.)); } /* -------------------------------------------------------------------------- */ +template +void Grid::uniformSetComponents(const Grid &vec) { + TAMAAS_ASSERT(vec.dataSize() == this->nb_components, + "Cannot set grid field with values of vector"); + +#pragma omp parallel for + for (UInt i = 0 ; i < dataSize() ; i++) { + data[i] = vec(i % nb_components); + } +} + +/* -------------------------------------------------------------------------- */ + template void Grid::printself(std::ostream & str) const { str << "Grid(" << dim << ", " << nb_components << ") {"; for (UInt i = 0 ; i < data.size() - 1 ; i++) { str << data[i] << ", "; } str << data[data.size()-1] << "}"; } /* -------------------------------------------------------------------------- */ #define GRID_SCALAR_OPERATOR_IMPL(op) \ template \ inline void Grid::operator op(const T & e) { \ _Pragma("omp parallel for") \ for (UInt i = 0 ; i < this->data.size() ; i++) { \ data[i] op e; \ } \ } \ GRID_SCALAR_OPERATOR_IMPL(+=); GRID_SCALAR_OPERATOR_IMPL(*=); GRID_SCALAR_OPERATOR_IMPL(-=); GRID_SCALAR_OPERATOR_IMPL(/=); GRID_SCALAR_OPERATOR_IMPL(=); #undef GRID_SCALAR_OPERATOR_IMPL /* -------------------------------------------------------------------------- */ /// Class instanciation #define GRID_INSTANCIATE_TYPE(type) template class Grid; \ template class Grid; \ template class Grid GRID_INSTANCIATE_TYPE(Real); GRID_INSTANCIATE_TYPE(UInt); GRID_INSTANCIATE_TYPE(complex); GRID_INSTANCIATE_TYPE(int); GRID_INSTANCIATE_TYPE(bool); GRID_INSTANCIATE_TYPE(unsigned long); #undef GRID_INSTANCIATE_TYPE __END_TAMAAS__ diff --git a/src/grid.hh b/src/grid.hh index 815a02a..422403d 100644 --- a/src/grid.hh +++ b/src/grid.hh @@ -1,304 +1,298 @@ /** * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __GRID_HH__ #define __GRID_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "array.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Multi-dimensional & multi-component array class */ template class Grid { /* -------------------------------------------------------------------------- */ /* Constructors */ /* -------------------------------------------------------------------------- */ public: /// Constructor by default (empty array) Grid(); /// Constructor Grid(const UInt n[dim], const Real L[dim], UInt nb_components); /// Destructor virtual ~Grid(); /* -------------------------------------------------------------------------- */ /* Iterator class */ /* -------------------------------------------------------------------------- */ struct iterator { /// constructor iterator(T * start, ptrdiff_t offset):ptr(start), offset(offset) {} /// destructor ~iterator() {} /// pre-increment iterator & operator++() { ptr += offset; return *this; } /// comparison bool operator<(const iterator & a) { return ptr < a.ptr; } /// increment with given offset iterator & operator+=(ptrdiff_t a) { ptr += a; return *this;} /// dereference iterator T * operator*() { return ptr; } /// needed for OpenMP range calculations Int operator-(const iterator & a) const { return (ptr - a.ptr)/offset; } private: T * ptr; ptrdiff_t offset; }; public: iterator begin(UInt n = 1) { return iterator(data.getPtr(), n); } iterator end(UInt n = 1) { return iterator(data.getPtr()+computeSize(), n); } public: /* -------------------------------------------------------------------------- */ /* Common operations */ /* -------------------------------------------------------------------------- */ /// Resize array void resize(const UInt n[dim]); - /// Compute size inline UInt computeSize() const; - /// Print void printself(std::ostream & str) const; - /// Get sizes const UInt * sizes() const {return n;} - /// Get total size UInt dataSize() const {return data.size(); } - /// Get number of points UInt getNbPoints() const {return dataSize() / getNbComponents(); } - /// Get lengths const Real * lengths() const {return L;} - /// Get number of components UInt getNbComponents() const {return nb_components;} /// Set number of components void setNbComponents(UInt n) {nb_components = n;} - /// Get internal data pointer (const) const T * getInternalData() const {return data.getPtr();} /// Get internal data pointer (non-const) T * getInternalData() {return data.getPtr();} + /// Set components + void uniformSetComponents(const Grid & vec); /* -------------------------------------------------------------------------- */ /* Access operators */ /* -------------------------------------------------------------------------- */ template inline T & operator()(T1... args); template inline const T & operator()(T1... args) const; inline T & operator()(UInt tuple[dim+1]); inline const T & operator()(UInt tuple[dim+1]) const; private: template inline UInt unpackOffset(UInt offset, UInt index, T1... rest) const; template inline UInt unpackOffset(UInt offset, UInt index) const; inline UInt computeOffset(UInt tuple[dim+1]) const; /* -------------------------------------------------------------------------- */ /* Arithmetic operators */ /* -------------------------------------------------------------------------- */ public: #define GRID_VEC_OPERATOR(op) \ template \ void operator op(const Grid & other) \ GRID_VEC_OPERATOR(+=); GRID_VEC_OPERATOR(*=); GRID_VEC_OPERATOR(-=); GRID_VEC_OPERATOR(/=); #undef GRID_VEC_OPERATOR #define GRID_SCALAR_OPERATOR(op) \ void operator op(const T & e) \ GRID_SCALAR_OPERATOR(+=); GRID_SCALAR_OPERATOR(*=); GRID_SCALAR_OPERATOR(-=); GRID_SCALAR_OPERATOR(/=); GRID_SCALAR_OPERATOR(=); #undef GRID_SCALAR_OPERATOR // = operator void operator=(const Grid & other); // = operator (not const input, otherwise gcc is confused) void operator=(Grid & other); // Copy data from another grid template void copy(const Grid & other); /* -------------------------------------------------------------------------- */ /* Member variables */ /* -------------------------------------------------------------------------- */ protected: UInt n[dim]; Real L[dim]; Array data; UInt nb_components; }; /* -------------------------------------------------------------------------- */ /* Inline function definitions */ /* -------------------------------------------------------------------------- */ template inline UInt Grid::computeSize() const { UInt size = 1; for (UInt i = 0 ; i < dim ; i++) size *= n[i]; size *= nb_components; return size; } /* -------------------------------------------------------------------------- */ template template inline UInt Grid::unpackOffset(UInt offset, UInt index, T1... rest) const { UInt size = sizeof...(T1); offset += index; offset *= (size == 1) ? this->nb_components : this->n[dim-size+1]; return unpackOffset(offset, rest...); // tail-rec bb } template template inline UInt Grid::unpackOffset(UInt offset, UInt index) const { return offset + index; } template template inline T & Grid::operator()(T1... args) { TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1, "operator() does not match dimension"); UInt offset = unpackOffset(0, args...); return this->data[offset]; } template template inline const T & Grid::operator()(T1... args) const { TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1, "operator() does not match dimension"); UInt offset = unpackOffset(0, args...); return this->data[offset]; } template inline UInt Grid::computeOffset(UInt tuple[dim+1]) const { UInt offset = 0; for (UInt i = 0 ; i < dim-1 ; i++) { offset += tuple[i]; offset *= n[i+1]; } offset += tuple[dim-1]; offset *= nb_components; return offset + tuple[dim]; } template inline T & Grid::operator()(UInt tuple[dim+1]) { UInt offset = computeOffset(tuple); return this->data[offset]; } template inline const T & Grid::operator()(UInt tuple[dim+1]) const { UInt offset = computeOffset(tuple); return this->data[offset]; } /* -------------------------------------------------------------------------- */ #define GRID_VEC_OPERATOR_IMPL(op) \ template \ template \ inline void Grid::operator op(const Grid & other) { \ TAMAAS_ASSERT(other.computeSize() == data.size(), "surface size does not match");\ _Pragma("omp parallel for") \ for (UInt i = 0 ; i < this->data.size() ; i++) { \ data[i] op other(i); \ } \ } \ GRID_VEC_OPERATOR_IMPL(+=); GRID_VEC_OPERATOR_IMPL(-=); GRID_VEC_OPERATOR_IMPL(*=); GRID_VEC_OPERATOR_IMPL(/=); #undef GRID_VEC_OPERATOR /* -------------------------------------------------------------------------- */ template void Grid::operator=(const Grid & other) { this->copy(other); } template void Grid::operator=(Grid & other) { this->copy(other); } /* -------------------------------------------------------------------------- */ template template void Grid::copy(const Grid & other) { resize(other.sizes()); #pragma omp parallel for for (UInt i = 0 ; i < data.size() ; i++) data[i] = other(i); } /* -------------------------------------------------------------------------- */ /* Stream operator */ /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Grid & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #endif // __GRID_HH__ diff --git a/src/surface.cpp b/src/surface.cpp index 9e6cd22..99c3547 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1,454 +1,445 @@ /** * * @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 #include #include #include #include /* -------------------------------------------------------------------------- */ #include "surface.hh" #include "surface_statistics.hh" /* -------------------------------------------------------------------------- */ #ifdef USING_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ std::string fname = "surface_height.pdf"; /* -------------------------------------------------------------------------- */ template void Surface::expandByLinearInterpolation() { UInt n = this->size(); Surface * old_s = new Surface(n,1.); int new_n = n*2; *old_s = *this; this->setGridSize(new_n); for (UInt i=0; i void Surface::expandByShannonInterpolation(int factor) { /* /!\ this->FFTTransform(); /!\ */ TAMAAS_EXCEPTION("Shannon interpolation needs old FFTTransform interface"); UInt n = this->size(); UInt N = std::pow(2,factor)*this->size(); std::cout << N << " " << factor << std::endl; Surface s = Surface(N,1.); for (UInt i = 0; i < n/2; ++i) for (UInt j = 0; j < n/2; ++j) { s(i,j) = this->at(i,j); } for (UInt i = n/2; i < n; ++i) for (UInt j = 0; j < n/2; ++j) { s(N-n+i,j) = this->at(i,j); } for (UInt i = n/2; i < n; ++i) for (UInt j = n/2; j < n; ++j) { s(N-n+i,N-n+j) = this->at(i,j); } for (UInt i = 0; i < n/2; ++i) for (UInt j = n/2; j < n; ++j) { s(i,N-n+j) = this->at(i,j); } // s.dumpToTextFile("temp.csv"); /* /!\ s.FFTITransform(); /!\ */ this->setGridSize(N); Real f = std::pow(2,factor); f = f*f; for (UInt i = 0; i < N*N; ++i) this->at(i) = f*s(i); // int new_n = 2*n; // Surface s(*this); // this->setGridSize(new_n); // // put the original signal with zeros for to be interpolated points // for (int i=0; iat(2*i,j*2) = s(i,j); // // do shannon interpolation // int rcut = 64; // //interpolate lines // for (int i=0; i= new_n) i2 -= new_n; // if (i2%2) continue; // Real d = fabs(k); // s[i+j*new_n] += s(i2,j)*sinc(d/2); // } // } // } // //interpolate columns // for (int i=0; i= new_n) j2 -= new_n; // if (j2%2) continue; // Real d = fabs(k); // s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2), // s[i+j2*new_n].imag()*sinc(d/2)); // } // } // } // // for (int i=0; i void Surface::expandByMirroring() { UInt n = this->size(); UInt new_n = n*2; Surface s(*this); this->setGridSize(new_n); for (UInt i=0; i n) i1 = new_n - i1; UInt j1 = j; if (j1 > n) j1 = new_n - j1; if (i == n || j == n) { this->at(i,j) = 0.; continue; } this->at(i,j) = s(i1,j1); } } } /* -------------------------------------------------------------------------- */ template void Surface::expandByPuttingZeros() { UInt n = this->size(); UInt new_n = n*2; Surface s(*this); this->setGridSize(new_n); for (UInt i=0; i< n; i++) { for (UInt j=0; j< n; j++) { this->at(2*i,j*2) = s(i,j); } } } /* -------------------------------------------------------------------------- */ template void Surface::contractByTakingSubRegion() { UInt n = this->size(); Surface s(*this); this->setGridSize(n/2); for (UInt i=0; i < n/2; i++) { for (UInt j=0; j < n/2; j++) { this->at(i,j) = s(i,j); } } } /* -------------------------------------------------------------------------- */ template void Surface::contractByIgnoringMidPoints() { UInt n = this->size(); Surface s(*this); this->setGridSize(n/2); for (UInt i=0; iat(i/2,j/2) = s(i,j); } } } /* -------------------------------------------------------------------------- */ // void Surface::dumpToTextFile(string filename, bool pr) { // cout << "Writing surface file " << filename << endl; // ofstream file(filename.c_str()); // file.precision(15); // if (pr) // { // for (int i=0; i::loadFromTextFile(std::string filename) { // std::string line; // std::ifstream myfile (filename.c_str()); // if (!myfile.is_open()){ // SURFACE_FATAL("cannot open file " << filename); // } // Surface * my_surface = NULL; // UInt n = 0; // if (myfile.is_open()) // { // while (! myfile.eof() ) // { // getline (myfile,line); // if (line[0] == '#') // continue; // if (line.size() == 0) // continue; // std::stringstream s; // s << line; // int i,j; // s >> i; // s >> j; // if (n < i) n = i; // if (n < j) n = j; // } // ++n; // std::cout << "read " << n << " x " << n << " grid points" << std::endl; // myfile.clear(); // forget we hit the end of file // myfile.seekg(0, std::ios::beg); // move to the start of the file // my_surface = new Surface(n); // while (! myfile.eof() ) // { // getline (myfile,line); // if (line[0] == '#') // continue; // if (line.size() == 0) // continue; // std::stringstream s; // s << line; // int i,j; // Real h,hIm; // s >> i; // s >> j; // s >> h; // s >> hIm; // my_surface->heights[i+j*n] = complex(h,hIm); // } // myfile.close(); // } // return my_surface; // } /* -------------------------------------------------------------------------- */ template void Surface::dumpToParaview(std::string filename) const { #ifdef USING_IOHELPER UInt n = this->size(); const Surface & data = *this; /* build 3D coordinates array to be usable with iohelper */ // cerr << "allocate coords for paraview" << endl; Real * coords = new Real[n*n*3]; //Real * coords = new Real[(2*n-1)*(2*n-1)*3]; //cerr << "allocate zcoords for paraview" << endl; Real * zcoords = new Real[n*n*3]; //Real * zcoords = new Real[(2*n-1)*(2*n-1)*3]; //cerr << "allocate zcoords2 for paraview" << endl; Real * zcoords2 = new Real[n*n*3]; // Real * zcoords2 = new Real[(2*n-1)*(2*n-1)*3]; int index = 0; // for (int i=-n+1; i Surface::Surface(UInt a, Real L) : Map2dSquare(a,L){ this->setGridSize(a); } /* -------------------------------------------------------------------------- */ template Surface::~Surface(){ } /* -------------------------------------------------------------------------- */ // void Surface::copy(Surface & s){ // this->n = s.n; // this->L = s.L; // this->heights = s.heights; // } /* -------------------------------------------------------------------------- */ template void Surface::recenterHeights(){ UInt n = this->size(); Real zmean = SurfaceStatistics::computeAverage(*this); for (UInt i = 0 ; i < n*n ; ++i) this->at(i) = this->at(i) - zmean; } /* -------------------------------------------------------------------------- */ template void Surface::generateWhiteNoiseSurface(Surface & sources, long random_seed){ + //set the random gaussian generator + std::mt19937 gen(random_seed); + std::normal_distribution<> gaussian_generator; - //set the random gaussian generator - std::mt19937 gen(random_seed); - std::normal_distribution<> gaussian_generator; - Surface * origins; - origins = new Surface(1024,1.); - UInt n = sources.size(); + UInt n = sources.size(); - for (UInt i = 0 ; i < 1024 ; ++i) - for (UInt j = 0 ; j < 1024 ; ++j){ - origins->at(i,j) = gaussian_generator(gen); - } - - for (UInt i = 0 ; i < n ; ++i) - for (UInt j = 0 ; j < n ; ++j){ - sources(i,j) = origins->at(i*1024/n,j*1024/n); - } - - delete origins; + for (UInt i = 0 ; i < n ; ++i) + for (UInt j = 0 ; j < n ; ++j){ + sources(i,j) = gaussian_generator(gen); + } } /* -------------------------------------------------------------------------- */ template class Surface; template class Surface; template class Surface; template class Surface; template class Surface; __END_TAMAAS__ diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 0e00325..20a86f1 100755 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -1,63 +1,70 @@ #!/bin/bash ######################################################################### function test_output { name=$( basename $1 .verified ) - echo "Testing $name ouput" + echo "Testing $name output" echo -n "------------ " diff $name.out $name.verified >/dev/null 2>&1 return $? } ######################################################################### function test_return { name=$( basename $1 .py ) echo "Testing $name.py" echo -n "------------ " PYTHONPATH=../python python $name.py 2>/dev/null > $name.out return $? } ######################################################################### tamaas_path=$(PYTHONPATH=../python python -c 'import tamaas; print tamaas' | cut -d "'" -f 4 ) +tamaas_lib_path=$(ldd ../python/_tamaas.so | grep -e 'libTamaas.so' | cut -d " " -f 3) echo "Tamaas test suite" echo "=================" echo "" echo "importing tamaas from $tamaas_path" +echo "importing library from $tamaas_lib_path" echo "" n_tests=0 passed_tests=0 for test_file in test_*.py; do ((n_tests++)) test_return $test_file result=$? if [ $result -ne 0 ]; then echo "failure" else echo "success" ((passed_tests++)) fi done for test_file in test_*.verified; do + # In case no .verified files exist + if [ $test_file = "test_*.verified" ]; then + break + fi + ((n_tests++)) test_output $test_file result=$? if [ $result -ne 0 ]; then echo "failure" else echo "success" ((passed_tests++)) fi done 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_fractal.verified b/tests/test_fractal.verified index fc7ca3a..913a051 100644 --- a/tests/test_fractal.verified +++ b/tests/test_fractal.verified @@ -1,22 +1,22 @@ rms: 1 , hurst: 0.8 , q0: 4 , q1: 4 , q2: 32 -white noise average spectrum is 454.48 +white noise average spectrum is 453.34 [N] 512 [rms] 1 [rmsSpectral] 1.0 - [rmsSlopeSpectral] 55.8689894263 - [rmsSlopeGeometric] 55.7920762876 + [rmsSlopeSpectral] 53.6031186888 + [rmsSlopeGeometric] 53.5324509297 [Hurst] 0.8 [k1] 4 [k2] 32 [moment A] 45.9409927858 [m0] 1.0 [analytic m0] 1.0 - [m2] 1560.67198976 + [m2] 1436.64716658 [analytic m2] 1700.03942712 - [m4] 13639157.3671 + [m4] 12040443.7542 [analytic m4] 15108734.8282 - [alpha] 5.59969365371 + [alpha] 5.83367528856 [analytic_alpha] 5.22769343815 [seed] 156 diff --git a/version.py b/version.py new file mode 100644 index 0000000..2f5705f --- /dev/null +++ b/version.py @@ -0,0 +1,33 @@ +import subprocess +import base64 + +def write_info_file(file_name): + header = '''#include \nstd::string tamaas_release_info = "''' + file_content = '' + + branch = subprocess.check_output(["git", "rev-parse", "--abbrev-ref", "HEAD"])[:-1] + commit = subprocess.check_output(["git", "rev-parse", branch])[:-1] + + file_content += branch + ' ' + commit + '\n\n' + + diff = subprocess.check_output(["git", "diff"]) + diff_cached = subprocess.check_output(["git", "diff", "--cached"]) + + file_content += "diff\n" + base64.b64encode(diff) + '\n\n' + file_content += "diff --cached\n" + base64.b64encode(diff_cached) + '\n' + + remotes = subprocess.check_output(['git', 'remote', '-v']) + file_content += "remotes\n" + remotes + '\n' + + try: + origin_branch_commit = subprocess.check_output(["git", "rev-parse", "origin/" + branch]) + file_content += "origin/" + branch + " commit: " + origin_branch_commit + except: + file_content += "origin/" + branch + " does not exist" + + file_content += '"' + + file_content = file_content.replace('\n', '\\\n') + + with open(file_name, 'w') as file: + file.write(header + file_content + ';\n')