diff --git a/src/bem_fft_base.cpp b/src/bem_fft_base.cpp index f92a785..f659097 100644 --- a/src/bem_fft_base.cpp +++ b/src/bem_fft_base.cpp @@ -1,270 +1,270 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" #include "surface_statistics.hh" #include "meta_functional.hh" #include "elastic_energy_functional.hh" #include "fft_plan_manager.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ BemFFTBase::BemFFTBase(Surface & p): BemInterface(p), functional(new MetaFunctional(*this)), surface_tractions(p.size(),p.getL()), surface_displacements(p.size(),p.getL()), surface_spectral_input(p.size(),p.getL()), surface_spectral_output(p.size(),p.getL()), surface_spectral_influence_disp(p.size(),p.getL()), true_displacements(p.size(),p.getL()), old_displacements(p.size(),p.getL()), gap(p.size(),p.getL()), e_star(std::nan("")), tmp_coeff(1./M_PI){ // this->setNumberOfThreads(this->nthreads); max_iterations = 5000; dump_freq = 100; this->functional->addFunctionalTerm(new ElasticEnergyFunctional(*this)); std::cerr << this << " is setting up the surfaces"; std::cerr << &p << " has size " << p.size() << std::endl; } /* -------------------------------------------------------------------------- */ BemFFTBase::~BemFFTBase() { - // delete functional; + FFTPlanManager::get().clean(); } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getTractions() const{ return surface_tractions; } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getDisplacements()const{ return surface_displacements; } /* -------------------------------------------------------------------------- */ const Surface & BemFFTBase::getSpectralInfluenceOverDisplacement(){ return surface_spectral_influence_disp; } /* -------------------------------------------------------------------------- */ void BemFFTBase::setEffectiveModulus(Real e_star){ this->e_star = e_star; this->computeSpectralInfluenceOverDisplacement(); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeSpectralInfluenceOverDisplacement(){ STARTTIMER("computeSpectralInfluenceOverDisplacement"); if (std::isnan(tmp_coeff)) SURFACE_FATAL("constant tmp_coeff is nan"); if (std::isnan(e_star)) SURFACE_FATAL("constant e_star is nan"); UInt n = surface.size(); surface_spectral_influence_disp(0) = 0.; // everything but the external lines and central lines #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { for (UInt j = 1; j < n/2; ++j) { Real d = sqrt(i*i+j*j); surface_spectral_influence_disp(i,j) = 1./d; surface_spectral_influence_disp(n-i,j) = 1./d; surface_spectral_influence_disp(i,n-j) = 1./d; surface_spectral_influence_disp(n-i,n-j) = 1./d; } } // external line (i or j == 0) #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { Real d = i; surface_spectral_influence_disp(i,0) = 1./d; surface_spectral_influence_disp(n-i,0) = 1./d; surface_spectral_influence_disp(0,i) = 1./d; surface_spectral_influence_disp(0,n-i) = 1./d; } //internal line (i or j == n/2) #pragma omp parallel for for (UInt i = 1; i < n/2; ++i) { Real d = sqrt(i*i+n/2*n/2); surface_spectral_influence_disp(i,n/2) = 1./d; surface_spectral_influence_disp(n-i,n/2) = 1./d; surface_spectral_influence_disp(n/2,i) = 1./d; surface_spectral_influence_disp(n/2,n-i) = 1./d; } { //i == n/2 j == n/2 Real d = sqrt(2*n/2*n/2); surface_spectral_influence_disp(n/2,n/2) = 1./d; //i == 0 j == n/2 d = n/2; surface_spectral_influence_disp(0,n/2) = 1./d; //i == n/2 j == 0 surface_spectral_influence_disp(n/2,0) = 1./d; } UInt size = n*n; Real L = surface.getL(); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { surface_spectral_influence_disp(i) *= tmp_coeff/this->e_star*L; } STOPTIMER("computeSpectralInfluenceOverDisplacement"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeDisplacementsFromTractions(){ STARTTIMER("computeDisplacementsFromTractions"); this->applyInfluenceFunctions(surface_tractions,surface_displacements); STOPTIMER("computeDisplacementFromTractions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeTractionsFromDisplacements(){ STARTTIMER("computeDisplacementsFromTractions"); this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions); STOPTIMER("computeDisplacementFromTractions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::applyInfluenceFunctions(Surface & input, Surface & output){ STARTTIMER("applyInfluenceFunctions"); FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); plan1.forwardTransform(); surface_spectral_output *= surface_spectral_influence_disp; FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); plan2.backwardTransform(); STOPTIMER("applyInfluenceFunctions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::applyInverseInfluenceFunctions(Surface & input, Surface & output){ STARTTIMER("applyInfluenceFunctions"); FFTransform & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output); plan1.forwardTransform(); surface_spectral_output /= surface_spectral_influence_disp; surface_spectral_output(0) = 0; FFTransform & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output); plan2.backwardTransform(); STOPTIMER("applyInfluenceFunctions"); } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeGaps() { UInt size = surface.size(); #pragma omp parallel for for (UInt i = 0 ; i < size*size ; i++) { gap(i) = true_displacements(i) - surface(i); } } /* -------------------------------------------------------------------------- */ void BemFFTBase::computeTrueDisplacements(){ this->true_displacements = this->surface_displacements; UInt n = surface.size(); UInt size = n*n; Real shift = 1e300; for (UInt i = 0; i < size; ++i) { if (surface_displacements(i) - shift - surface(i) < 0.){ shift = surface_displacements(i) - surface(i); } } for (UInt i = 0; i < size; ++i) { true_displacements(i) = surface_displacements(i) - shift; } } /* -------------------------------------------------------------------------- */ void BemFFTBase::addFunctional(Functional *new_funct) { this->functional->addFunctionalTerm(new_funct); } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/fft_plan_manager.cpp b/src/fft_plan_manager.cpp index 1a7c39d..b825fe5 100644 --- a/src/fft_plan_manager.cpp +++ b/src/fft_plan_manager.cpp @@ -1,140 +1,152 @@ /** * * @author Guillaume Anciaux * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fft_plan_manager.hh" #include "fftransform_fftw.hh" #if defined(USE_CUFFT) #include "fftransform_cufft.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ FFTPlanManager::FFTPlanManager(){ } /* -------------------------------------------------------------------------- */ FFTPlanManager::~FFTPlanManager(){ + clean(); } /* -------------------------------------------------------------------------- */ FFTPlanManager & FFTPlanManager::get(){ if (FFTPlanManager::singleton == NULL) FFTPlanManager::singleton = new FFTPlanManager(); return *FFTPlanManager::singleton; } /* -------------------------------------------------------------------------- */ FFTransform & FFTPlanManager::createPlan(Grid & input, GridHermitian & output){ auto index = std::make_pair(const_cast(input.getInternalData()), const_cast(output.getInternalData())); auto it = plans.find(index); auto end = plans.end(); if (it == end) { #if defined(USE_CUFFT) plans[index] = new FFTransformCUFFT(input.size(),input,output); #else plans[index] = new FFTransformFFTW(input,output); #endif } return *plans[index]; } /* -------------------------------------------------------------------------- */ +void FFTPlanManager::clean() { + auto it = plans.begin(); + auto end = plans.end(); + for (; it != end ; ++it) { + delete it->second; + plans.erase(it->first); + } +} + +/* -------------------------------------------------------------------------- */ + template<> void FFTPlanManager::destroyPlan(const Real * some) { auto it = plans.begin(); auto end = plans.end(); for (; it != end ; ++it) { if (it->first.first == some) { delete it->second; plans.erase(it->first); } } } /* -------------------------------------------------------------------------- */ template<> void FFTPlanManager::destroyPlan(const complex * some) { auto it = plans.begin(); auto end = plans.end(); for (; it != end ; ++it) { if (it->first.second == some) { delete it->second; plans.erase(it->first); } } } /* -------------------------------------------------------------------------- */ /// FFT should not be used for any other type #define DESTROY(type) template<> void FFTPlanManager::destroyPlan( \ __attribute__((unused)) const type * some){} /// This list of types is defined in map_2d.cpp DESTROY(int); DESTROY(UInt); DESTROY(bool); DESTROY(unsigned long); #undef DESTROY /* -------------------------------------------------------------------------- */ void FFTPlanManager::destroyPlan(Grid & input, GridHermitian & output){ auto index = std::make_pair(const_cast(input.getInternalData()), const_cast(output.getInternalData())); auto it = plans.find(index); auto end = plans.end(); if (it != end){ delete (plans[index]); plans.erase(index); } } /* -------------------------------------------------------------------------- */ FFTPlanManager * FFTPlanManager::singleton = NULL; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/fft_plan_manager.hh b/src/fft_plan_manager.hh index ed8d35d..f7d6942 100644 --- a/src/fft_plan_manager.hh +++ b/src/fft_plan_manager.hh @@ -1,81 +1,84 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef FFT_PLAN_MANAGER_H #define FFT_PLAN_MANAGER_H /* -------------------------------------------------------------------------- */ #include #include "fftransform.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ class FFTPlanManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: FFTPlanManager(); public: virtual ~FFTPlanManager(); public: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ /// Get singleton instance static FFTPlanManager & get(); /// Create/retrieve a plan from two surfaces FFTransform & createPlan(Grid & input, GridHermitian & output); + /// Remove all plans + void clean(); + /// Destroy a plan from two surfaces void destroyPlan(Grid & input, GridHermitian & output); /// Destroy any plan containing a given data pointer template void destroyPlan(const T * some); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: std::map,FFTransform *> plans; static FFTPlanManager * singleton; }; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ /* -------------------------------------------------------------------------- */ #endif /* FFT_PLAN_MANAGER_H */ diff --git a/src/grid.cpp b/src/grid.cpp index bb96711..813b9ed 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -1,113 +1,114 @@ /** * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template Grid::Grid(): data() { memset(this->n, 0, dim * sizeof(UInt)); memset(this->L, 0, dim * sizeof(Real)); this->nb_components = 1; } /* -------------------------------------------------------------------------- */ template Grid::Grid(const UInt n[dim], const Real L[dim], UInt nb_components): data() { memcpy(this->n, n, dim * sizeof(UInt)); memcpy(this->L, L, dim * sizeof(Real)); this->nb_components = nb_components; data.resize(this->computeSize()); } /* -------------------------------------------------------------------------- */ template Grid::~Grid() {} /* -------------------------------------------------------------------------- */ template void Grid::resize(const UInt *n) { - memcpy(this->n, n, dim * sizeof(UInt)); + if (n != this->n) + memcpy(this->n, n, dim * sizeof(UInt)); UInt size = this->computeSize(); data.resize(size); data.assign(size, T(0.)); } /* -------------------------------------------------------------------------- */ template void Grid::printself(std::ostream & str) const { str << "Grid(" << dim << ", " << nb_components << ") {"; for (UInt i = 0 ; i < data.size() - 1 ; i++) { str << data[i] << ", "; } str << data[data.size()-1] << "}"; } /* -------------------------------------------------------------------------- */ #define GRID_SCALAR_OPERATOR_IMPL(op) \ template \ inline void Grid::operator op(const T & e) { \ _Pragma("omp parallel for") \ for (UInt i = 0 ; i < this->data.size() ; i++) { \ data[i] op e; \ } \ } \ GRID_SCALAR_OPERATOR_IMPL(+=); GRID_SCALAR_OPERATOR_IMPL(*=); GRID_SCALAR_OPERATOR_IMPL(-=); GRID_SCALAR_OPERATOR_IMPL(/=); GRID_SCALAR_OPERATOR_IMPL(=); #undef GRID_SCALAR_OPERATOR_IMPL /* -------------------------------------------------------------------------- */ /// Class instanciation #define GRID_INSTANCIATE_TYPE(type) template class Grid; \ template class Grid; \ template class Grid GRID_INSTANCIATE_TYPE(Real); GRID_INSTANCIATE_TYPE(UInt); GRID_INSTANCIATE_TYPE(complex); GRID_INSTANCIATE_TYPE(int); GRID_INSTANCIATE_TYPE(bool); GRID_INSTANCIATE_TYPE(unsigned long); #undef GRID_INSTANCIATE_TYPE __END_TAMAAS__ diff --git a/src/surface_generator_filter.cpp b/src/surface_generator_filter.cpp index 4dd4514..04173a2 100644 --- a/src/surface_generator_filter.cpp +++ b/src/surface_generator_filter.cpp @@ -1,147 +1,150 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "surface_generator_filter.hh" #include "fft_plan_manager.hh" #include "fftransform.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ SurfaceGeneratorFilter::SurfaceGeneratorFilter(){ h_coeff = NULL; sources = NULL; surface_spectral = NULL; white_noise_FFT = NULL; h_coeff_power = NULL; } /* -------------------------------------------------------------------------- */ SurfaceGeneratorFilter::~SurfaceGeneratorFilter(){ delete white_noise_FFT; delete surface_spectral; } /* -------------------------------------------------------------------------- */ void SurfaceGeneratorFilter::Init(){ SurfaceGenerator::Init(); } /* -------------------------------------------------------------------------- */ void SurfaceGeneratorFilter::applyFilterOnSource(){ std::cerr << "resize surface" << std::endl; if (surface_spectral == NULL) surface_spectral = new SurfaceComplex(surface->size(),surface->getL()); else surface_spectral->setGridSize(surface->size()); std::cerr << "do full product in fourier space" << std::endl; *surface_spectral = *white_noise_FFT; *surface_spectral *= *h_coeff_power; std::cerr << "fourier transform filter spectrum" << std::endl; FFTransform & plan = FFTPlanManager::get().createPlan(*surface, *surface_spectral); plan.backwardTransform(); + FFTPlanManager::get().destroyPlan(*surface, *surface_spectral); } /* -------------------------------------------------------------------------- */ void SurfaceGeneratorFilter::generateWhiteNoiseFFT(){ white_noise_FFT = new SurfaceComplex(sources->size(), sources->getL()); FFTransform & plan = FFTPlanManager::get().createPlan(*sources, *white_noise_FFT); plan.forwardTransform(); + FFTPlanManager::get().destroyPlan(*sources, *white_noise_FFT); // white_noise_FFT->dumpToParaview("white noise FFT",NO_REPLICATION); } /* -------------------------------------------------------------------------- */ void SurfaceGeneratorFilter::generateHFFT(){ FFTransform & plan = FFTPlanManager::get().createPlan(*h_coeff, *h_coeff_power); // fourier transform it plan.forwardTransform(); // make it real so that applicable filter is constructed h_coeff_power->makeItRealByAbs(); + FFTPlanManager::get().destroyPlan(*h_coeff, *h_coeff_power); } /* -------------------------------------------------------------------------- */ Surface & SurfaceGeneratorFilter::buildSurface(){ if (surface == NULL) SURFACE_FATAL("Init function was not called"); int n = surface->size(); if (white_noise_FFT == NULL){ // create the random source sources = new Surface(n,1.); Surface::generateWhiteNoiseSurface(*sources,random_seed); fprintf(stderr,"build FFT of white noise\n"); generateWhiteNoiseFFT(); delete (sources); } fprintf(stderr,"build filter coefficients\n"); computeFilterCoefficients(); if (h_coeff_power == NULL) { h_coeff_power = new SurfaceComplex(h_coeff->size(), h_coeff->getL()); generateHFFT(); } // s->setHeights(h_coeff_power); // h_coeff_power->dumpToParaview("h_coeff_power",NO_REPLICATION); fprintf(stderr,"apply filter on white noise\n"); applyFilterOnSource(); delete h_coeff_power; fprintf(stderr,"recenter arround mean the surface\n"); surface->recenterHeights(); constrainRMS(); // s->setHeights(heights); //s->SetHeights(sources); //s->SetHeights(h); return *surface; } /* -------------------------------------------------------------------------- */ void SurfaceGeneratorFilter::clearSource(){ if (white_noise_FFT != NULL){ delete white_noise_FFT; white_noise_FFT = NULL; } } __END_TAMAAS__ diff --git a/src/surface_statistics.hh b/src/surface_statistics.hh index 03d3bae..de20989 100644 --- a/src/surface_statistics.hh +++ b/src/surface_statistics.hh @@ -1,722 +1,721 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __SURFACE_STATISTICS_HH__ #define __SURFACE_STATISTICS_HH__ /* -------------------------------------------------------------------------- */ #include "surface.hh" #include "fftransform.hh" #include "fft_plan_manager.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class SurfaceStatistics { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SurfaceStatistics(){}; virtual ~SurfaceStatistics(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline static Real computeMaximum(const Surface & s); inline static Real computeMinimum(const Surface & s); inline static UInt computeAverage(const Surface & s, UInt epsilon=0){SURFACE_FATAL("TO IMPLEMENT");}; inline static unsigned long computeAverage(const Surface & s, unsigned long epsilon=0){SURFACE_FATAL("TO IMPLEMENT");}; inline static Real computeAverage(const Surface & s, Real epsilon= -1.*std::numeric_limits::max()){ return computeAverage(s,epsilon); } template inline static T computeAverage(const Surface & s, T epsilon= -1.*std::numeric_limits::max()); inline static Real computeStdev(const Surface & s, Real epsilon=-1. * std::numeric_limits::max()); inline static Real computeSkewness(const Surface & surface, Real avg, Real std); inline static Real computeKurtosis(const Surface & surface, Real avg, Real std); inline static Real computeRMSSlope(const Surface & surface); inline static Real computeSpectralRMSSlope(Surface & surface); inline static Real computeSpectralMeanCurvature(Surface & surface); inline static Real computeSpectralStdev(Surface & surface); // inline static void computeStatistics(Surface & surface); inline static std::vector computeMoments(Surface & surface); inline static Real computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst, Real rms, Real L); //! compute perimeter inline static Real computePerimeter(const Surface & surface, Real eps = 1e-10); //! get real contact area inline static Real computeContactArea(const Surface & surface); //! get real contact area inline static Real computeContactAreaRatio(const Surface & surface); template inline static std::vector computeDistribution(const Surface & surface, unsigned int n_bins, Real min_value=std::numeric_limits::max(), Real max_value=-1*std::numeric_limits::max()); inline static std::vector computeSpectralDistribution(const Surface & surface, unsigned int n_bins, UInt cutoff=std::numeric_limits::max()); //! compute average amplitude inline static Real computeAverageAmplitude(const Surface & surface); inline static Real computeAverageAmplitude(const SurfaceComplex & surface); //! sum all heights inline static Real computeSum(const Surface & surface); //! compute auto correlation inline static Surface computeAutocorrelation(Surface & surface); //! compute powerspectrum inline static SurfaceComplex computePowerSpectrum(Surface & surface); }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "surface_statistics_inline_impl.cc" /* -------------------------------------------------------------------------- */ template std::vector SurfaceStatistics::computeDistribution(const Surface & surface, unsigned int n_bins, Real min_value,Real max_value){ if (n_bins == 0) SURFACE_FATAL("badly set number of bins"); Real min_heights; Real max_heights; std::vector bins; bins.resize(n_bins); bins.assign(n_bins,0); min_heights = min_value; max_heights = max_value; if (max_value == -1.*std::numeric_limits::max() ) max_heights = computeMaximum(surface); if (min_value == std::numeric_limits::max() ) min_heights = computeMinimum(surface); int counter = 0; const int n = surface.size(); for (int i = 0 ; i < n*n ; ++i) { Real value = surface(i); if (!consider_zeros && value < zero_threshold) continue; if (value > max_heights || value < min_heights) SURFACE_FATAL("out of bound value: " << min_heights << " ::max(); #pragma omp parallel for reduction(min: min_non_zero) for (int i = 0 ; i < n*n ; ++i) { Real value = surface(i); if (value < min_non_zero && value > 1e-14) min_non_zero = value; } return bins; } /* -------------------------------------------------------------------------- */ std::vector SurfaceStatistics::computeSpectralDistribution(const Surface & surface, unsigned int n_bins, UInt cutoff){ if (n_bins == 0) SURFACE_FATAL("badly set number of bins"); const int n = surface.size(); Real min_value = SurfaceStatistics::computeMinimum(surface); Real max_value = SurfaceStatistics::computeMaximum(surface); std::vector distrib_fft(n_bins); const int N = n_bins; Real L = max_value - min_value; #pragma omp parallel for for (int k = 0 ; k < N ; ++k) distrib_fft[k] = 0.; #pragma omp parallel for for (int k = 1 ; k < N/2 ; ++k) { for (int i = 0 ; i < n*n ; ++i) { Real value = surface(i); Real angle = -2.*M_PI*k*(value-min_value)/L; complex phase = std::polar(1.,angle); distrib_fft[k] += phase; distrib_fft[N-k] += phase; } } for (int i = 0 ; i < n*n ; ++i) { Real value = surface(i); Real angle = -2.*M_PI*N/2*(value-min_value)/L; complex phase = std::polar(1.,angle); distrib_fft[N/2] += phase; } #pragma omp parallel for for (int k = 1 ; k < N ; ++k) distrib_fft[k] /= n*n; distrib_fft[0] = std::polar(1.,0.); #pragma omp parallel for for (int k = 1 ; k < N/2 ; ++k) { if (k >= (int)cutoff) { distrib_fft[k] = 0.; distrib_fft[N-k] = 0.; } } //proceed to inverse fourier transform #pragma omp parallel for for (int k = 0 ; k < N ; ++k) { distrib_fft[k] = std::conj(distrib_fft[k]); } fftw_complex * in = (fftw_complex*)&distrib_fft[0]; fftw_plan p = fftw_plan_dft_1d(N, in, in ,FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(p); fftw_destroy_plan(p); std::vector distrib; distrib.resize(N); #pragma omp parallel for for (int k = 0 ; k < N ; ++k) { distrib[k] = std::conj(distrib_fft[k]); } return distrib; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeAverageAmplitude(const Surface & surface){ Real average_amplitude = 0; - UInt n = surface.size(); -#pragma omp parallel for collapse(2), reduction(+: average_amplitude) + UInt n = surface.dataSize(); +#pragma omp parallel for reduction(+: average_amplitude) for (UInt i = 0 ; i < n ; ++i) - for (UInt j = 0 ; j < n ; ++j) - average_amplitude += std::abs(surface(i,j)); + average_amplitude += std::abs(surface(i)); - average_amplitude *= 1./n/n; + average_amplitude *= 1./n; return average_amplitude; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeAverageAmplitude(const SurfaceComplex & surface){ Real average_amplitude = 0; - UInt n = surface.size(); -#pragma omp parallel for collapse(2), reduction(+: average_amplitude) + UInt n = surface.dataSize(); +#pragma omp parallel for reduction(+: average_amplitude) for (UInt i = 0 ; i < n ; ++i) - for (UInt j = 0 ; j < n ; ++j) - average_amplitude += std::abs(surface(i,j)); + average_amplitude += std::abs(surface(i)); - average_amplitude *= 1./n/n; + average_amplitude *= 1./n; return average_amplitude; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeContactArea(const Surface & surface) { Real L = surface.getL(); return SurfaceStatistics::computeContactAreaRatio(surface)*L*L; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeContactAreaRatio(const Surface & surface) { UInt area = 0; UInt n = surface.size(); #pragma omp parallel for collapse(2), reduction(+:area) for (UInt i = 0; i 1e-10) area++; } } return area/Real(n*n); } /* -------------------------------------------------------------------------- */ Surface SurfaceStatistics::computeAutocorrelation(Surface & surface){ SurfaceComplex correlation_heights = computePowerSpectrum(surface); Surface acf(correlation_heights.size(), correlation_heights.getL()); FFTransform & plan = FFTPlanManager::get().createPlan(acf, correlation_heights); plan.backwardTransform(); + FFTPlanManager::get().destroyPlan(acf, correlation_heights); return acf; } /* -------------------------------------------------------------------------- */ SurfaceComplex SurfaceStatistics::computePowerSpectrum(Surface & surface){ SurfaceComplex power_spectrum(surface.size(), surface.getL()); FFTransform & plan = FFTPlanManager::get().createPlan(surface, power_spectrum); plan.forwardTransform(); - UInt n = power_spectrum.size(); - Real factor = 1./(n*n); + UInt n = power_spectrum.dataSize(); + UInt factor = power_spectrum.size()*power_spectrum.size(); -#pragma omp parallel for collapse(2) +#pragma omp parallel for for (UInt i = 0 ; i < n ; ++i) - for (UInt j = 0 ; j < n ; ++j) - power_spectrum(i,j) *= factor; + power_spectrum(i) /= factor; power_spectrum.makeItRealBySquare(); + FFTPlanManager::get().destroyPlan(surface, power_spectrum); return power_spectrum; } /* -------------------------------------------------------------------------- */ template T SurfaceStatistics::computeAverage(const Surface & surface, T epsilon){ Real avg = 0.0; UInt cpt = 0; UInt n = surface.size(); #pragma omp parallel for reduction(+: avg, cpt) for (UInt i = 0; i < n*n; i++) { Real val = surface(i); if (std::abs(val) > epsilon){ avg += val; ++cpt; } } avg /= cpt; return avg; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeSkewness(const Surface & surface, Real avg, Real std){ Real skw = 0.0; UInt n = surface.size(); #pragma omp parallel for reduction(+: skw) for (UInt i = 0; i < n*n; i++) { Real tmp = surface(i) - avg; skw += (tmp*tmp*tmp); } skw /= (n*std*std*std); return skw; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeKurtosis(const Surface & surface, Real avg, Real std){ Real kur=0.0; UInt n = surface.size(); #pragma omp parallel for reduction(+: kur) for (UInt i = 0; i < n*n; i++) { Real tmp = surface.size() - avg; kur += (tmp*tmp*tmp*tmp); } kur /= (n*std*std*std*std); return kur-3; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeMaximum(const Surface & surface){ Real _max = -1. * std::numeric_limits::max(); UInt n = surface.size(); #pragma omp parallel for reduction(max: _max) for (UInt i = 0; i < n*n; i++) { _max = std::max(surface(i),_max); } return _max; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeMinimum(const Surface & surface){ Real _min = std::numeric_limits::max(); UInt n = surface.size(); #pragma omp parallel for reduction(min: _min) for (UInt i = 0; i < n*n; i++) { _min = std::min(surface(i),_min); } return _min; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeRMSSlope(const Surface & surface){ UInt n = surface.size(); Surface grad(n,1.); for (UInt i = 0, m = 0; i < n; i++) { for (UInt j = 0; j < n; j++) { UInt k = i+n*j; // xl = ((i==0)?(k+n-1):(k-1)); UInt xh = ((i==(n-1))?(k-n+1):(k+1)); // yl = ((j==0)?(k+n*n-n):(k-n)); UInt yh = ((j==(n-1))?(k-n*n+n):(k+n)); Real gx = surface(xh) - surface(k); Real gy = surface(yh) - surface(k); grad(m++) = gx*gx + gy*gy; } } Real avg = SurfaceStatistics::computeAverage(grad); Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */ Real scale_factor = surface.getL()/n; return slp/scale_factor; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeSpectralRMSSlope(Surface & surface){ SurfaceComplex power = computePowerSpectrum(surface); UInt n = surface.size(); Real rms_slopes = 0; // /!\ commented code was causing race conditions with dummy variable in GridHermitian #pragma omp parallel for collapse(2), reduction(+:rms_slopes) for (UInt i = 1 ; i < n/2 ; ++i) for (UInt j = 1 ; j < n/2 ; ++j){ rms_slopes += 2.*((i*i)+(j*j))*power(i,j).real(); //rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real(); rms_slopes += 2.*((i*i)+(j*j))*power(n-i,j).real(); //rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real(); } // external line (i or j == 0) #pragma omp parallel for reduction(+:rms_slopes) for (UInt i = 1 ; i < n/2 ; ++i){ rms_slopes += 1.*(i*i)*power(i,0).real(); rms_slopes += 1.*(i*i)*power(n-i,0).real(); rms_slopes += 2.*(i*i)*power(0,i).real(); //rms_slopes += 1.*(i*i)*power(0,n-i).real(); } //internal line (i or j == n/2) #pragma omp parallel for reduction(+:rms_slopes) for (UInt i = 1; i < n/2; ++i) { rms_slopes += 0.5*(i*i+n*n/4)*power(i,n/2).real(); rms_slopes += 0.5*(i*i+n*n/4)*power(n-i,n/2).real(); rms_slopes += 1.*(i*i+n*n/4)*power(n/2,i).real(); //rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real(); } { //i == n/2 j == n/2 rms_slopes += 0.25*(n*n/2)*power(n/2,n/2).real(); //i == 0 j == n/2 rms_slopes += 1.*(n*n/4)*power(0,n/2).real(); //i == n/2 j == 0 rms_slopes += 1.*(n*n/4)*power(n/2,0).real(); } rms_slopes = 2.*M_PI/n*sqrt(rms_slopes); Real scale_factor = surface.getL()/n; return rms_slopes/scale_factor; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeStdev(const Surface & surface, Real epsilon){ UInt n = surface.size(); UInt cpt = 0; Real std = 0.0; Real avg = 0.0; #pragma omp parallel for reduction(+:avg, std, cpt) for (UInt i = 0; i < n*n; i++) { Real val = surface(i); if (val > epsilon){ avg += val; std += val*val; ++cpt; } } avg /= cpt; return sqrt(std/cpt - avg*avg); } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeSpectralStdev(Surface & surface){ SurfaceComplex power = computePowerSpectrum(surface); Real rms_heights = 0; UInt n = surface.size(); #pragma omp parallel for collapse(2), reduction(+:rms_heights) for (UInt i = 0 ; i < n ; ++i) for (UInt j = 0 ; j < n ; ++j){ rms_heights += power(i,j).real(); } rms_heights = sqrt(rms_heights); return rms_heights; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeSpectralMeanCurvature(Surface & surface){ SurfaceComplex power = computePowerSpectrum(surface); UInt n = surface.size(); Real rms_curvatures = 0; #pragma omp parallel for collapse(2), reduction(+:rms_curvatures) for (UInt i = 1 ; i < n/2 ; ++i) for (UInt j = 1 ; j < n/2 ; ++j){ Real coeff = ((i*i)+(j*j)); coeff *= coeff; rms_curvatures += 1.*coeff*power(i,j).real(); rms_curvatures += 1.*coeff*power(n-i,n-j).real(); rms_curvatures += 1.*coeff*power(n-i,j).real(); rms_curvatures += 1.*coeff*power(i,n-j).real(); } // external line (i or j == 0) #pragma omp parallel for reduction(+:rms_curvatures) for (UInt i = 1 ; i < n/2 ; ++i){ Real coeff = i*i; coeff *= coeff; rms_curvatures += 4./3.*coeff*power(i,0).real(); rms_curvatures += 4./3.*coeff*power(n-i,0).real(); rms_curvatures += 4./3.*coeff*power(0,i).real(); rms_curvatures += 4./3.*coeff*power(0,n-i).real(); } //internal line (i or j == n/2) #pragma omp parallel for reduction(+:rms_curvatures) for (UInt i = 1; i < n/2; ++i) { Real coeff = (i*i+n*n/4); coeff *= coeff; rms_curvatures += 0.5*coeff*power(i,n/2).real(); rms_curvatures += 0.5*coeff*power(n-i,n/2).real(); rms_curvatures += 0.5*coeff*power(n/2,i).real(); rms_curvatures += 0.5*coeff*power(n/2,n-i).real(); } { //i == n/2 j == n/2 rms_curvatures += 0.25*(n*n/2)*(n*n/2)*power(n/2,n/2).real(); //i == 0 j == n/2 rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(0,n/2).real(); //i == n/2 j == 0 rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(n/2,0).real(); } rms_curvatures *= 9./4.; rms_curvatures = M_PI/n/n * sqrt(rms_curvatures); return rms_curvatures; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computePerimeter(const Surface & surface, Real eps) { int perimeter = 0; UInt n = surface.size(); #pragma omp parallel for collapse(2), reduction(+:perimeter) for (UInt i = 0 ; i < n ; ++i) { for (UInt j = 0 ; j < n ; ++j) { if ( j > 0 && ( ( surface(i,(j-1)) < eps && surface(i,j) > eps ) || ( surface(i,(j-1)) > eps && surface(i,j) < eps ) ) ) perimeter++; if ( i > 0 && ( ( surface(j,(i-1)) < eps && surface(j,i) > eps ) || ( surface(j,(i-1)) > eps && surface(j,i) < eps )) ) perimeter++; } } return perimeter/Real(n); } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeSum(const Surface & surface){ Real sum = 0.; UInt n = surface.size(); #pragma omp parallel for reduction(+: sum) for (UInt i = 0 ; i < n*n ; ++i) sum += surface(i); return sum; } /* -------------------------------------------------------------------------- */ std::vector SurfaceStatistics::computeMoments(Surface & surface){ SurfaceComplex power = computePowerSpectrum(surface); Real L = surface.getL(); UInt n = surface.size(); Real k_factor = 2.*M_PI/L; std::vector > moment; moment.resize(6); for (int i = 0; i < 5; i++) moment[i].resize(5); for (UInt i = 0; i < n/2; i++) { Real i0 = pow(Real(i)*k_factor,0); Real i2 = pow(Real(i)*k_factor,2); Real i4 = pow(Real(i)*k_factor,4); Real im0 = pow(-Real(i)*k_factor,0); Real im2 = pow(-Real(i)*k_factor,2); Real im4 = pow(-Real(i)*k_factor,4); for (UInt j = 0; j < n/2; j++) { Real Phi = power(i,j).real(); Real Phi_10; Real Phi_01; Real Phi_11; if (i==0 && j==0) { Phi_10 = 0; Phi_11 = 0; Phi_01 = 0; } else if (i==0) { Phi_10 = 0; Phi_11 = 0; Phi_01 = power(i,(n-j)).real(); } else if (j==0) { Phi_10 = power(n-i,j).real(); Phi_11 = 0; Phi_01 = 0; } else { Phi_10 = power(n-i,j).real(); Phi_11 = power(n-i,(n-j)).real(); Phi_01 = power(i,(n-j)).real(); } Real j0 = pow(Real(j)*k_factor,0); Real j2 = pow(Real(j)*k_factor,2); Real j4 = pow(Real(j)*k_factor,4); Real jm0 = pow(-Real(j)*k_factor,0); Real jm2 = pow(-Real(j)*k_factor,2); Real jm4 = pow(-Real(j)*k_factor,4); moment[0][0] += i0*j0*Phi + im0*j0*Phi_10 + i0*jm0*Phi_01 + im0*jm0*Phi_11; moment[0][2] += i0*j2*Phi + im0*j2*Phi_10 + i0*jm2*Phi_01 + im0*jm2*Phi_11; moment[2][0] += i2*j0*Phi + im2*j0*Phi_10 + i2*jm0*Phi_01 + im2*jm0*Phi_11; moment[2][2] += i2*j2*Phi + im2*j2*Phi_10 + i2*jm2*Phi_01 + im2*jm2*Phi_11; moment[0][4] += i0*j4*Phi + im0*j4*Phi_10 + i0*jm4*Phi_01 + im0*jm4*Phi_11; moment[4][0] += i4*j0*Phi + im4*j0*Phi_10 + i4*jm0*Phi_01 + im4*jm0*Phi_11; } } std::vector mean_moments; //m2 Real mean_m2 = 0.5*(moment[0][2]+moment[2][0]); mean_moments.push_back(mean_m2); //m4 Real mean_m4 = (3*moment[2][2]+moment[0][4]+moment[4][0])/3.; mean_moments.push_back(mean_m4); return mean_moments; } /* -------------------------------------------------------------------------- */ Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst,Real A, Real L){ Real xi = k2/k1; Real exp = - 2.*hurst; Real _k1 = 2.*M_PI/L*k1; Real Cq; switch (order){ case 0: Cq = 2.*M_PI;break; case 2: Cq = M_PI;break; case 4: Cq = 3./4.*M_PI;break; default: SURFACE_FATAL("invalid order"); } exp = order - 2.*hurst; Real res = A * Cq * pow(_k1,exp) / exp * ( pow(xi, exp) - 1 ); return res; } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ #endif /* __SURFACE_STATISTICS_HH__ */