Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/core/surface_statistics.hh
cb641f5bb816features/netcdf-mpi

Last Change
Raw File

surface_statistics.hh
View Options

/**
* @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 <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef SURFACE_STATISTICS_HH
#define SURFACE_STATISTICS_HH
/* -------------------------------------------------------------------------- */
#include "surface.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @deprecated should move code to Statistics class instead
* @brief Statistics calculations on surfaces
*/
class SurfaceStatistics {
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
static Real computeMaximum(const Surface<Real>& s);
static Real computeMinimum(const Surface<Real>& s);
static Real
computeAverage(const Surface<Real>& s,
Real epsilon = std::numeric_limits<Real>::lowest()) {
return computeAverage<Real>(s, epsilon);
}
template <typename T>
static T computeAverage(const Surface<T>& s,
T epsilon = std::numeric_limits<T>::lowest());
static Real
computeStdev(const Grid<Real, 2>& s,
Real epsilon = std::numeric_limits<Real>::lowest());
static Real computeSkewness(const Surface<Real>& surface, Real avg,
Real std);
static Real computeKurtosis(const Surface<Real>& surface, Real avg,
Real std);
static Real computeRMSSlope(const Surface<Real>& surface);
static Real computeSpectralRMSSlope(Surface<Real>& surface);
static Real computeSpectralMeanCurvature(Surface<Real>& surface);
static Real computeSpectralStdev(Surface<Real>& surface);
// static void computeStatistics(Surface<Real> & surface);
static std::vector<Real> computeMoments(Surface<Real>& surface);
static Real computeAnalyticFractalMoment(UInt order, UInt k1, UInt k2,
Real hurst, Real rms, Real L);
//! compute perimeter
static Real computePerimeter(const Surface<Real>& surface,
Real eps = 1e-10);
//! get real contact area
static Real computeContactArea(const Surface<Real>& surface);
//! get real contact area
static Real computeContactAreaRatio(const Surface<Real>& surface);
template <bool consider_zeros = false>
static std::vector<int>
computeDistribution(const Surface<Real>& surface, unsigned int n_bins,
Real min_value = std::numeric_limits<Real>::max(),
Real max_value = std::numeric_limits<Real>::lowest());
static std::vector<Complex>
computeSpectralDistribution(const Surface<Real>& surface, unsigned int n_bins,
UInt cutoff = std::numeric_limits<UInt>::max());
//! compute average amplitude
static Real computeAverageAmplitude(const Surface<Real>& surface);
static Complex
computeAverageAmplitude(const SurfaceComplex<Real>& surface);
//! sum all heights
static Real computeSum(const Surface<Real>& surface);
//! compute auto correlation
static Surface<Real> computeAutocorrelation(Surface<Real>& surface);
//! compute powerspectrum
static SurfaceComplex<Real>
computePowerSpectrum(Surface<Real>& surface);
};
/* -------------------------------------------------------------------------- */
template <bool consider_zeros>
std::vector<int>
SurfaceStatistics::computeDistribution(const Surface<Real>& 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<int> 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<Real>::max())
max_heights = computeMaximum(surface);
if (min_value == std::numeric_limits<Real>::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 < 1e-14)
continue;
if (value > max_heights || value < min_heights)
SURFACE_FATAL("out of bound value: "
<< min_heights << " <? " << value << " <? " << max_heights
<< " : bounds probably badly chosen");
Real ind = (value - min_heights) / (max_heights - min_heights);
unsigned int index = ind * n_bins;
if (index == n_bins)
--index;
++bins[index];
++counter;
}
if (consider_zeros && counter != n * n)
SURFACE_FATAL("problem in counting points");
// Real offset = (max_heights - min_heights)/n_bins;
// Real x = 0;
// Real f1 = 0, f2 = 0;
// for (unsigned int i = 0 ; i < n_bins ; ++i)
// {
// x = i*offset+min_heights;
// f1 = bins[i]/Real(n*n);
// f2 = bins[i]/Real(counter);
//}
Real min_non_zero = std::numeric_limits<Real>::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;
}
/* -------------------------------------------------------------------------- */
template <typename T>
T SurfaceStatistics::computeAverage(const Surface<T>& 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;
}
} // namespace tamaas
#endif /* __SURFACE_STATISTICS_HH__ */
c4science · Help