diff --git a/python/clusters.i b/python/clusters.i index eb4f209..f234044 100644 --- a/python/clusters.i +++ b/python/clusters.i @@ -1,17 +1,11 @@ /* -------------------------------------------------------------------------- */ /* Wrapping contact clusters classes */ /* -------------------------------------------------------------------------- */ %ignore operator <<; -%include "contact_area.hh" -%include "contact_cluster_collection.hh" %include "flood_fill.hh" -%include "std_vector.i" %include "std_list.i" namespace tamaas { -%template(VectorClusters) ::std::vector<::tamaas::ContactCluster>; -%template(ListClusters) ::std::list<::tamaas::Cluster>; + %template(ListClusters) ::std::list<::tamaas::Cluster>; } -%include "contact_cluster.hh" - diff --git a/python/surface.i b/python/surface.i index addcdef..0d6c495 100644 --- a/python/surface.i +++ b/python/surface.i @@ -1,388 +1,378 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ %{ #include <complex> -#include "map_2d.hh" #include "surface.hh" #include "types.hh" #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include <numpy/arrayobject.h> #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 "flood_fill.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<double>; //%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<double>, comlpexp); %define %my_pointer_class(NAME,TYPE) %template(SmartPointer##NAME) SmartPointer< TYPE >; %typemap(out) TYPE & { #define TOREPLACE SWIGTYPE_p_SmartPointerT_##NAME##_t SmartPointer<TYPE> * ptr = new SmartPointer<TYPE>($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(unsigned_int,tamaas::UInt) // dunno why this doesn't work %my_pointer_class(vectorReal,std::vector<tamaas::Real>) %apply (unsigned int IN_ARRAY1[ANY]) {(tamaas::UInt n[ANY])}; /* -------------------------------------------------------------------------- */ /* Wrapping Surface class */ /* -------------------------------------------------------------------------- */ %include "tamaas.hh" %ignore tamaas::Grid::iterator::operator++; %ignore tamaas::Grid::begin; %ignore tamaas::Grid::end; %ignore *::operator=; %include "grid.hh" %include "grid_hermitian.hh" namespace tamaas { %template (Grid1dReal) Grid<Real, 1>; %template (Grid2dReal) Grid<Real, 2>; %template (Grid2dBool) Grid<bool, 2>; %template (Grid2dComplex) Grid<Complex, 2>; %template (Grid2dInt) Grid<Int, 2>; %template (Grid2dUInt) Grid<UInt, 2>; %template (GridHermitian2dReal) GridHermitian<Real, 2>; } %include "types.hh" %include "grid_hermitian.hh" -%include "map_2d.hh" -namespace tamaas { -%template(Map2dComplex) Map2d<Complex>; -%template(Map2dInt) Map2d<int>; -%template(Map2dReal) Map2d<Real>; -} - %include "map_2d_square.hh" namespace tamaas { %template(Map2dSquareComplex) Map2dSquare<Complex>; %template(Map2dSquareInt) Map2dSquare<int>; %template(Map2dSquareReal) Map2dSquare<Real>; } %include "surface.hh" %include "surface_complex.hh" namespace tamaas { %template(SurfaceReal) Surface<Real>; %template(SurfaceInt) Surface<int>; %template(SurfaceRealComplex) SurfaceComplex<Real>; } %inline %{ __BEGIN_TAMAAS__ template <int required_type> 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 <typename T, UInt d> class GridForPython : public Grid<T, d> { public: GridForPython(T * wrapped_memory, UInt n[d], Real L[d], UInt nb_components): Grid<T, d>() { 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<T, d>() { 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 <typename T> class SurfaceForPython : public Surface<T>{ public: SurfaceForPython(T * wrapped_memory, UInt n, Real L) : Surface<T>(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<T>(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<DATA_TYPE *>($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<DATA_TYPE, DIM> & { 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 ("<class 'tamaas.GridForPython" + std::string("DATA_TYPE") + "'>" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } - else if ("<class 'tamaas.Grid" + std::string("DATA_TYPE") + "'>" == __type) { + else if ("<class 'tamaas.Grid" + std::string("DIM") + "d" + std::string("DATA_TYPE") + "'>" == __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<DATA_TYPE *>($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 ("<class 'tamaas.SurfaceForPython" + std::string("DATA_TYPE") + "'>" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else if ("<class 'tamaas.Surface" + std::string("DATA_TYPE") + "'>" == __type) { SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0); } else if ("<class 'tamaas.SurfaceComplex" + std::string("DATA_TYPE") + "'>" == __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>, Real, NPY_DOUBLE) %surface_typemaps(tamaas::SurfaceComplex<Real>, std::complex<Real>, NPY_COMPLEX128 ) %surface_typemaps(tamaas::Surface<unsigned long>, unsigned long, NPY_UINT ) %surface_typemaps(tamaas::Surface<UInt>, UInt, NPY_UINT ) %surface_typemaps(tamaas::Surface<int>, int, NPY_INT ) %grid_typemaps(tamaas::Grid1dReal, Real, NPY_DOUBLE, 1) %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2) %grid_typemaps(tamaas::Grid2dBool, bool, NPY_BOOL, 2) %grid_typemaps(tamaas::GridHermitian2dReal, tamaas::Complex, NPY_COMPLEX128, 2) %grid_typemaps(tamaas::Grid2dUInt, UInt, NPY_UINT, 2) %grid_typemaps(tamaas::SurfaceComplex<Real>, tamaas::Complex, NPY_COMPLEX128, 2) namespace tamaas { %template(SurfaceForPythonReal) SurfaceForPython< Real >; } %include "surface_timer.hh" diff --git a/python/surface_statistics.i b/python/surface_statistics.i index 651d118..89b3dfa 100644 --- a/python/surface_statistics.i +++ b/python/surface_statistics.i @@ -1,46 +1,38 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Wrapping surface statistics /* -------------------------------------------------------------------------- */ //%include "surface.i" %include "surface_statistics.hh" %template(computeDistribution) tamaas::SurfaceStatistics::computeDistribution<false>; %template(computeDistributionWithZeros) tamaas::SurfaceStatistics::computeDistribution<true>; %template(vecInt) std::vector<int>; //%include "contact_area.hh" //%include "contact_cluster_collection.hh" -namespace tamaas{ -%extend ContactArea -{ - Surface<int> & getSurface() { - return *$self; - } -} -} diff --git a/src/SConscript b/src/SConscript index f8078e6..8aa0cbe 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,78 +1,73 @@ Import('main_env') import os env = main_env # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_filter.cpp surface_generator_filter_fft.cpp isopowerlaw.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 flood_fill.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/array.hh b/src/array.hh index 06c748a..ea7d052 100644 --- a/src/array.hh +++ b/src/array.hh @@ -1,137 +1,134 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ - #ifndef __ARRAY__HH__ #define __ARRAY__HH__ /* -------------------------------------------------------------------------- */ -#include <cstring> #include "tamaas.hh" - /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <typename T> class Array { public: /// Copy constructor (deep) Array(const Array & v): _data(NULL), _size(0), wrapped(false) { this->resize(v._size); - memcpy(this->_data,v.getPtr(),this->_size*sizeof(T)); + std::copy(v.getPtr(), v.getPtr()+v._size, this->_data); } Array(): _data(NULL), _size(0), wrapped(false) {} Array(UInt size): _data(NULL), _size(0), wrapped(false) { this->resize(size); this->wrapped = false; } Array(T * data, UInt size): _data(NULL), _size(0), wrapped(false) { this->wrapMemory(data,size); } virtual ~Array(){ if (wrapped == false && _data != NULL){ delete [] _data; } } Array & operator=(const Array & v){ this->wrapped = false; if (v.size() != this->size()) this->resize(v.size()); memcpy(this->_data,v.getPtr(),this->_size*sizeof(T)); return *this; } void wrapMemory(T * data, UInt size) { this->_data = data; this->_size = size; this->wrapped = true; } void assign(UInt size, const T & value){ this->resize(size); for (UInt i = 0; i < size; i++) { _data[i] = value; } } const T * getPtr() const{ return _data; } T * getPtr() { return _data; } void setPtr(T * new_ptr) { _data = new_ptr;} void resize(UInt size){ if (size == 0 && wrapped == false) { if (_data != NULL) { delete [] this->_data; this->_data = NULL; } this->_size = 0; return; } if (size == this->_size) return; if (wrapped == true) TAMAAS_EXCEPTION("cannot resize wrapped surface"); if (_data != NULL) delete [] _data; this->_data = new T[size]; this->_size = size; } inline T & operator[](UInt i){ TAMAAS_ASSERT(i < _size, "memory overflow"); return _data[i]; } inline const T & operator[](UInt i)const{ TAMAAS_ASSERT(i < _size, "memory overflow"); return _data[i]; } inline UInt size() const{return _size;} private: T * _data; UInt _size; bool wrapped; }; __END_TAMAAS__ /* -------------------------------------------------------------------------- */ #endif /* __ARRAY__HH__ */ diff --git a/src/cluster_grow.cpp b/src/cluster_grow.cpp deleted file mode 100644 index 0be39b2..0000000 --- a/src/cluster_grow.cpp +++ /dev/null @@ -1,107 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#include "cluster_grow.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - -ClusterGrow::ClusterGrow(Map2d<int> & area, - int initial_state, - int to_check_state, - int contact_state, - int no_contact_state): - current_state (initial_state), - to_check_state (to_check_state), - contact_state (contact_state), - no_contact_state(no_contact_state), - area(area){ - -} - -/* -------------------------------------------------------------------------- */ - -ClusterGrow::~ClusterGrow(){ - -} - -/* -------------------------------------------------------------------------- */ - - -void ClusterGrow::addIndicesToExplore(std::vector<std::pair<UInt,UInt> > & - indices_to_add){ - - auto it = indices_to_add.begin(); - auto end = indices_to_add.end(); - - for( ; it != end ; ++it){ - int & val = area(*it); - - if (val == contact_state) { - - index_to_explore.push(*it); - - val = to_check_state; - } - } -} -/* -------------------------------------------------------------------------- */ - -template <bool periodic> -void ClusterGrow::grow(UInt i, UInt j){ - - if (area(i,j) != contact_state) return; - - ++current_state; - - area(i,j) = current_state; - - std::vector<std::pair<UInt,UInt> > v = area.getNeighborIndexes<periodic>(i,j); - this->addIndicesToExplore(v); - - while (!index_to_explore.empty()){ - - std::pair<UInt,UInt> ind = index_to_explore.front(); - - int & val = area(ind); - - if (val == contact_state || val == to_check_state){ - val = current_state; - std::vector<std::pair<UInt,UInt> > v = area.getNeighborIndexes<periodic>(ind); - this->addIndicesToExplore(v); - } - - index_to_explore.pop(); - - } -} - -/* -------------------------------------------------------------------------- */ - -template void ClusterGrow::grow<true>(UInt i, UInt j); -template void ClusterGrow::grow<false>(UInt i, UInt j); - -__END_TAMAAS__ diff --git a/src/cluster_grow.hh b/src/cluster_grow.hh deleted file mode 100644 index 31b6d7c..0000000 --- a/src/cluster_grow.hh +++ /dev/null @@ -1,105 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#ifndef __CLUSTER_GROW_HH__ -#define __CLUSTER_GROW_HH__ -/* -------------------------------------------------------------------------- */ -#include "contact_area.hh" -/* -------------------------------------------------------------------------- */ -#include <queue> -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - - -class ClusterGrow { - - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - -public: - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - - ClusterGrow(Map2d<int> & area, - int initial_state, - int to_check_state = -1, - int contact_state = 1, - int no_contact_state = 0); - virtual ~ClusterGrow(); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -public: - - //! start grow of a cluster from a given point - template <bool periodic> - void grow(UInt i, UInt j); - -private: - - void addIndicesToExplore(std::vector<std::pair<UInt,UInt> > & indices_to_add); - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ - -public: - - UInt getNbClusters(){ - UInt tmp = current_state - 2; - return tmp; - }; - - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ - -private: - - int current_state; - - const int to_check_state ; - const int contact_state ; - const int no_contact_state; - - Map2d<int> & area; - - std::queue<std::pair<UInt,UInt> > index_to_explore; -}; - -/* -------------------------------------------------------------------------- */ - -__END_TAMAAS__ - - -#endif /* __CLUSTER_GROW_HH__ */ diff --git a/src/contact_area.cpp b/src/contact_area.cpp deleted file mode 100644 index fbdd3bb..0000000 --- a/src/contact_area.cpp +++ /dev/null @@ -1,146 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#include"contact_area.hh" -#include <iostream> -#include "cluster_grow.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - -ContactArea::ContactArea(UInt n, Real L): Surface<int>(n,L) { - //num_cluster = 0; -} - -/* -------------------------------------------------------------------------- */ - -template <> -ContactArea::ContactArea(const Map2dSquare<Real> & m) : - Surface<int>(m.size(),m.getL()){ - - UInt n = this->size(); - - for (UInt i = 0 ; i < n*n ; ++i) - this->at(i) = (m(i) > 0.); - -} - -/* -------------------------------------------------------------------------- */ - - -ContactArea::~ContactArea() {} - -/* -------------------------------------------------------------------------- */ - -void ContactArea::detectContactClusters() { - - // starting state was push as 2 by vlad - ClusterGrow clg(*this,2); - UInt N = this->size(); - - for (UInt i = 0; i < N; ++i) { - for (UInt j = 0; j < N; ++j) { - clg.grow<true>(i,j); - } - } - - for (UInt i = 0; i < N; ++i) { - for (UInt j = 0; j < N; ++j) { - if (this->at(i,j) == -1) - SURFACE_FATAL("problem " << i << " " << j << " is -1"); - } - } - - for (UInt i = 0; i < clg.getNbClusters(); ++i) - colors.push_back(i+3); -} -/* -------------------------------------------------------------------------- */ - - -void ContactArea::replaceColors(int from, int to) { - // Repalce color "from" by the color "to" - UInt n = this->size(); - for (UInt i = 0; i < n; i++) - for (UInt j = 0; j < n; j++) - if (this->at(i,j) == from) - this->at(i,j) = to; -} -/* -------------------------------------------------------------------------- */ - -void ContactArea::randomizeClusterID(UInt seed) { - std::cout << "/ Randomize colors of contact clusters." << std::endl; - UInt shift=1000, tryid; - bool found; - std::vector<UInt> cid; - UInt num_cluster = colors.size(); - cid.resize(num_cluster); - while (shift < num_cluster) - shift*=10; - srand(seed); - for (int i = 0; i < int(colors.size()); i++) { - tryid=-1; - do { - tryid = shift+int((num_cluster+10)*random()/Real(RAND_MAX)); - found = true; - for (int j = 0; j < i; j++) - if (cid[j] == tryid) - found = false; - } - while (!found); - cid[i]=tryid; - replaceColors(colors[i],tryid); - colors[i]=tryid; - } - // Shift all colors back to 1 - for (int i = 0; i < int(colors.size()); i++) { - replaceColors(colors[i],colors[i]-shift+1); - colors[i]-=shift-1; - } -} - -/* -------------------------------------------------------------------------- */ - -std::vector<ContactCluster> ContactArea::getContactClusters() { - - std::vector<ContactCluster> clusters; - UInt nb_colors = this->colors.size(); - - for (UInt i = 0; i < nb_colors; i++) { - // std::cout << "Extracting contact cluster " << i << std::endl; - ContactCluster c(colors[i], *this); - //std::cout << "Compute statistics of contact cluster " << i << std::endl; - c.computeStatistics(); - clusters.push_back(c); - } - - return clusters; -} - -/* -------------------------------------------------------------------------- */ - - - -__END_TAMAAS__ diff --git a/src/contact_area.hh b/src/contact_area.hh deleted file mode 100644 index 6528c4d..0000000 --- a/src/contact_area.hh +++ /dev/null @@ -1,70 +0,0 @@ -//==================================================================// -// Percolation detection and characterization of contact clusters // -// To be integrated in FFT-BEM code // -// // -// V.A. Yastrebov, 2014 // -// Centre des Materiaux, CNRS, MINES ParisTech // -// // -// G. Anciaux, 2014 // -// EPFL, ENAC, IIC, LSMS // -//==================================================================// - -#ifndef __CONTACT_AREA__ -#define __CONTACT_AREA__ -/* -------------------------------------------------------------------------- */ -#include "map_2d.hh" -#include <vector> -#include "contact_cluster.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - - -class ContactArea : public Surface<int> { - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - -public: - - ContactArea(UInt n, Real L); - - template <typename T> - ContactArea(const Map2dSquare<T> & n); - ~ContactArea(); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - -public: - - //! TODO - void detectContactClusters(); - //! TODO - void replaceColors(int from, int to); - //! TODO - void randomizeClusterID(UInt seed); - //! generate contact clusters - std::vector<::tamaas::ContactCluster> getContactClusters(); - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ - - UInt getNumClusters(){return colors.size();}; - - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ - -private: - std::vector<int> colors; - -}; - -__END_TAMAAS__ - -#endif diff --git a/src/contact_cluster.cpp b/src/contact_cluster.cpp deleted file mode 100644 index 85bb6c2..0000000 --- a/src/contact_cluster.cpp +++ /dev/null @@ -1,405 +0,0 @@ -//==================================================================// -// Percolation detection and characterization of contact clusters // -// To be integrated in FFT-BEM code // -// // -// V.A. Yastrebov, 2014 // -// Centre des Materiaux, CNRS, MINES ParisTech // -// // -// G. Anciaux, 2014 // -// EPFL, ENAC, IIC, LSMS // -//==================================================================// - - -#include "contact_cluster.hh" -#include <iostream> -#include <cmath> -#include "contact_area.hh" -#include "cluster_grow.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - -void ContactCluster::computeBoundingBox(ContactArea & contact_area) { - - int N = contact_area.size(); - - int i0 = -1; - int j0 = -1; - int i1 = -1; - int j1 = -1; - - for (int i=0; i<N; i++) { - for (int j=0; j<N; j++) { - if (contact_area(i,j) == this->id && i0 < 0) i0 = i; - if (contact_area(j,i) == this->id && j0 < 0) j0 = i; - if (contact_area(N-i-1,N-j-1) == this->id && i1 < 0) i1 = N-i-1; - if (contact_area(N-j-1,N-i-1) == this->id && j1 < 0) j1 = N-i-1; - - // if all extremum are found then break - if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0) break; - } - // if all extremum are found then break - if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0) break; - } - - if (i0 < 0 || j0 < 0) SURFACE_FATAL("bounds not found " << id << " " << contact_area(0,0)); - if (i1 < 0 || j1 < 0) SURFACE_FATAL("bounds not found"); - - lb_corner[0] = i0; - lb_corner[1] = j0; - - rt_corner[0] = i1; - rt_corner[1] = j1; - - bbox_size[0] = i1-i0+1; - bbox_size[1] = j1-j0+1; - -} - - -/* -------------------------------------------------------------------------- */ - -void ContactCluster::checkLines(UInt dir,ContactArea & contact_area){ - - UInt N = contact_area.size(); - UInt i0 = lb_corner[dir]; - - start[dir] = -1; - stop[dir] = -1; - - for (UInt i = 0; i < N; i++) { - - bool fail = false; - for (UInt j = 0; j < N; j++){ - int val = -1; - if (dir == 0) val = contact_area(i,j); - if (dir == 1) val = contact_area(j,i); - - if (val == id) { - fail = true; - break; - } - } - - if (i0 == 0 && !fail && start[dir] < 0) { - start[dir] = i; - if (stop[dir] >= 0) SURFACE_FATAL("Error occured at "<<__LINE__); - } - if (i0 == 0 && fail && start[dir] >= 0 && stop[dir] < 0) - stop[dir] = i; - } - -} - -/* -------------------------------------------------------------------------- */ - - -void ContactCluster::minimizeBoundingBox(ContactArea & contact_area) { - - this->split_cluster = false; - UInt N = contact_area.size(); - if (bbox_size[0] != N && bbox_size[1] != N) return; - - // Minimize bounding box - this->split_cluster = true; - - // Check if there are vertical lines that do not contact "id" - this->checkLines(0,contact_area); - // Check if there are horizontal lines that do not contact "id" - this->checkLines(1,contact_area); - - if (start[0] >= 0) bbox_size[0] -= (stop[0]-start[0])-1; - if (start[1] >= 0) bbox_size[1] -= (stop[1]-start[1])-1; - -} - -/* -------------------------------------------------------------------------- */ - -void ContactCluster::extractCluster(ContactArea & contact_area) { - - UInt Nx = this->size(0); - UInt Ny = this->size(1); - UInt i0 = lb_corner[0]; - UInt j0 = lb_corner[1]; - UInt N = contact_area.size(); - - mass_center[0] = 0; - mass_center[1] = 0; - - for (UInt i = 0; i < Nx; i++) { - for (UInt j = 0; j < Ny; j++) { - if (i+i0 >= N || j+j0 >= N) - SURFACE_FATAL(" i0=" << i0 - << ",j0=" << j0 - << ", i=" << i - << ", j=" << j); - - int val = contact_area(i+i0,j+j0); - if (val == id) { - this->at(i,j) = 1; - mass_center[0] += i+i0; - mass_center[1] += j+j0; - } - } - } -} - -/* -------------------------------------------------------------------------- */ - -void ContactCluster::extractClusterPeriodic(ContactArea & contact_area) { - - // if the cluster is split by periodicity - mass_center[0] = 0; - mass_center[1] = 0; - UInt Nx = this->size(0); - UInt Ny = this->size(1); - UInt N = contact_area.size(); - UInt i0 = lb_corner[0]; - UInt j0 = lb_corner[1]; - - UInt Nxmax = start[0]<0?Nx:N; - UInt Nymax = start[1]<0?Ny:N; - - for (UInt i = 0; i < Nxmax; i++) { - for (UInt j = 0; j < Nymax; j++) { - if (i+i0>=0 && i+i0<N && - j+j0>=0 && j+j0<N) { - - int val = contact_area(i+i0,j+j0); - if (val == id) { - if (start[0] >= 0 && start[1] < 0) { - mass_center[1] += j+j0; - if ((int)i < start[0]) { - this->at(i,j) = 1; - mass_center[0]+=i; - } - else { - this->at(i-(stop[0]-start[0])+1,j) = 1; - mass_center[0] -= (N-i); - } - } - else if (start[1] >= 0 && start[0] < 0) { - mass_center[0]+=i+i0; - if ((int)j < start[1]) { - this->at(i,j) = 1; - mass_center[1]+=j; - } - else { - this->at(i,j-(stop[1]-start[1])+1) = 1; - mass_center[1] -= (N-j); - } - } - else if (start[1] >= 0 && start[0] >= 0) { - if ((int)j < start[1] && (int)i < start[0]) { - this->at(i,j) = 1; - mass_center[0] += i; - mass_center[1] += j; - } - else if ((int)j < start[1] && (int)i >= start[0]) { - this->at(i-(stop[0]-start[0])+1,j) = 1; - mass_center[1]+=j; - mass_center[0]-=(N-i); - } - else if ((int)j >= start[1] && (int)i >= start[0]) { - this->at(i-(stop[0]-start[0])+1,j-(stop[1]-start[1])+1) = 1; - mass_center[0] -= (N-i); - mass_center[1] -= (N-j); - } - else if ((int)j >= start[1] && (int)i < start[0]) { - this->at(i,j-(stop[1]-start[1])+1) = 1; - mass_center[0] += i; - mass_center[1] -= (N-j); - } - } - else if (start[0] < 0 && start[1] < 0) { - this->at(i,j) = 1; - mass_center[0] += i+i0; - mass_center[1]+=j+j0; - } - } - } - else { - std::cout<<"Error i0="<<i0<<",j0="<<j0 - <<",i="<<i<<",j="<<j - <<",Nxmax="<<Nxmax << ",Nymax" << Nymax - <<std::endl; - abort(); - } - } - } -} - -/* -------------------------------------------------------------------------- */ - - -ContactCluster::ContactCluster(int color, ContactArea & contact_area) - : Map2d(0,1.) { - - concave = false; - id = 0; - lb_corner[0] = 0.; - lb_corner[1] = 0.; - rt_corner[0] = 0.; - rt_corner[1] = 0.; - start[0] = 0.; - stop[1] = 0.; - bbox_size[0] = 0.; - bbox_size[1] = 0.; - A = 0.; - P = 0.; - IP = 0.; - mass_center[0] = 0; - mass_center[1] = 0; - Ps = 0.; - Cd = 0.; - //RI = 0.; - comp = 0.; - split_cluster = 0.; - - this->id = color; - - this->computeBoundingBox(contact_area); - - if (bbox_size[0] < 0 || bbox_size[1] < 0) { - SURFACE_FATAL("Something wrong with the detection of the bounding box of a contact cluster"); - } - - this->minimizeBoundingBox(contact_area); - this->setGridSize(bbox_size); - - if (!split_cluster) this->extractCluster(contact_area); - else this->extractClusterPeriodic(contact_area); -} - -/* -------------------------------------------------------------------------- */ - -void ContactCluster::computeStatistics() { - - bool start = false; - // UInt start_i, start_j; - A=0; - P=0; - IP=0; - - // UInt Nmax = bbox_size[0] > bbox_size[1] ? bbox_size[0]:bbox_size[1]; - UInt intersection[2] = {0,0}; - - for (UInt i = 0; i < bbox_size[0]; i++) { - intersection[0] = 0; - intersection[1] = 0; - for (UInt j = 0; j < bbox_size[1]; j++) { - // Area - A += this->at(i,j); - - // Perimeter at edges - if ( (i==0 || i== bbox_size[0] - 1) && this->at(i,j) == 1) { - if (!start) { - start=true; - // start_i=i; - //start_j=j; - } - P++; - } - if ( (j==0 || j== bbox_size[0] - 1) && - this->at(i,j) == 1) { - if (!start) { - start=true; - //start_i=i; - //start_j=j; - } - P++; - } - // Perimeter along lines - if (i>0 && this->at(i,j) != this->at(i-1,j)) { - if (!start) { - start=true; - //start_i=i; - //start_j=j; - } - P++; - } - // Perimeter along columns - if (j>0 && this->at(i,j) != this->at(i,j-1)) { - P++; - } - // Interpixel perimeter along lines - if (i>0 && this->at(i,j) == this->at(i-1,j) && this->at(i,j) > 0) - IP++; - // Interpixel perimeter along columns - if (j>0 && this->at(i,j) == this->at(i,j-1) && this->at(i,j) > 0) - IP++; - } - } - - for (UInt k = 0; k < 2; ++k) { - mass_center[k] /= Real(A); - //if (mass_center[k] < 0) { - // mass_center[k] += N; - //} - } - - this->comp = P/sqrt(Real(A)); - Cd = (A-P/4.)/(A-sqrt(1.*A)); - if (intersection[0] > 2 || intersection[1] > 2) - concave=true; - - this->findHoles(); -} - -/* -------------------------------------------------------------------------- */ - -void ContactCluster::findHoles() { - - UInt Nx = this->size(0); - UInt Ny = this->size(1); - - ClusterGrow clg(*this,id,-1,0,id); - - for (UInt i0 = 0; i0 < Nx; ++i0) { - for (UInt j0 = 0; j0 < Ny; ++j0) { - - // check if the point touch the border and is not in contact - if ( (i0 == 0 || j0 == 0 || - i0 == Nx-1 || j0 == Ny-1) && - this->at(i0,j0) == 0) { - - clg.grow<true>(i0,j0); - } - } - } - // find the holes - std::map<UInt,Real> holes_areas; - for (UInt i=0; i<Nx; i++) - for (UInt j=0; j<Ny; j++) - if (this->at(i,j) > id){ - if (holes_areas.count(this->at(i,j)) == 0) - holes_areas[this->at(i,j)] = 0; - holes_areas[this->at(i,j)] += 1; - } - - // this->RI = (Nx*Ny-Ai-A)/Real(A); -} - - -/* -------------------------------------------------------------------------- */ - - -void ContactCluster::printself(std::ostream & stream, int indent) const { - - stream << mass_center[0] << " " << mass_center[1] << " " - << bbox_size[0] << " " << bbox_size[1] << " " - << sqrt(A) << " " - << P << " " - << comp << " " - << IP << " " - << Cd << " " - << this->getNbHoles() << " " - // << RI << " " - << id; -} - -/* -------------------------------------------------------------------------- */ - - - -__END_TAMAAS__ diff --git a/src/contact_cluster.hh b/src/contact_cluster.hh deleted file mode 100644 index eaa1d79..0000000 --- a/src/contact_cluster.hh +++ /dev/null @@ -1,167 +0,0 @@ -//==================================================================// -// Percolation detection and characterization of contact clusters // -// To be integrated in FFT-BEM code // -// // -// V.A. Yastrebov, 2014 // -// Centre des Materiaux, CNRS, MINES ParisTech // -// // -// G. Anciaux, 2014 // -// EPFL, ENAC, IIC, LSMS // -//==================================================================// - -#ifndef __CONTACT_CLUSTER__ -#define __CONTACT_CLUSTER__ -/* -------------------------------------------------------------------------- */ -#include <map> -#include "surface.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - -class ContactArea; -class ContactClusterCollection; -/* -------------------------------------------------------------------------- */ - - -/* Class Contact Cluster class - - based on the publication: - - [1] Ernesto Bribiesca, - "An easy measure of compactness for 2D and 3D shapes", - Pattern Recognition 41 (2008) 543--554. - doi:10.1016/j.patcog.2007.06.029 - -*/ - -class ContactCluster : public Map2d<int> { - - - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - - friend class ContactClusterCollection; - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - -public: - - ContactCluster() : Map2d(0,1.) {} - ~ContactCluster() {} - ContactCluster(int color, ContactArea & contact_area); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - -public: - - void computeBoundingBox(ContactArea & contact_area); - void minimizeBoundingBox(ContactArea & contact_area); - void extractCluster(ContactArea & contact_area); - void extractClusterPeriodic(ContactArea & contact_area); - void computeStatistics(); - void findHoles(); - - void checkLines(UInt dir,ContactArea & contact_area); - - - - //! function to print the contain of the class - virtual void printself(std::ostream & stream, int indent = 0) const; - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ - -public: - - std::vector<Real> getHolesAreas() const{ - std::vector<Real> areas; - for (auto & x: holes_areas){ - areas.push_back(x.second); - } - return areas; - } - UInt getA(){return A;}; - UInt getP(){return P;}; - UInt getNbHoles() const {return holes_areas.size();}; - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ - -private: - - /* If the shape is for sure concave then concave=true, - however it may be concave even if concave=false. - */ - //! concavity flag - bool concave; - - //! color id of the cluster - int id; - - //! coordinates of the left bottom corner of the bounding box (i0,j0) - UInt lb_corner[2]; - //! coordinates of the right top corner of the bounding box (i1,j1) - UInt rt_corner[2]; - - //! from clusters connected through periodicity start stop is identifying inner part - int start[2]; - int stop[2]; - - - //! size of the bounding box (Nx,Ny) - UInt bbox_size[2]; - - //! contact area (A) - UInt A; - - //! contact perimeter computed discretely (P) - UInt P; - - //! interpixel perimeter (IP), the same as contact perimeter in notions of [1] - int IP; - - //! the holes areas - std::map<UInt,Real> holes_areas; - - //! center of mass (ic,jc) - Real mass_center[2]; - - //! corrected contact perimeter (Ps) - Real Ps; - - //! discrete compactness (Cd=(n-P/4)/(n-sqrt(n)) ) from [1] - Real Cd; - - /* Ratio of the internal (closed gap) area to the contact area - $$RI = (A0-Ai-A)/Real(A)$$, - where Ai is the area within the bounding box that surrounds contact cluster - */ - // Real RI; - - //! compactness (comp = Ps/sqrt(A)) - Real comp; - - - //! TODO - bool split_cluster; - -}; - -/* -------------------------------------------------------------------------- */ - -inline std::ostream & operator <<(std::ostream & stream, const ContactCluster & _this) -{ - _this.printself(stream); - return stream; -} - - -__END_TAMAAS__ - -#endif diff --git a/src/contact_cluster_collection.cpp b/src/contact_cluster_collection.cpp deleted file mode 100644 index b66aa92..0000000 --- a/src/contact_cluster_collection.cpp +++ /dev/null @@ -1,234 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -/* -------------------------------------------------------------------------- */ -#include "contact_cluster_collection.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - - -ContactClusterCollection::ContactClusterCollection(ContactArea & area): - area(area){ - - clusters = area.getContactClusters(); - - TotalArea = 0; - TotalPerimeter = 0; - NumberOfClustersWithHoles = 0; - - for (UInt i = 0; i < clusters.size(); ++i) { - TotalArea += clusters[i].getA(); - TotalPerimeter += clusters[i].getP(); - NumberOfClustersWithHoles += clusters[i].getNbHoles(); - } - -} - -/* -------------------------------------------------------------------------- */ - -ContactClusterCollection::~ContactClusterCollection(){ - -} - - -/* -------------------------------------------------------------------------- */ - -void ContactClusterCollection::dumpAllClusters(const std::string & prefix){ - - for (UInt i = 0; i < clusters.size(); i++) { - std::stringstream sstr; - sstr << prefix << "-" << i << ".csv"; - clusters[i].dumpToTextFile(sstr.str()); - } - -} - -/* -------------------------------------------------------------------------- */ - - - -void ContactClusterCollection::dumpClusterInfoToFile(const std::string & filename){ - - std::ofstream clusterData(filename.c_str()); - // Output info on every cluster separately - clusterData << "# N=" << area.size() << std::endl; - clusterData << "# Number of clusters=" << clusters.size() << std::endl; - clusterData << "# [1]Num [2]ic [3]jc [4]Nx [5]Ny [6]sqrt(A) [7]P [8]comp [9]IP [10]Cd [11]holes [12]RI [13]color [14]TotalNumberOfClusters [15]TotalArea [16]TotalPerimeter [17]NumberOfClustersWithHoles" << std::endl; - for (UInt i = 0; i < clusters.size(); i++) { - clusterData << i+1 << " " << clusters[i] << " " - << clusters.size() << " " - << TotalArea << " " - << TotalPerimeter << " " - << NumberOfClustersWithHoles << std::endl; - } - -} - -/* -------------------------------------------------------------------------- */ - -Real ContactClusterCollection::getTotalArea(){ - return TotalArea; -} - -/* -------------------------------------------------------------------------- */ - - -Real ContactClusterCollection::getTotalPerimeter(){ - return TotalPerimeter; -} - -/* -------------------------------------------------------------------------- */ - -Real ContactClusterCollection::getNbClustersWithHole(){ - return NumberOfClustersWithHoles; -} - -/* -------------------------------------------------------------------------- */ - -Real ContactClusterCollection::getNbClusters(){ - return clusters.size(); -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getAreas(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - // std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].A); - } - return res; -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getMassCenters(UInt i){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].mass_center[i]); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getPerimeters(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].P); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getCompactnesses(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].comp); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getInterPixelPerimeters(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].IP); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getDiscretePerimeters(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].Cd); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getNbHoles(){ - - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].getNbHoles()); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - -// std::vector<Real> ContactClusterCollection::getRIs(){ -// std::vector<Real> res; -// for (UInt i = 0; i < clusters.size(); i++) { -// std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; -// res.push_back(clusters[i].RI); -// } -// return res; - -// } - -/* -------------------------------------------------------------------------- */ - -std::vector<Real> ContactClusterCollection::getColors(){ - std::vector<Real> res; - for (UInt i = 0; i < clusters.size(); i++) { - std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl; - res.push_back(clusters[i].id); - } - return res; - -} - -/* -------------------------------------------------------------------------- */ - - - -__END_TAMAAS__ diff --git a/src/contact_cluster_collection.hh b/src/contact_cluster_collection.hh deleted file mode 100644 index 54d5592..0000000 --- a/src/contact_cluster_collection.hh +++ /dev/null @@ -1,97 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#ifndef __CONTACT_CLUSTER_COLLECTION_HH__ -#define __CONTACT_CLUSTER_COLLECTION_HH__ -/* -------------------------------------------------------------------------- */ -#include "contact_area.hh" -#include "contact_cluster.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - - -class ContactClusterCollection { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - - ContactClusterCollection(ContactArea & area); - virtual ~ContactClusterCollection(); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -public: - - /// dump information about clusters to file - void dumpClusterInfoToFile(const std::string & filename); - - Real getTotalArea(); - Real getTotalPerimeter(); - Real getNbClustersWithHole(); - Real getNbClusters(); - - - std::vector<Real> getAreas(); - std::vector<Real> getMassCenters(UInt i); - std::vector<Real> getPerimeters(); - std::vector<Real> getCompactnesses(); - std::vector<Real> getInterPixelPerimeters(); - std::vector<Real> getDiscretePerimeters(); - std::vector<Real> getNbHoles(); - //std::vector<Real> getRIs(); - std::vector<Real> getColors(); - - /// write one file per cluster - void dumpAllClusters(const std::string & prefix); - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -private: - - - //! the vector of clusters - std::vector<ContactCluster> clusters; - - //! the input area where to extract all the clusters - ContactArea & area; - - //! the total area - UInt TotalArea; - //! the total perimeter - UInt TotalPerimeter; - //! the number of clusters with holes - UInt NumberOfClustersWithHoles; - -}; - -/* -------------------------------------------------------------------------- */ - -__END_TAMAAS__ - -#endif /* __CONTACT_CLUSTER_COLLECTION_HH__ */ diff --git a/src/grid.cpp b/src/grid.cpp index ba1fc54..437f290 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -1,127 +1,128 @@ /** * * @author Lucas Frérot <lucas.frerot@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include <cstring> #include <complex> +#include <algorithm> /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template <typename T, UInt dim> Grid<T, dim>::Grid(): data() { - memset(this->n, 0, dim * sizeof(UInt)); - memset(this->L, 0, dim * sizeof(Real)); + std::fill(this->n, this->n+dim, 0); + std::fill(this->L, this->L+dim, 0); this->nb_components = 1; } /* -------------------------------------------------------------------------- */ template <typename T, UInt dim> Grid<T, dim>::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)); + std::copy(n, n+dim, this->n); + std::copy(L, L+dim, this->L); this->nb_components = nb_components; this->resize(this->n); } /* -------------------------------------------------------------------------- */ template <typename T, UInt dim> Grid<T, dim>::~Grid() {} /* -------------------------------------------------------------------------- */ template <typename T, UInt dim> void Grid<T, dim>::resize(const UInt *n) { if (n != this->n) - memcpy(this->n, n, dim * sizeof(UInt)); + std::copy(n, n+dim, this->n); UInt size = this->computeSize(); data.resize(size); data.assign(size, T(0.)); } /* -------------------------------------------------------------------------- */ template <typename T, UInt dim> void Grid<T, dim>::uniformSetComponents(const Grid<T, 1> &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 <typename T, UInt dim> void Grid<T, dim>::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 <typename T, UInt dim> \ inline void Grid<T, dim>::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<type, 1>; \ template class Grid<type, 2>; \ template class Grid<type, 3> 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/isopowerlaw.cpp b/src/isopowerlaw.cpp index 6ae78e3..1d9de73 100644 --- a/src/isopowerlaw.cpp +++ b/src/isopowerlaw.cpp @@ -1,57 +1,109 @@ /** * * @author Lucas Frérot <lucas.frerot@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "isopowerlaw.hh" #include "fftransform.hh" +#include <map> /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <UInt dim> void Isopowerlaw<dim>::computeFilter(GridHermitian<Real, dim> & filter_coefficients) const { const Real l[dim] = {1.}; Grid<Real, dim> wavevectors(filter_coefficients.sizes(), l, dim); FFTransform<Real, dim>::template computeFrequencies<true>(wavevectors); auto end = wavevectors.end(dim); #pragma omp parallel for for (auto it = wavevectors.begin(dim) ; it < end ; ++it) { VectorProxy<Real> q(*it, dim); UInt i = wavevectors.getNbPoints() - (end - it); filter_coefficients(i) = (*this)(q); } } +/* -------------------------------------------------------------------------- */ + template <UInt dim> -Real Isopowerlaw<dim>::analyticalRMS() const { +Real Isopowerlaw<dim>::rmsHeights() const { return std::sqrt(M_PI * ( (hurst + 1)/hurst * q1*q1 - 1./hurst * std::pow(q1, 2*(hurst+1)) * std::pow(q2, -2*hurst) - q0*q0 )); } +/* -------------------------------------------------------------------------- */ + +/** + * Analytical moments, cf. Yastrebov et al. (2015) + * "From infinitesimal to full contact between rough surfaces: Evolution + * of the contact area", appendix A + */ +template <> +std::vector<Real> Isopowerlaw<2>::moments() const { + std::map<UInt, Real> T; + T[0] = 2*M_PI; T[2] = M_PI; T[4] = 3 * M_PI / 4.; + Real xi = q0 / q1; + Real zeta = q2 / q1; + + std::vector<Real> moments; + + using namespace std; + for (int q : {0, 2, 4}) { + Real m = T[q] * pow(q1, q-2*hurst) * ((1-pow(xi, q+2))/(q+2) + + (pow(zeta, q-2*hurst)-1)/(q-2*hurst)); + moments.push_back(m); + } + + return moments; +} + +template <> +std::vector<Real> Isopowerlaw<1>::moments() const { + TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces"); + return std::vector<Real>(); +} + +/* -------------------------------------------------------------------------- */ + +template <UInt dim> +Real Isopowerlaw<dim>::alpha() const { + std::vector<Real> m = moments(); + return m[0] * m[2] / (m[1]*m[1]); +} + +/* -------------------------------------------------------------------------- */ + +template <UInt dim> +Real Isopowerlaw<dim>::rmsSlopes() const { + std::vector<Real> m = moments(); + return std::sqrt(M_PI * m[1] / 2); +} + + template class Isopowerlaw<1>; template class Isopowerlaw<2>; __END_TAMAAS__ diff --git a/src/isopowerlaw.hh b/src/isopowerlaw.hh index c5ad09e..412f1db 100644 --- a/src/isopowerlaw.hh +++ b/src/isopowerlaw.hh @@ -1,85 +1,98 @@ /** * * @author Lucas Frérot <lucas.frerot@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef ISOPOWERLAW_H #define ISOPOWERLAW_H /* -------------------------------------------------------------------------- */ #include "filter.hh" #include "grid_hermitian.hh" #include "types.hh" #include "tamaas.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /// Class representing an isotropic power law spectrum template <UInt dim> class Isopowerlaw : public Filter<dim> { public: /// Default constructor Isopowerlaw() {}; /// Destructor virtual ~Isopowerlaw() {}; public: /// Compute filter coefficients virtual void computeFilter(GridHermitian<Real, dim> & filter_coefficients) const; /// Compute a point of the PSD inline Real operator()(const VectorProxy<Real> & q_vec) const; /// Analytical rms of heights - Real analyticalRMS() const; + Real rmsHeights() const; + + /// Analytical moments + std::vector<Real> moments() const; + + /// Analytical Nayak's parameter + Real alpha() const; + + /// Analytical RMS slopes + Real rmsSlopes() const; public: TAMAAS_ACCESSOR(q0, UInt, Q0); TAMAAS_ACCESSOR(q1, UInt, Q1); TAMAAS_ACCESSOR(q2, UInt, Q2); TAMAAS_ACCESSOR(hurst, Real, Hurst); protected: UInt q0, q1, q2; Real hurst; }; +/* -------------------------------------------------------------------------- */ +/* Inline implementations */ +/* -------------------------------------------------------------------------- */ + template<UInt dim> inline Real Isopowerlaw<dim>::operator()(const VectorProxy<Real> & q_vec) const { const Real C = 1.; const Real q = q_vec.l2norm(); Real val; if (q < q0) val = 0; else if (q > q2) val = 0; else if (q < q1) val = C; else val = C * std::pow((q/q1), -2*(hurst+1)); return std::sqrt(val); } __END_TAMAAS__ #endif // ISOPOWERLAW_H diff --git a/src/map_2d.cpp b/src/map_2d.cpp deleted file mode 100644 index 1636479..0000000 --- a/src/map_2d.cpp +++ /dev/null @@ -1,141 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#include "map_2d.hh" -#include "fft_plan_manager.hh" -#include <fstream> -#include <sstream> -#include <fftw3.h> -#include <omp.h> -/* -------------------------------------------------------------------------- */ -#define TIMER -#include "surface_timer.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - - -template <typename T> -void Map2d<T>::rescaleHeights(Real C){ - /* compute actual RMS */ - for (UInt i = 0 ; i < this->n[0] ; ++i) - for (UInt j = 0 ; j < this->n[1] ; ++j){ - this->at(i,j) *= C; - } -} - - -/* -------------------------------------------------------------------------- */ - -template <typename T> -void Map2d<T>::setGridSize(UInt n0, UInt n1){ - const UInt n[] = {n0, n1}; - this->resize(n); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -void Map2d<T>::setGridSize(const UInt n[2]){ - this->resize(n); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -void Map2d<T>::operator=(const T & e){ - this->Grid<T, 2>::operator=(e); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -Map2d<T>::Map2d(UInt n0, UInt n1, Real L0, Real L1) { - setGridSize(n0,n1); - this->L[0] = L0; - this->L[1] = L1; -} -/* -------------------------------------------------------------------------- */ -template <typename T> -Map2d<T>::Map2d(const UInt * n, const Real *L) { - setGridSize(n); - this->L[0] = L[0]; - this->L[1] = L[1]; -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -Map2d<T>::~Map2d() -{} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -void Map2d<T>::dumpToTextFile(std::string filename) { - - UInt n1 = this->size(0); - UInt n2 = this->size(1); - std::cout << "Writing surface file " << filename << std::endl; - std::ofstream file(filename.c_str()); - if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename); - - file.precision(15); - for (UInt i = 0; i < n1; i++) { - for (UInt j=0; j< n2; j++) { - file << i << " " << j << " "; - this->writeValue(file,this->at(i,j)); - file << std::endl; - } - file << std::endl; - } - file.close(); -} - -/* -------------------------------------------------------------------------- */ -template <typename T> -void Map2d<T>::writeValue(std::ostream & os, const T & val) { - os << val; -} - -/* -------------------------------------------------------------------------- */ -template <> -void Map2d<Complex>::writeValue(std::ostream & os, const Complex & val) { - - os << std::scientific << val.real(); - os << " " << val.imag(); -} - -/* -------------------------------------------------------------------------- */ - -template class Map2d<Complex>; -template class Map2d<Real>; -template class Map2d<int>; -template class Map2d<unsigned int>; -template class Map2d<bool>; -template class Map2d<unsigned long>; - -__END_TAMAAS__ diff --git a/src/map_2d.hh b/src/map_2d.hh deleted file mode 100644 index 41643ef..0000000 --- a/src/map_2d.hh +++ /dev/null @@ -1,416 +0,0 @@ -/** - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * - * @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 <http://www.gnu.org/licenses/>. - * - */ -/* -------------------------------------------------------------------------- */ - -#ifndef __MAP_2D_HH__ -#define __MAP_2D_HH__ -/* -------------------------------------------------------------------------- */ -#include <vector> -#include <iostream> -#include <fftw3.h> -#include "tamaas.hh" -#include <fstream> -#include "grid.hh" -/* -------------------------------------------------------------------------- */ - -__BEGIN_TAMAAS__ - -template <typename T> -class Map2d : public Grid<T, 2> { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.); - Map2d(const UInt * n, const Real * L); - - virtual ~Map2d(); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -public: - - template <bool periodic> - std::vector<std::pair<UInt,UInt> > getNeighborIndexes(UInt i, UInt j); - template <bool periodic> - std::vector<std::pair<UInt,UInt> > getNeighborIndexes(std::pair<UInt,UInt> i); - - // //! return an iterator over the neighbors - // iterator beginNeighborhood(); - // //! return the end iterator over the neighbors - // iterator endNeighborhood(); - - //compute compute component by component product and put result in object - template <typename T1, typename T2> - void fullProduct(Map2d<T1> & s1,Map2d<T2> &s2); - //! set the grid size - void setGridSize(const UInt n[2]); - //! set the grid size - void setGridSize(UInt n0,UInt n1); - //! rescale height profile - void rescaleHeights(Real factor); - //! Write heights array to text file usable with plot programs - void dumpToTextFile(std::string filename); - //! set to scalar value - void operator=(const T & val); - //! copy method - template <template<typename> class surf, typename T2> - void copy(const surf<T2> & s); - // copy the map passed as parameter - template <template<typename> class surf, typename T2> - void operator=(const surf<T2> & s); - // temporary hack - void operator=(const Map2d<Complex> & s); - // copy the map passed as const parameter - template <template<typename> class surf, typename T2> - void operator=(surf<T2> & s); - -/* -------------------------------------------------------------------------- */ - -private: - - //! write an entry to a stringstream - void writeValue(std::ostream & os, const T & val); - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ - -public: - - // access to the i,j-th value in the surface - inline T & operator()(UInt i,UInt j); - // access to the i,j-th value in the surface - inline T & operator()(std::pair<UInt,UInt> i); - // access to the i,j-th value in the surface - inline const T & operator()(UInt i,UInt j)const; - // access to the i,j-th value in the surface - inline const T & operator()(std::pair<UInt,UInt> i)const; - // access to the i,j-th value in the surface - inline T & at(UInt i,UInt j); - // access to the i-th value in the array - inline T & at(UInt i); - // access to the i,j-th value in the surface - inline T & at(std::pair<UInt,UInt> i); - // access to the i,j-th value in the surface - inline const T & at(UInt i,UInt j) const; - // access to the i,j-th value in the surface - inline const T & at(std::pair<UInt,UInt> i) const; - // access to the i,j-th value in the surface - inline const T & at(UInt i) const; - // access to the i-th value in the array - inline T & operator()(int i); - // access to the i-th value in the array - // inline T & operator[](std::pair<UInt,UInt> i); - // access to the i-th value in the array - inline const T & operator()(int i) const; - //! get lateral number of discretization points - inline UInt size(UInt dir)const {return this->n[dir];}; - //! get both lateral number of discretization points - inline const UInt* sizes()const {return this->n;}; - //! get the lateral size of the sample - inline const Real * getLs() const {return this->L;}; - //! get the lateral size of the sample - inline Real getL(UInt dir) const {return this->L[dir];}; - //! get the lateral size of the sample - inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];}; - //! get the lateral size of the sample - inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;}; - - const T * getInternalData() const{return this->data.getPtr();}; - - //! Write heights array to text file usable with plot programs - template <bool consider_zeros=true> - void dumpToTextFile(std::string filename) const; -}; - -/* -------------------------------------------------------------------------- */ -//#define _at(I,J,N) I+J*N[1] -#define _at(I,J,N) (I)*(N[1])+(J) -/* -------------------------------------------------------------------------- */ -template <typename T> -inline T & Map2d<T>::at(UInt i,UInt j){ - return this->Grid<T, 2>::operator()(i, j, 0); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline const T & Map2d<T>::at(UInt i,UInt j)const { - return this->Grid<T, 2>::operator()(i, j, 0); -} -/* -------------------------------------------------------------------------- */ -#undef _at -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline T & Map2d<T>::at(UInt i){ - return this->Grid<T, 2>::operator()(i); -} - -/* -------------------------------------------------------------------------- */ -template <typename T> -inline const T & Map2d<T>::at(UInt i) const{ - return this->Grid<T, 2>::operator()(i); -} - -/* -------------------------------------------------------------------------- */ - - -template <typename T> -inline T & Map2d<T>::operator()(UInt i,UInt j){ - return this->at(i,j); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline T & Map2d<T>::at(std::pair<UInt,UInt> i){ - return this->at(i.first,i.second); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline const T & Map2d<T>::at(std::pair<UInt,UInt> i)const { - return this->at(i.first,i.second); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline T & Map2d<T>::operator()(std::pair<UInt,UInt> i){ - return this->at(i.first,i.second); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline const T & Map2d<T>::operator()(UInt i,UInt j)const{ - return this->at(i,j); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline const T & Map2d<T>::operator()(std::pair<UInt,UInt> i)const { - return this->at(i.first,i.second); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline T & Map2d<T>::operator()(int i){ - return this->Grid<T, 2>::operator()(i); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -inline const T & Map2d<T>::operator()(int i) const{ - return this->Grid<T, 2>::operator()(i); -} - -/* -------------------------------------------------------------------------- */ -// template <typename T> -// inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){ -// return this->at(c.first,c.second); -// } - -/* -------------------------------------------------------------------------- */ - -template <typename T> -template <bool periodic> -std::vector<std::pair<UInt,UInt> > Map2d<T>::getNeighborIndexes(UInt i, UInt j){ - - std::vector<std::pair<UInt,UInt> > vres; - - - std::pair<int,int> res; - - for (int ii = -1; ii < 2; ++ii) { - res.first = ii+static_cast<int>(i); - if (res.first < 0 && !periodic) continue; - if (res.first >= (int)this->n[0] && !periodic) continue; - if (res.first < 0 && periodic) res.first += this->n[0]; - if (res.first >= (int)this->n[0] && periodic) res.first -= this->n[0]; - - for (int jj = -1; jj < 2; ++jj) { - - res.second = jj+static_cast<int>(j); - if (res.second < 0 && !periodic) continue; - if (res.second >= (int)this->n[1] && !periodic) continue; - if (res.second < 0 && periodic) res.second += this->n[1]; - if (res.second >= (int)this->n[1] && periodic) res.second -= this->n[1]; - - if (ii == 0 && jj == 0) continue; - - { - UInt i = res.first; - UInt j = res.second; - - if (i+j*this->n[0] >= this->n[0]*this->n[1]) SURFACE_FATAL("out of bound access " - << i << " " << j << " " - << this->n[0] << " " << this->n[1] << " " - << i+j*this->n[0]); - } - vres.push_back(res); - } - } - return vres; -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -template <bool periodic> -std::vector<std::pair<UInt,UInt> > Map2d<T> -::getNeighborIndexes(std::pair<UInt,UInt> i){ - return this->getNeighborIndexes<periodic>(i.first,i.second); -} - -/* -------------------------------------------------------------------------- */ -template <typename T> -template <typename T1, typename T2> -void Map2d<T>::fullProduct(Map2d<T1> & s1,Map2d<T2> & s2){ - - auto n1 = s1.sizes(); - auto n2 = s2.sizes(); - - if (n1[0] != n2[0] || - n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices"); - - if (this->n[0] != n1[0] || - this->n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices") - - for (UInt i = 0 ; i < this->n[0] ; ++i) - for (UInt j = 0 ; j < this->n[1] ; ++j){ - this->at(i,j) = s1(i,j)* s2(i,j); - } -} -/* -------------------------------------------------------------------------- */ -template <typename T> -template <bool consider_zeros> -void Map2d<T>::dumpToTextFile(std::string filename) const { - - UInt m = this->n[0]; - UInt n = this->n[1]; - std::cout << "Writing surface file " << filename << std::endl; - std::ofstream file(filename.c_str()); - if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename); - - file.precision(15); - for (UInt i = 0; i < m; i++) { - for (UInt j=0; j< n; j++) { - if (!consider_zeros && - std::norm(this->at(i,j)) < zero_threshold*zero_threshold) - continue; - - file << i << " " << j << " "; - file << std::scientific << this->at(i,j); - file << std::endl; - } - file << std::endl; - } - file.close(); -} - - -// template <typename T> -// class Map2d<T>::iterator { - -// public: - -// iterator(std::vector<T> & data, -// std::queue<std::pair<UInt,UInt> > & indices): -// data(data),indices(indices){ -// } - -// inline void operator++(){ -// indices.pop(); -// } - -// inline T & operator *(){ -// return data[indices.front]; -// } - -// bool operator != (Map2d<T>::iterator & it){ -// return (indices.front == it.indices.front); -// } - -// private: - -// std::vector<T> & data; -// std::queue<std::pair<UInt,UInt> > indices; - -// bool finished; -// }; - - -/* -------------------------------------------------------------------------- */ -template <typename T> -template <template<typename> class surf, typename T2> -void Map2d<T>::copy(const surf<T2> & s){ - //this->setGridSize(s.sizes()); - //UInt m = this->sizes()[0]; - //UInt n = this->sizes()[1]; - //#pragma omp parallel for - //for (UInt i = 0; i < m*n; i++) { - // this->at(i) = s(i); - //} - this->Grid<T, 2>::copy(s); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -template <template<typename> class surf, typename T2> -void Map2d<T>::operator=(surf<T2> & s){ - this->copy(s); -} -/* -------------------------------------------------------------------------- */ -template <typename T> -template <template<typename> class surf, typename T2> -void Map2d<T>::operator=(const surf<T2> & s){ - this->copy(s); -} - -/* -------------------------------------------------------------------------- */ - -template <typename T> -void Map2d<T>::operator=(const Map2d<Complex> &){ - SURFACE_FATAL("Cannot copy a complex surface into incompatible type"); -} -template <> -inline void Map2d<Complex>::operator=(const Map2d<Complex> & s){ - this->copy(s); -} - -/* -------------------------------------------------------------------------- */ - -__END_TAMAAS__ - -#endif /* __MAP_2D_HH__ */ diff --git a/src/map_2d_square.cpp b/src/map_2d_square.cpp index 0c9dac7..e1eaece 100644 --- a/src/map_2d_square.cpp +++ b/src/map_2d_square.cpp @@ -1,106 +1,106 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "map_2d_square.hh" #include <fstream> /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <typename T> void Map2dSquare<T>::loadFromTextFile(std::string filename, int each) { std::string line; std::ifstream myfile (filename.c_str()); 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; UInt 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 n = int(n/Real(each)); setGridSize(n); while (! myfile.eof() ) { getline (myfile,line); if (line == "") continue; std::stringstream s; s << line; int i,j; s >> i; s >> j; T val = this->readValue(s); if (i % each == 0 && j % each == 0) { - this->at(i/each,j/each) = val; + (*this)(i/each,j/each) = val; } } myfile.close(); } } /* -------------------------------------------------------------------------- */ template <typename T> T Map2dSquare<T>::readValue(std::stringstream & sstr) { SURFACE_FATAL("To implement"); } /* -------------------------------------------------------------------------- */ template <> Complex Map2dSquare<Complex>::readValue(std::stringstream & sstr) { Real h = 0.,hIm = 0.; sstr >> h; sstr >> hIm; return Complex(h,hIm); } /* -------------------------------------------------------------------------- */ template class Map2dSquare<Complex>; template class Map2dSquare<Real>; template class Map2dSquare<int>; template class Map2dSquare<unsigned long>; __END_TAMAAS__ diff --git a/src/map_2d_square.hh b/src/map_2d_square.hh index 37282c1..3f6eb22 100644 --- a/src/map_2d_square.hh +++ b/src/map_2d_square.hh @@ -1,81 +1,93 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef __MAP_2D_SQUARE_HH__ #define __MAP_2D_SQUARE_HH__ /* -------------------------------------------------------------------------- */ -#include "map_2d.hh" +#include "grid.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <typename T> -class Map2dSquare : public Map2d<T> { +class Map2dSquare : public Grid<T, 2> { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - Map2dSquare(UInt n, Real L = 1.): Map2d<T>(n,n,L,L){}; - virtual ~Map2dSquare(){}; + Map2dSquare(UInt n, Real L = 1.): Grid<T, 2>() { + UInt dims[2] = {n, n}; + this->resize(dims); + this->L[0] = L; + this->L[1] = L; + } + + virtual ~Map2dSquare() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: //! get lateral number of discretization points - inline UInt size() const{return Map2d<T>::size(0);}; - //! load all n-th point of heights array from text file usable with plot programs - void loadFromTextFile(std::string filename, int each=1); + inline UInt size() const{return this->n[0];}; //! set the grid size - void setGridSize(UInt n){Map2d<T>::setGridSize(n,n);}; + void setGridSize(UInt n) { + UInt dims[2] = {n, n}; + this->resize(dims); + } + + void loadFromTextFile(std::string filename, int each); private: //! read an entry from a stringstream T readValue(std::stringstream & sstr); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! get the lateral size of the sample - inline Real getL()const {return Map2d<T>::getL(0);}; + inline Real getL() const {return this->L[0];}; //! get the lateral size of the sample - inline void setL(Real L){Map2d<T>::setL(L,L);}; + inline void setL(Real L) { + this->L[0] = L; + this->L[1] = L; + } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; __END_TAMAAS__ #endif /* __MAP_2D_SQUARE_HH__ */ diff --git a/src/surface.cpp b/src/surface.cpp index 99c3547..7aaca73 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1,445 +1,445 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include <iostream> #include <fstream> #include <sstream> #include <string> #include <cmath> #include <random> /* -------------------------------------------------------------------------- */ #include "surface.hh" #include "surface_statistics.hh" /* -------------------------------------------------------------------------- */ #ifdef USING_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ std::string fname = "surface_height.pdf"; /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::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<n; i++) { for (UInt j=0; j<n; j++) { (*this)(i,j) = (*old_s)(i,j); } } delete old_s; } Real sinc(Real x){ return sin(M_PI*x)/(M_PI*x); } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::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); + s(i,j) = (*this)(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); + s(N-n+i,j) = (*this)(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); + s(N-n+i,N-n+j) = (*this)(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(i,N-n+j) = (*this)(i,j); } // s.dumpToTextFile<false,true>("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); + (*this)(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; i<n; i++) // for (int j=0; j<n; j++) -// this->at(2*i,j*2) = s(i,j); +// (*this)(2*i,j*2) = s(i,j); // // do shannon interpolation // int rcut = 64; // //interpolate lines // for (int i=0; i<new_n; i++) { // for (int j=0; j<new_n; j++) { // if (i%2 == 0) continue; // // Real loop on local neighbors to interpolate // for (int k=-rcut; k<rcut+1; k++) { // int i2 = i + k; // if (i2 < 0) i2 += new_n; // if (i2 >= 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; i++) { // for (int j=0; j<new_n; j++) { // if (j%2 == 0) continue; // // Real loop on local neighbors to interpolate // for (int k=-rcut; k<rcut+1; k++) { // int j2 = j + k; // if (j2 < 0) j2 += new_n; // if (j2 >= 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<n; i++) { // // for (int j=0; j<n; j++) { // // s[2*i+j*2*new_n].real() = 0; // // s[2*i+j*2*new_n].im = 0; // // } // // } } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::expandByMirroring() { UInt n = this->size(); UInt new_n = n*2; Surface s(*this); this->setGridSize(new_n); for (UInt i=0; i<new_n; i++) { for (UInt j=0; j<new_n; j++) { UInt i1 = i; if (i1 > 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.; + (*this)(i,j) = 0.; continue; } - this->at(i,j) = s(i1,j1); + (*this)(i,j) = s(i1,j1); } } } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::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); + (*this)(2*i,j*2) = s(i,j); } } } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::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); + (*this)(i,j) = s(i,j); } } } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::contractByIgnoringMidPoints() { UInt n = this->size(); Surface s(*this); this->setGridSize(n/2); for (UInt i=0; i<n; i++) { for (UInt j=0; j<n; j++) { if (i%2) continue; if (j%2) continue; - this->at(i/2,j/2) = s(i,j); + (*this)(i/2,j/2) = s(i,j); } } } /* -------------------------------------------------------------------------- */ // void Surface<T>::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<n; i++) { // for (int j=0; j<n; j++) { // file << i << " " << j << " "; // if (heights[i+j*n].real() < 1e-10) // file << scientific << -1e10 << " " << 0; // else // file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im; // file << endl; // } // file << endl; // } // } // else // { // for (int i=0; i<n; i++) { // for (int j=0; j<n; j++) { // file << i << " " << j << " "; // file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im; // file << endl; // } // } // } // file.close(); // } /* -------------------------------------------------------------------------- */ // Surface * Surface<T>::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 <typename T> void Surface<T>::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<n; i++) { // for (int j=-n+1; j<n; j++) { for (int i=0; i<n; i++) { for (int j=0; j<n; j++) { coords[3*index] = 1.*i/n; coords[3*index+1] = 1.*j/n; coords[3*index+2] = 0; zcoords[3*index] = 0; zcoords[3*index+1] = 0; zcoords2[3*index] = 0; zcoords2[3*index+1] = 0; int i1 = (i%n); int j1 = (j%n); if (i1 < 0) i1 += n; if (j1 < 0) j1 += n; zcoords[3*index+2] = data(i1,j1).real(); zcoords2[3*index+2] = data(i1,j1).imag(); ++index; } } iohelper::DumperParaview dumper; dumper.setPoints(coords, 3, index, filename.c_str()); dumper.setConnectivity(NULL, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE); dumper.addNodeDataField("data.re",zcoords,3,index); dumper.addNodeDataField("data.im",zcoords2,3,index); dumper.setPrefix("./paraview"); dumper.init(); dumper.dump(); delete [] coords; delete [] zcoords; delete [] zcoords2; #endif //USING_IOHELPER } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template <typename T> Surface<T>::Surface(UInt a, Real L) : Map2dSquare<T>(a,L){ this->setGridSize(a); } /* -------------------------------------------------------------------------- */ template <typename T> Surface<T>::~Surface(){ } /* -------------------------------------------------------------------------- */ // void Surface<T>::copy(Surface & s){ // this->n = s.n; // this->L = s.L; // this->heights = s.heights; // } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::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; + (*this)(i) = (*this)(i) - zmean; } /* -------------------------------------------------------------------------- */ template <typename T> void Surface<T>::generateWhiteNoiseSurface(Surface<Real> & sources, long random_seed){ //set the random gaussian generator std::mt19937 gen(random_seed); std::normal_distribution<> gaussian_generator; UInt n = sources.size(); for (UInt i = 0 ; i < n ; ++i) for (UInt j = 0 ; j < n ; ++j){ sources(i,j) = gaussian_generator(gen); } } /* -------------------------------------------------------------------------- */ template class Surface<Real>; template class Surface<unsigned long>; template class Surface<UInt>; template class Surface<int>; template class Surface<bool>; __END_TAMAAS__ diff --git a/src/surface.hh b/src/surface.hh index 5ef47fa..d2f49b3 100644 --- a/src/surface.hh +++ b/src/surface.hh @@ -1,92 +1,92 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __SURFACE_HH__ #define __SURFACE_HH__ /* -------------------------------------------------------------------------- */ #include <string> #include <math.h> #include <iostream> #include <limits> #include <complex> /* -------------------------------------------------------------------------- */ #include "map_2d_square.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <typename T> class Surface : public Map2dSquare<T> { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Surface(UInt a, Real L); ~Surface(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ //! Write VTU file that is a set of points void dumpToParaview(std::string filename) const; //! create a new surface with linear interpolated points void expandByLinearInterpolation(); //! create a new surface with linear interpolated points void expandByShannonInterpolation(int factor); //! create a new surface with mirroring void expandByMirroring(); //! create a new surface with mirroring void expandByPuttingZeros(); //! contract the surface by ignoring extra data void contractByIgnoringMidPoints(); //! contract the surface by taking left bottom part of it void contractByTakingSubRegion(); //recenter the heights to zero value void recenterHeights(); // set all the members of the map to a given value - using Map2d<T>::operator=; + using Grid<T, 2>::operator=; static void generateWhiteNoiseSurface(Surface<Real> & surface, long random_seed); virtual const Surface<T> & getWrappedNumpy(){return *this;}; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #include "surface_complex.hh" #endif /* __SURFACE_HH__ */ diff --git a/src/surface_generator.hh b/src/surface_generator.hh index 09329ea..07332bd 100644 --- a/src/surface_generator.hh +++ b/src/surface_generator.hh @@ -1,64 +1,64 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef SURFACE_GENERATOR_H #define SURFACE_GENERATOR_H /* -------------------------------------------------------------------------- */ #include "grid.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /// Class generating random surfaces template<UInt dim> class SurfaceGenerator { public: /// Default constructor SurfaceGenerator(); /// Default destructor virtual ~SurfaceGenerator(); public: /// Build surface profile (array of heights) virtual Grid<Real, dim> & buildSurface() = 0; /// Set surface sizes void setSizes(const UInt n[dim]); - /// Get random seed - long & getRandomSeed() { return random_seed; } + + TAMAAS_ACCESSOR(random_seed, long, RandomSeed); /// Maintaining old interface friend class SurfaceGeneratorFilterFFT; protected: Grid<Real, dim> grid; long random_seed; }; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #endif diff --git a/src/surface_generator_filter_fft.hh b/src/surface_generator_filter_fft.hh index 099d565..3febbb8 100644 --- a/src/surface_generator_filter_fft.hh +++ b/src/surface_generator_filter_fft.hh @@ -1,80 +1,80 @@ /** * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef SURFACE_GENERATOR_FILTER_FFT_H #define SURFACE_GENERATOR_FILTER_FFT_H /* -------------------------------------------------------------------------- */ #include "surface_generator_filter.hh" #include "isopowerlaw.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ typedef Grid<Real, 2> Grid2dReal; /// Legacy class to maintain interface class SurfaceGeneratorFilterFFT { public: //! can construct a generator with default values SurfaceGeneratorFilterFFT(); //! destructor virtual ~SurfaceGeneratorFilterFFT(); public: //! get high frequency cut UInt & getQ2(){return pl_filter.getQ2();} //! get low frequency cut UInt & getQ0(){return pl_filter.getQ0();} //! get roll frequency cut UInt & getQ1(){return pl_filter.getQ1();} //! get hurst Real & getHurst(){return pl_filter.getHurst();} /// compute analytical RMS - Real analyticalRMS() const {return pl_filter.analyticalRMS(); } + Real analyticalRMS() const {return pl_filter.rmsHeights(); } /// Swig does not wrap the numpy otherwise Grid2dReal & buildSurface(); /// Get gridsize UInt & getGridSize() { return grid_size; } /// Get random seed long & getRandomSeed() { return generator.getRandomSeed(); } /// Get RMS Real & getRMS() { return rms; } /// Contrain RMS Real constrainRMS(); /// Init void Init() {} protected: SurfaceGeneratorFilter<2> generator; Isopowerlaw<2> pl_filter; UInt grid_size; Real rms; }; __END_TAMAAS__ #endif diff --git a/src/tamaas.hh b/src/tamaas.hh index 643a99c..b254164 100644 --- a/src/tamaas.hh +++ b/src/tamaas.hh @@ -1,105 +1,105 @@ /** * * @author Lucas Frérot <lucas.frerot@epfl.ch> * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef TAMAAS_HH #define TAMAAS_HH /* -------------------------------------------------------------------------- */ #include <exception> #include <string> #include <complex> #include <iostream> /* -------------------------------------------------------------------------- */ /// Namespace macros #define __BEGIN_TAMAAS__ namespace tamaas { #define __END_TAMAAS__ } /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ /// Common types definitions typedef double Real; typedef unsigned int UInt; typedef int Int; typedef std::complex<Real> Complex; static const Real zero_threshold = 1e-14; /* -------------------------------------------------------------------------- */ /// initialize tamaas (0 threads => let OMP_NUM_THREADS decide) void initialize(UInt num_threads = 0); /// cleanup tamaas void finalize(); /* -------------------------------------------------------------------------- */ class Exception : public std::exception { public: Exception(const std::string & mesg) throw() { msg = mesg; } virtual const char* what() const throw() { return msg.c_str(); } virtual ~Exception() throw() {} private: std::string msg; }; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ /* -------------------------------------------------------------------------- */ #define TAMAAS_EXCEPTION(mesg) { \ std::stringstream sstr; \ sstr \ << __FILE__ \ << ":" << __LINE__ << ":FATAL: " \ << mesg << std::endl; \ std::cerr.flush(); \ throw ::tamaas::Exception(sstr.str()); \ } #define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg) #if defined(TAMAAS_DEBUG) #define TAMAAS_ASSERT(cond, reason) do { \ if (!(cond)) { \ TAMAAS_EXCEPTION(#cond " assert failed: " << reason); \ } } while(0) #define TAMAAS_DEBUG_EXCEPTION(reason) TAMAAS_EXCEPTION(reason) #else -#define TAMAAS_ASSERT(mesg, reason) +#define TAMAAS_ASSERT(cond, reason) #define TAMAAS_DEBUG_EXCEPTION(reason) #endif #define TAMAAS_ACCESSOR(var, type, name) type & get##name() { return var; } \ void set##name(const type & new_var) { var = new_var; } /* -------------------------------------------------------------------------- */ #endif // TAMAAS_HH diff --git a/src/types.hh b/src/types.hh index bba9353..f8ac28a 100644 --- a/src/types.hh +++ b/src/types.hh @@ -1,106 +1,106 @@ /** * * @author Lucas Frérot <lucas.frerot@epfl.ch> * * @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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #ifndef __TYPES_HH__ #define __TYPES_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include <cstddef> /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template <typename T, UInt dim> class Proxy : public Grid<T, dim> { public: /// Default constructor Proxy() : Grid<T, dim>() {} /// Copy constructor Proxy(Proxy & p) : Grid<T, dim>() { if (p.n != this->n) - memcpy(this->n, p.n, dim * sizeof(UInt)); + std::copy(p.n, p.n+dim, this->n); this->nb_components = p.nb_components; this->data.wrapMemory(p.getInternalData(), p.dataSize()); } /// Destructor virtual ~Proxy() {} public: /// Set pointer void setPointer(T * ptr) {this->data.setPtr(ptr);} /// Get increment ptrdiff_t increment() const {return this->data.size();} }; // Forward declarations template <typename T> class MatrixProxy; template <typename T> class VectorProxy : public Proxy<T, 1> { public: /// Wrap on memory VectorProxy(T * mem, UInt n); /// Copy constructor VectorProxy(VectorProxy & v): Proxy<T, 1>(v) {} /// Destructor virtual ~VectorProxy(); public: /// Access const T & operator()(UInt i) const; /// Access (non-const) T & operator()(UInt i); /// Matrix-vector multiplication void mul(const MatrixProxy<T> & mat, const VectorProxy<T> & vec); /// Dot product T dot(const VectorProxy<T> & vec); /// L2 norm T l2norm() const; /// = operator using Grid<T, 1>::operator=; }; /* -------------------------------------------------------------------------- */ template <typename T> class MatrixProxy : public Proxy<T, 2> { public: /// Wrap on memory MatrixProxy(T * mem, UInt n, UInt m); /// Copy constructor MatrixProxy(MatrixProxy & v): Proxy<T, 2>(v) {} /// Destructor virtual ~MatrixProxy(); public: /// Access const T & operator()(UInt i, UInt j) const; /// Access (non-const) T & operator()(UInt i, UInt j); }; __END_TAMAAS__ #endif // __TYPES_HH__ diff --git a/tests/test_surface.py b/tests/test_surface.py index f7efd80..8fe12ca 100644 --- a/tests/test_surface.py +++ b/tests/test_surface.py @@ -1,147 +1,135 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot <lucas.frerot@epfl.ch> # # @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 <http://www.gnu.org/licenses/>. # ----------------------------------------------------------------------------- from __future__ import print_function, division import sys import numpy as np import tamaas as tm def compute_mse(values, mean): return np.mean((values - mean)**2) / values.size def main(): tm.initialize(1) # Spectrum parameters in wavelengths disctretization_size = 1. lambda_s = 8. lambda_r = 64. lambda_l = 64. box_size = 512 # Independent parameters Hurst = 0.8 hrms = 1. # Spectrum parameters in wavenumbers k_l = int(box_size // lambda_l) k_r = int(box_size // lambda_r) k_s = int(box_size // lambda_s) N = int(box_size // disctretization_size) - # Analytical moments, cf. Yastrebov et al. (2015) - # "From infinitesimal to full contact between rough surfaces: Evolution - # of the contact area", appendix A - xi = k_l / k_r - zeta = k_s / k_r - H = Hurst - T = {0: 2*np.pi, - 2: np.pi, - 4: 3*np.pi/4} - - m = {} - for q in [0, 2, 4]: - m[q] = T[q] * k_r**(q-2*H)*((1-xi**(q+2))/(q+2) + (zeta**(q-2*H)-1)/(q-2*H)) - - # Some spectral parameters - alpha = m[0]*m[4]/m[2]**2 - rms_slopes = np.sqrt(np.pi * m[2] / 2) - # Output setup realizations = 100 hrms_out = np.zeros(realizations) alpha_out = np.zeros(realizations) rms_slopes_out = {"spectral": np.zeros(realizations), "geometric": np.zeros(realizations)} + iso = tm.Isopowerlaw2D() + iso.getQ0().assign(k_l) + iso.getQ1().assign(k_r) + iso.getQ2().assign(k_s) + iso.getHurst().assign(Hurst) + sg = tm.SurfaceGeneratorFilter2D() + sg.setFilter(iso) + sg.setSizes([N, N]) + + # Analytical moments + m = iso.moments() + alpha = iso.alpha() + rms_slopes = iso.rmsSlopes() + for i in range(realizations): # Surface generator setup - sg = tm.SurfaceGeneratorFilterFFT() - sg.getGridSize().assign(N) - sg.getHurst().assign(Hurst) - sg.getRMS().assign(hrms) - sg.getQ0().assign(k_l) - sg.getQ1().assign(k_r) - sg.getQ2().assign(k_s) sg.getRandomSeed().assign(i) - sg.Init() # Generating surface - surface = sg.buildSurface() / sg.analyticalRMS() + surface = tm.convertGrid(sg.buildSurface()) / iso.rmsHeights() # Computing RMS heights # We know surface.mean() = 0, so we get an unbiased estimator here hrms_out[i] = np.sqrt(1 / (N**2-1) * np.sum(surface**2)) # Computing spectrum bandwidth surface /= box_size**2 # Normalize so that the computed PSD does not contain the box size moments = tm.SurfaceStatistics.computeMoments(surface) m0 = hrms_out[i]**2 / box_size**4 m2 = moments[0] m4 = moments[1] alpha_out[i] = m0*m4/m2**2 # Computing RMS slopes rms_slopes_out["spectral"][i] = tm.SurfaceStatistics.computeSpectralRMSSlope(surface) rms_slopes_out["geometric"][i] = tm.SurfaceStatistics.computeRMSSlope(surface) # print(rms_slopes_out) # print(rms_slopes) # Computing Mean Square Error print("""========= Analytical values ========= hrms = {} alpha = {} m0, m2, m4 = {} rms_slopes = {}""".format(hrms, alpha, m, rms_slopes)) print("========= MSE Calculations =========") data_association = [ (hrms_out, hrms, "hrms"), (alpha_out, alpha, "alpha"), (rms_slopes_out["spectral"], rms_slopes, "RMS Slopes Spectral"), (rms_slopes_out["geometric"], rms_slopes, "RMS Slopes Geometric") ] # Computing and printing MSE values mse_vals = [] for out, analytical, string in data_association: mse = np.sqrt(compute_mse(out, analytical)) / analytical print("Normalized error (MSE) for {}: {}".format(string, mse)) mse_vals.append(mse) # Checking tolerances tols = [5e-3, 5e-3, 1e-1, 1e-1] for mse, tol in zip(mse_vals, tols): if mse >= tol: print("{} >= {}".format(mse, tol)) return 1 tm.finalize() return 0 if __name__ == "__main__": sys.exit(main())