diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp index 5567cab..90784b6 100644 --- a/src/surface/isopowerlaw.cpp +++ b/src/surface/isopowerlaw.cpp @@ -1,102 +1,101 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "isopowerlaw.hh" -#include "fftransform.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { template void Isopowerlaw::computeFilter( GridHermitian& filter_coefficients) const { Filter::computeFilter( [this] CUDA_LAMBDA(Complex & coeff, VectorProxy q) { coeff = (*this)(q); }, filter_coefficients); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::rmsHeights() const { return std::sqrt(moments()[0]); } /* -------------------------------------------------------------------------- */ /** * 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 Isopowerlaw<2>::moments() const { std::map T; T[0] = 2 * M_PI; T[2] = M_PI; T[4] = 3 * M_PI / 4.; Real q1 = static_cast(this->q1); // shadowing this->q1 Real xi = q0 / q1; Real zeta = q2 / q1; Real A = std::pow(q1, 2 * (hurst + 1)); std::vector moments; moments.reserve(3); for (UInt q : {0, 2, 4}) { Real m = A * T[q] * std::pow(q1, q - 2. * hurst) * ((1 - std::pow(xi, q + 2)) / (q + 2.) + (std::pow(zeta, q - 2. * hurst) - 1) / (q - 2. * hurst)); moments.push_back(m); } return moments; } template <> std::vector Isopowerlaw<1>::moments() const { TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces"); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::alpha() const { auto m = moments(); return m[0] * m[2] / (m[1] * m[1]); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::rmsSlopes() const { // 2 pi because moments are computed with k instead of q // and slopes should be computed with q return 2 * M_PI * std::sqrt(2 * moments()[1]); } template class Isopowerlaw<1>; template class Isopowerlaw<2>; } // namespace tamaas diff --git a/src/surface/isopowerlaw.hh b/src/surface/isopowerlaw.hh index 6d0ab6d..b12ebd8 100644 --- a/src/surface/isopowerlaw.hh +++ b/src/surface/isopowerlaw.hh @@ -1,95 +1,94 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef ISOPOWERLAW_H #define ISOPOWERLAW_H /* -------------------------------------------------------------------------- */ #include "filter.hh" #include "grid_hermitian.hh" #include "static_types.hh" -#include "tamaas.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /// Class representing an isotropic power law spectrum template class Isopowerlaw : public Filter { public: /// Compute filter coefficients void computeFilter(GridHermitian& filter_coefficients) const override; /// Compute a point of the PSD CUDA_LAMBDA inline Real operator()(const VectorProxy& q_vec) const; /// Analytical rms of heights Real rmsHeights() const; /// Analytical moments std::vector 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 CUDA_LAMBDA inline Real Isopowerlaw:: operator()(const VectorProxy& 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 + dim)); return std::sqrt(val); } } // namespace tamaas #endif // ISOPOWERLAW_H diff --git a/src/surface/regularized_powerlaw.hh b/src/surface/regularized_powerlaw.hh index 44b677b..a3723ef 100644 --- a/src/surface/regularized_powerlaw.hh +++ b/src/surface/regularized_powerlaw.hh @@ -1,78 +1,77 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef REGULARIZED_POWERLAW_H #define REGULARIZED_POWERLAW_H /* -------------------------------------------------------------------------- */ #include "filter.hh" #include "grid_hermitian.hh" #include "static_types.hh" -#include "tamaas.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /// Class representing an isotropic power law spectrum template class RegularizedPowerlaw : public Filter { public: /// Compute filter coefficients void computeFilter(GridHermitian& filter_coefficients) const override; /// Compute a point of the PSD CUDA_LAMBDA inline Real operator()(const VectorProxy& q_vec) const; public: TAMAAS_ACCESSOR(q1, UInt, Q1); TAMAAS_ACCESSOR(q2, UInt, Q2); TAMAAS_ACCESSOR(hurst, Real, Hurst); protected: UInt q1, q2; Real hurst; }; /* -------------------------------------------------------------------------- */ /* Inline implementations */ /* -------------------------------------------------------------------------- */ template CUDA_LAMBDA inline Real RegularizedPowerlaw:: operator()(const VectorProxy& q_vec) const { const Real C = 1.; const Real q = q_vec.l2norm(); Real val; if (q > q2) val = 0; else val = C * std::pow(1 + (q / q1) * (q / q1), -(hurst + 1)); return std::sqrt(val); } } // namespace tamaas #endif // REGULARIZED_POWERLAW_H diff --git a/src/surface/surface_generator_filter.cpp b/src/surface/surface_generator_filter.cpp index f03392b..ab3496e 100644 --- a/src/surface/surface_generator_filter.cpp +++ b/src/surface/surface_generator_filter.cpp @@ -1,84 +1,73 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "surface_generator_filter.hh" -#include "fft_plan_manager.hh" -#include "fftransform.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ template void SurfaceGeneratorFilter::applyFilterOnSource() { - GridHermitian fft_white_noise(this->filter_coefficients.sizes(), - 1); - FFTransform& wn_plan = - FFTPlanManager::get().createPlan(this->white_noise, fft_white_noise); - wn_plan.forwardTransform(); + GridHermitian fft_white_noise(filter_coefficients.sizes(), 1); + engine.forward(this->white_noise, fft_white_noise); fft_white_noise *= filter_coefficients; - - FFTransform& surface_plan = - FFTPlanManager::get().createPlan(this->grid, fft_white_noise); - surface_plan.backwardTransform(); - - FFTPlanManager::get().destroyPlan(this->white_noise, fft_white_noise); - FFTPlanManager::get().destroyPlan(this->grid, fft_white_noise); + engine.backward(this->grid, fft_white_noise); } /* -------------------------------------------------------------------------- */ template Grid& SurfaceGeneratorFilter::buildSurface() { if (this->grid.dataSize() == 0) TAMAAS_EXCEPTION("the size of the grid is zero, did you call setSizes() ?"); if (this->filter == nullptr) TAMAAS_EXCEPTION("filter is null, did you call setFilter() ?"); // Resizing arrays this->white_noise.resize(this->grid.sizes()); auto hermitian_dim = GridHermitian::hermitianDimensions(this->grid.sizes()); this->filter_coefficients.resize(hermitian_dim); /// Generating white noise with gaussian distribution this->generateWhiteNoise>(); /// Applying filters filter->computeFilter(this->filter_coefficients); this->applyFilterOnSource(); /// Normalizing surface for hrms independent of number of points this->grid *= std::sqrt(this->grid.dataSize()); return this->grid; } /* -------------------------------------------------------------------------- */ template class SurfaceGeneratorFilter<1>; template class SurfaceGeneratorFilter<2>; /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/surface/surface_generator_filter.hh b/src/surface/surface_generator_filter.hh index 2b36967..32f2f85 100644 --- a/src/surface/surface_generator_filter.hh +++ b/src/surface/surface_generator_filter.hh @@ -1,79 +1,81 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef SURFACE_GENERATOR_FILTER_H #define SURFACE_GENERATOR_FILTER_H /* -------------------------------------------------------------------------- */ #include "filter.hh" +#include "fft_engine.hh" #include "surface_generator.hh" #include "tamaas.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { template class SurfaceGeneratorFilter : public SurfaceGenerator { public: /// Default constructor SurfaceGeneratorFilter() = default; public: /// Build surface with Hu & Tonder algorithm Grid& buildSurface() override; /// Set filter object void setFilter(std::shared_ptr> new_filter) { setSpectrum(new_filter); } /// Set spectrum void setSpectrum(std::shared_ptr> spectrum) { filter = spectrum; } protected: /// Apply filter coefficients on white noise void applyFilterOnSource(); /// Generate white noise with given distribution template void generateWhiteNoise(); protected: std::shared_ptr> filter = nullptr; ///< not owned GridHermitian filter_coefficients; Grid white_noise; + FFTEngine engine; }; /* -------------------------------------------------------------------------- */ /* Template implementations */ /* -------------------------------------------------------------------------- */ template template void SurfaceGeneratorFilter::generateWhiteNoise() { random_engine gen(this->random_seed); T distribution; for (auto& noise : white_noise) noise = distribution(gen); } } // namespace tamaas #endif diff --git a/src/surface/surface_generator_filter_fft.cpp b/src/surface/surface_generator_filter_fft.cpp index 1765439..e69d4b3 100644 --- a/src/surface/surface_generator_filter_fft.cpp +++ b/src/surface/surface_generator_filter_fft.cpp @@ -1,55 +1,55 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "surface_generator_filter_fft.hh" #include "isopowerlaw.hh" #include "logger.hh" -#include "surface_statistics.hh" +#include "statistics.hh" +#include /* -------------------------------------------------------------------------- */ namespace tamaas { extern template class Grid; SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT() { pl_filter = std::make_shared>(); generator.setFilter(pl_filter); } Real SurfaceGeneratorFilterFFT::constrainRMS() { Logger().get(LogLevel::warning) << "Statistical bias in using the generated RMS as " - "normalizing factor" - << std::endl; - Real hrms_grid = SurfaceStatistics::computeStdev(generator.grid); + "normalizing factor!\n"; + Real hrms_grid = Statistics<2>::computeRMSHeights(generator.grid); generator.grid *= rms / hrms_grid; return rms; } Grid& SurfaceGeneratorFilterFFT::buildSurface() { generator.setSizes({{grid_size, grid_size}}); return generator.buildSurface(); } } // namespace tamaas diff --git a/src/surface/surface_generator_filter_fft.hh b/src/surface/surface_generator_filter_fft.hh index 0b02e53..d730ddd 100644 --- a/src/surface/surface_generator_filter_fft.hh +++ b/src/surface/surface_generator_filter_fft.hh @@ -1,73 +1,74 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef SURFACE_GENERATOR_FILTER_FFT_H #define SURFACE_GENERATOR_FILTER_FFT_H /* -------------------------------------------------------------------------- */ #include "isopowerlaw.hh" #include "surface_generator_filter.hh" +#include /* -------------------------------------------------------------------------- */ namespace tamaas { /// Legacy class to maintain interface class SurfaceGeneratorFilterFFT { public: //! can construct a generator with default values 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->rmsHeights(); } /// Swig does not wrap the numpy otherwise Grid& 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; std::shared_ptr> pl_filter = nullptr; UInt grid_size; Real rms; }; } // namespace tamaas #endif diff --git a/src/surface/surface_generator_random_phase.cpp b/src/surface/surface_generator_random_phase.cpp index 137b1a1..b9af274 100644 --- a/src/surface/surface_generator_random_phase.cpp +++ b/src/surface/surface_generator_random_phase.cpp @@ -1,72 +1,67 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "surface_generator_random_phase.hh" -#include "fft_plan_manager.hh" -#include "fftransform.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ template Grid& SurfaceGeneratorRandomPhase::buildSurface() { if (this->grid.dataSize() == 0) TAMAAS_EXCEPTION("the size of the grid is zero, did you call setSizes() ?"); if (this->filter == nullptr) TAMAAS_EXCEPTION("filter is null, did you call setFilter() ?"); // Resizing arrays auto hermitian_dim = GridHermitian::hermitianDimensions(this->grid.sizes()); this->white_noise.resize(hermitian_dim); this->filter_coefficients.resize(hermitian_dim); /// Generating white noise with uniform distribution for phase this->template generateWhiteNoise>(); /// Applying filters this->filter->computeFilter(this->filter_coefficients); Loop::loop([](Complex& coeff, Real& phase) { coeff *= thrust::polar(1., 2 * M_PI * phase); }, this->filter_coefficients, this->white_noise); - FFTPlanManager::get() - .createPlan(this->grid, this->filter_coefficients) - .backwardTransform(); - FFTPlanManager::get().destroyPlan(this->grid, this->filter_coefficients); + this->engine.backward(this->grid, this->filter_coefficients); /// Normalizing surface for hrms independent of number of points this->grid *= this->grid.dataSize(); return this->grid; } /* -------------------------------------------------------------------------- */ template class SurfaceGeneratorRandomPhase<1>; template class SurfaceGeneratorRandomPhase<2>; /* -------------------------------------------------------------------------- */ } // namespace tamaas