diff --git a/.clang-tidy b/.clang-tidy
index 76c907d..d2fdbe2 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -1,4 +1,4 @@
-Checks: 'modernize-use-*, performance-*, cppcoreguidelines-*'
+Checks: '-*, modernize-use-*, -modernize-use-trailing-return-type*, performance-*, readability-*, -readability-braces-around-statements, -readability-implicit-bool-conversion'
AnalyzeTemporaryDtors: false
HeaderFilterRegex: 'src/.*'
FormatStyle: file
diff --git a/src/core/fftw/fftw_engine.hh b/src/core/fftw/fftw_engine.hh
index e5a2634..3273aba 100644
--- a/src/core/fftw/fftw_engine.hh
+++ b/src/core/fftw/fftw_engine.hh
@@ -1,104 +1,104 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 FFTW_ENGINE_HH
#define FFTW_ENGINE_HH
/* -------------------------------------------------------------------------- */
#include "fft_engine.hh"
#include "fftw/interface.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
class FFTWEngine : public FFTEngine {
protected:
using plan_t = std::pair, fftw::plan>;
using complex_t = fftw::helper::complex;
/// Perform forward (R2C) transform
template
void forwardImpl(const Grid& real,
GridHermitian& spectral);
/// Perform backward (C2R) transform
template
void backwardImpl(Grid& real,
const GridHermitian& spectral);
/// Return the plans pair for a given transform signature
plan_t& getPlans(key_t key);
public:
/// Initialize with flags
explicit FFTWEngine(unsigned int flags = FFTW_ESTIMATE) noexcept
- : _flags(flags), plans() {}
+ : _flags(flags) {}
void forward(const Grid& real,
GridHermitian& spectral) override {
forwardImpl(real, spectral);
}
void forward(const Grid& real,
GridHermitian& spectral) override {
forwardImpl(real, spectral);
}
void backward(Grid& real,
GridHermitian& spectral) override {
backwardImpl(real, spectral);
}
void backward(Grid& real,
GridHermitian& spectral) override {
backwardImpl(real, spectral);
}
unsigned int flags() const { return _flags; }
/// Cast to FFTW complex type
static auto cast(Complex* data) { return reinterpret_cast(data); }
static auto cast(const Complex* data) {
return const_cast(reinterpret_cast(data));
}
protected:
unsigned int _flags; ///< FFTW flags
std::map plans; ///< plans corresponding to signatures
};
/* -------------------------------------------------------------------------- */
template
void FFTWEngine::forwardImpl(const Grid& real,
GridHermitian& spectral) {
auto& plans = getPlans(make_key(real, spectral));
fftw::execute(plans.first, const_cast(real.getInternalData()),
cast(spectral.getInternalData()));
}
template
void FFTWEngine::backwardImpl(Grid& real,
const GridHermitian& spectral) {
auto& plans = getPlans(make_key(real, spectral));
fftw::execute(plans.second, cast(spectral.getInternalData()),
real.getInternalData());
// Normalize
real *= (1. / real.getNbPoints());
}
} // namespace tamaas
#endif
diff --git a/src/core/grid_hermitian.hh b/src/core/grid_hermitian.hh
index 8d9375f..1a3a496 100644
--- a/src/core/grid_hermitian.hh
+++ b/src/core/grid_hermitian.hh
@@ -1,161 +1,159 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 GRID_HERMITIAN_HH
#define GRID_HERMITIAN_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "tamaas.hh"
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
static complex dummy;
/**
* @brief Multi-dimensional, multi-component herimitian array
*
* This class represents an array of hermitian data, meaning it has dimensions
* of: n1 * n2 * n3 * ... * (nx / 2 + 1)
*
* However, it acts as a fully dimensioned array, returning a dummy reference
* for data outside its real dimension, allowing one to write full (and
* inefficient) loops without really worrying about the reduced dimension.
*
* It would however be better to just use the true dimensions of the surface
* for all intents and purposes, as it saves computation time.
*/
template
class GridHermitian : public Grid, dim> {
public:
using value_type = complex;
/* ------------------------------------------------------------------------ */
/* Constructors */
/* ------------------------------------------------------------------------ */
public:
GridHermitian() = default;
GridHermitian(const GridHermitian& o) = default;
GridHermitian(GridHermitian&& o) noexcept = default;
using Grid::Grid;
using Grid, dim>::operator=;
/* ------------------------------------------------------------------------ */
/* Access operators */
/* ------------------------------------------------------------------------ */
template
inline complex& operator()(T1... args);
template
inline const complex& operator()(T1... args) const;
static std::array hermitianDimensions(std::array n) {
n[dim - 1] /= 2;
n[dim - 1] += 1;
return n;
}
template ::value>>
static std::vector hermitianDimensions(std::vector n) {
n.back() /= 2;
n.back() += 1;
return n;
}
private:
template
inline void packTuple(UInt* t, UInt i, T1... args) const;
template
inline void packTuple(UInt* t, UInt i) const;
};
/* -------------------------------------------------------------------------- */
/* Inline function definitions */
/* -------------------------------------------------------------------------- */
template
template
inline void GridHermitian::packTuple(UInt* t, UInt i,
T1... args) const {
*t = i;
packTuple(t + 1, args...);
}
template
template
inline void GridHermitian::packTuple(UInt* t, UInt i) const {
*t = i;
}
template
template
inline complex& GridHermitian::operator()(T1... args) {
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs <= dim + 1, "Too many arguments");
UInt last_index = std::get(std::forward_as_tuple(args...));
if (nargs != 1 && last_index >= this->n[dim - 1]) {
std::array tuple = {{0}};
packTuple(tuple.data(), args...);
for (UInt i = 0; i < dim; i++) {
if (tuple[i] && i != dim - 1)
tuple[i] = this->n[i] - tuple[i];
else if (tuple[i] && i == dim - 1)
tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
}
dummy = conj(this->Grid, dim>::operator()(tuple));
return dummy;
}
- else
return Grid, dim>::operator()(args...);
}
template
template
inline const complex& GridHermitian::operator()(T1... args) const {
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs <= dim + 1, "Too many arguments");
UInt last_index = std::get(std::forward_as_tuple(args...));
if (nargs != 1 && last_index >= this->n[dim - 1]) {
std::array tuple = {{0}};
packTuple(tuple.data(), args...);
for (UInt i = 0; i < dim; i++) {
if (tuple[i] && i != dim - 1)
tuple[i] = this->n[i] - tuple[i];
else if (tuple[i] && i == dim - 1)
tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
}
dummy = conj(this->Grid, dim>::operator()(tuple));
return dummy;
}
- else
return Grid, dim>::operator()(args...);
}
} // namespace tamaas
#endif // GRID_HERMITIAN_HH
diff --git a/src/core/loops/loop_utils.hh b/src/core/loops/loop_utils.hh
index 78bb12e..5fac93a 100644
--- a/src/core/loops/loop_utils.hh
+++ b/src/core/loops/loop_utils.hh
@@ -1,103 +1,103 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 LOOP_UTILS_HH
#define LOOP_UTILS_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "tamaas.hh"
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace detail {
template
UInt loopSize(Grids&&... grids) {
return (only_points)
? std::get<0>(std::forward_as_tuple(grids...)).getNbPoints()
: std::get<0>(std::forward_as_tuple(grids...)).dataSize();
}
template
-void areAllEqual(bool, UInt) {}
+void areAllEqual(bool /*unused*/ /*unused*/, UInt /*unused*/ /*unused*/) {}
template
void areAllEqual(bool result, UInt prev, UInt current) {
if (!(result && prev == current))
TAMAAS_EXCEPTION("Cannot loop over ranges that do not have the same size!");
}
template
void areAllEqual(bool result, UInt prev, UInt current, Sizes... rest) {
areAllEqual(result && prev == current, current, rest...);
}
} // namespace detail
template
void checkLoopSize(Ranges&&... ranges) {
detail::areAllEqual(true, (ranges.end() - ranges.begin())...);
}
namespace detail {
template
struct reduction_helper;
template
struct reduction_helper
: public thrust::plus {
ReturnType init() const {
return ReturnType(0);
}
};
template
struct reduction_helper
: public thrust::multiplies {
ReturnType init() const {
return ReturnType(1);
}
};
template
struct reduction_helper
: public thrust::minimum {
ReturnType init() const {
return std::numeric_limits::max();
}
};
template
struct reduction_helper
: public thrust::maximum {
ReturnType init() const {
return std::numeric_limits::lowest();
}
};
} // namespace detail
} // namespace tamaas
#endif // __LOOP_UTILS_HH
diff --git a/src/core/mpi_interface.hh b/src/core/mpi_interface.hh
index 46add03..41df35d 100644
--- a/src/core/mpi_interface.hh
+++ b/src/core/mpi_interface.hh
@@ -1,244 +1,246 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 MPI_INTERFACE_HH
#define MPI_INTERFACE_HH
/* -------------------------------------------------------------------------- */
#include "static_types.hh"
#include "tamaas.hh"
#include
#ifdef TAMAAS_USE_MPI
#include
#endif
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/// Contains mock mpi functions
namespace mpi_dummy {
struct comm {
static comm world;
};
struct sequential {
static void enter(){};
static void exit(){};
};
struct sequential_guard {
sequential_guard() { sequential::enter(); }
~sequential_guard() { sequential::exit(); }
};
enum class thread : int { single, funneled, serialized, multiple };
inline bool initialized() { return false; }
inline bool finalized() { return false; }
-inline int init(int*, char***) { return 0; }
-inline int init_thread(int*, char***, thread, thread* provided) {
+inline int init(int* /*unused*/, char*** /*unused*/) { return 0; }
+inline int init_thread(int* /*unused*/, char*** /*unused*/, thread /*unused*/,
+ thread* provided) {
*provided = thread::funneled;
return 0;
}
inline int finalize() { return 0; }
-inline int rank(comm = comm::world) { return 0; }
-inline int size(comm = comm::world) { return 1; }
+inline int rank(comm /*unused*/ = comm::world) { return 0; }
+inline int size(comm /*unused*/ = comm::world) { return 1; }
template
-inline decltype(auto) reduce(T&& val, comm = comm::world) {
+inline decltype(auto) reduce(T&& val, comm /*unused*/ = comm::world) {
return std::forward(val);
}
template
-inline decltype(auto) allreduce(T&& val, comm = comm::world) {
+inline decltype(auto) allreduce(T&& val, comm /*unused*/ = comm::world) {
return std::forward(val);
}
template
inline decltype(auto) gather(const T* send, T* recv, int count,
- comm = comm::world) {
+ comm /*unused*/ = comm::world) {
if (send == recv)
return;
thrust::copy_n(send, count, recv);
}
template
-inline decltype(auto) bcast(T*, int, comm = comm::world) {}
+inline decltype(auto) bcast(T* /*unused*/, int /*unused*/,
+ comm /*unused*/ = comm::world) {}
} // namespace mpi_dummy
/* -------------------------------------------------------------------------- */
#ifdef TAMAAS_USE_MPI
/// Contains real mpi functions
namespace mpi_impl {
/// MPI_Comm wrapper
struct comm {
MPI_Comm _comm;
operator MPI_Comm() { return _comm; }
static comm& world();
};
struct sequential {
static void enter() { comm::world()._comm = MPI_COMM_SELF; }
static void exit() { comm::world()._comm = MPI_COMM_WORLD; }
};
struct sequential_guard {
sequential_guard() { sequential::enter(); }
~sequential_guard() { sequential::exit(); }
};
/// MPI Thread level
enum class thread : int {
single = MPI_THREAD_SINGLE,
funneled = MPI_THREAD_FUNNELED,
serialized = MPI_THREAD_SERIALIZED,
multiple = MPI_THREAD_MULTIPLE
};
template
struct type_trait;
#define TYPE(t, mpi_t) \
template <> \
struct type_trait { \
static constexpr MPI_Datatype value = mpi_t; \
}
TYPE(double, MPI_DOUBLE);
TYPE(int, MPI_INT);
TYPE(unsigned int, MPI_UNSIGNED);
TYPE(long double, MPI_LONG_DOUBLE);
TYPE(long, MPI_LONG);
TYPE(unsigned long, MPI_UNSIGNED_LONG);
TYPE(::thrust::complex, MPI_CXX_DOUBLE_COMPLEX);
TYPE(::thrust::complex, MPI_CXX_LONG_DOUBLE_COMPLEX);
TYPE(bool, MPI_CXX_BOOL);
#undef TYPE
template
struct operation_trait;
#define OPERATION(op, mpi_op) \
template <> \
struct operation_trait { \
static constexpr MPI_Op value = mpi_op; \
}
OPERATION(plus, MPI_SUM);
OPERATION(min, MPI_MIN);
OPERATION(max, MPI_MAX);
OPERATION(times, MPI_PROD);
#undef OPERATION
inline bool initialized() {
int has_init;
MPI_Initialized(&has_init);
return has_init != 0;
}
inline bool finalized() {
int has_final;
MPI_Finalized(&has_final);
return has_final != 0;
}
inline int init(int* argc, char*** argv) { return MPI_Init(argc, argv); }
inline int init_thread(int* argc, char*** argv, thread required,
thread* provided) {
return MPI_Init_thread(argc, argv, static_cast(required),
reinterpret_cast(provided));
}
inline int finalize() { return MPI_Finalize(); }
inline int rank(comm communicator = comm::world()) {
int rank;
MPI_Comm_rank(communicator, &rank);
return rank;
}
inline int size(comm communicator = comm::world()) {
int size;
MPI_Comm_size(communicator, &size);
return size;
}
template
inline decltype(auto) reduce(T val, comm communicator = comm::world()) {
MPI_Reduce(&val, &val, 1, type_trait::value, operation_trait::value, 0,
communicator);
return val;
}
template ::value or
std::is_same::value>>
inline decltype(auto) allreduce(T val, comm communicator = comm::world()) {
MPI_Allreduce(&val, &val, 1, type_trait::value, operation_trait::value,
communicator);
return val;
}
template
inline decltype(auto) allreduce(const StaticVector& v,
comm communicator = comm::world()) {
Vector res;
MPI_Allreduce(v.begin(), res.begin(), n, type_trait::value,
operation_trait::value, communicator);
return res;
}
template
inline decltype(auto) allreduce(const StaticMatrix& v,
comm communicator = comm::world()) {
Matrix res;
MPI_Allreduce(v.begin(), res.begin(), n * m, type_trait::value,
operation_trait::value, communicator);
return res;
}
template
inline decltype(auto) gather(const T* send, T* recv, int count,
comm communicator = comm::world()) {
MPI_Gather(send, count, type_trait::value, recv, count,
type_trait::value, 0, communicator);
}
template
inline decltype(auto) bcast(T* buffer, int count,
comm communicator = comm::world()) {
MPI_Bcast(buffer, count, type_trait::value, 0, communicator);
}
} // namespace mpi_impl
namespace mpi = mpi_impl;
#else
namespace mpi = mpi_dummy;
#endif // TAMAAS_USE_MPI
} // namespace tamaas
/* -------------------------------------------------------------------------- */
#endif // MPI_INTERFACE_HH
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 80ca0bc..d84c453 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,213 +1,211 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 "statistics.hh"
#include "fft_engine.hh"
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
template
Real Statistics::computeRMSHeights(Grid& surface) {
return std::sqrt(surface.var());
}
template
Real Statistics::computeSpectralRMSSlope(Grid& surface) {
const auto h_size =
GridHermitian::hermitianDimensions(surface.sizes());
auto wavevectors =
FFTEngine::template computeFrequencies(h_size);
wavevectors *= 2 * M_PI; // need q for slopes
const auto psd = computePowerSpectrum(surface);
const Real rms_slope_mean = Loop::reduce(
[] CUDA_LAMBDA(VectorProxy q, const Complex& psd_val) {
// Checking if we're in the zone that does not have hermitian symmetry
if (std::abs(q.back()) < 1e-15)
return q.l2squared() * psd_val.real();
- else
- return 2 * q.l2squared() * psd_val.real();
+ return 2 * q.l2squared() * psd_val.real();
},
range>(wavevectors), psd);
return std::sqrt(rms_slope_mean);
}
/* -------------------------------------------------------------------------- */
template
GridHermitian
Statistics::computePowerSpectrum(Grid& surface) {
const auto h_size =
GridHermitian::hermitianDimensions(surface.sizes());
GridHermitian psd(h_size, surface.getNbComponents());
FFTEngine::makeEngine()->forward(surface, psd);
Real factor = 1. / surface.getGlobalNbPoints();
// Squaring the fourier transform of surface and normalizing
Loop::loop(
[factor] CUDA_LAMBDA(Complex & c) {
c *= factor;
c *= conj(c);
},
psd);
return psd;
}
/* -------------------------------------------------------------------------- */
template
Grid
Statistics::computeAutocorrelation(Grid& surface) {
Grid acf(surface.sizes(), surface.getNbComponents());
auto psd = computePowerSpectrum(surface);
FFTEngine::makeEngine()->backward(acf, psd);
acf *= acf.getGlobalNbPoints();
return acf;
}
/* -------------------------------------------------------------------------- */
template
Real Statistics::contact(const Grid& tractions,
UInt perimeter) {
Real points = 0;
UInt nc = tractions.getNbComponents();
switch (nc) {
case 1:
points = Loop::reduce(
[] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions);
break;
case 2:
points = Loop::reduce(
[] CUDA_LAMBDA(VectorProxy t) -> Real {
return t.back() > 0;
},
range>(tractions));
break;
case 3:
points = Loop::reduce(
[] CUDA_LAMBDA(VectorProxy t) -> Real {
return t.back() > 0;
},
range>(tractions));
break;
default:
TAMAAS_EXCEPTION("Invalid number of components in traction");
}
auto area = points / tractions.getGlobalNbPoints();
if (dim == 1)
perimeter = 0;
// Correction from Yastrebov et al. (Trib. Intl., 2017)
// 10.1016/j.triboint.2017.04.023
- return area -
- (M_PI - 1 + std::log(2)) / (24. * tractions.getGlobalNbPoints()) * perimeter;
+ return area - (M_PI - 1 + std::log(2)) /
+ (24. * tractions.getGlobalNbPoints()) * perimeter;
}
/* -------------------------------------------------------------------------- */
namespace {
template
class moment_helper {
public:
moment_helper(const std::array& exp) : exponent(exp) {}
CUDA_LAMBDA Complex operator()(VectorProxy q,
const Complex& phi) const {
Real mul = 1;
for (UInt i = 0; i < dim; ++i)
mul *= std::pow(q(i), exponent[i]);
// Do not duplicate everything from hermitian symmetry
if (std::abs(q.back()) < 1e-15)
return mul * phi;
- else
- return 2 * mul * phi;
+ return 2 * mul * phi;
}
private:
std::array exponent;
};
} // namespace
template <>
std::vector Statistics<1>::computeMoments(Grid& surface) {
constexpr UInt dim = 1;
std::vector moments(3);
const auto psd = computePowerSpectrum(surface);
auto wavevectors =
FFTEngine::template computeFrequencies(psd.sizes());
// we don't multiply by 2 pi because moments are computed with k
moments[0] = Loop::reduce(moment_helper{{{0}}},
range(wavevectors), psd)
.real();
moments[1] = Loop::reduce(moment_helper{{{2}}},
range(wavevectors), psd)
.real();
moments[2] = Loop::reduce(moment_helper{{{4}}},
range(wavevectors), psd)
.real();
return moments;
}
template <>
std::vector Statistics<2>::computeMoments(Grid& surface) {
constexpr UInt dim = 2;
std::vector moments(3);
const auto psd = computePowerSpectrum(surface);
auto wavevectors =
FFTEngine::template computeFrequencies(psd.sizes());
// we don't multiply by 2 pi because moments are computed with k
moments[0] = Loop::reduce(moment_helper{{{0, 0}}},
range(wavevectors), psd)
.real();
auto m02 = Loop::reduce(moment_helper{{{0, 2}}},
range(wavevectors), psd)
.real();
auto m20 = Loop::reduce(moment_helper{{{2, 0}}},
range(wavevectors), psd)
.real();
moments[1] = 0.5 * (m02 + m20);
auto m22 = Loop::reduce(moment_helper{{{2, 2}}},
range(wavevectors), psd)
.real();
auto m40 = Loop::reduce(moment_helper{{{4, 0}}},
range(wavevectors), psd)
.real();
auto m04 = Loop::reduce(moment_helper{{{0, 4}}},
range(wavevectors), psd)
.real();
moments[2] = (3 * m22 + m40 + m04) / 3.;
return moments;
}
template struct Statistics<1>;
template struct Statistics<2>;
} // namespace tamaas
diff --git a/src/model/adhesion_functional.hh b/src/model/adhesion_functional.hh
index fa4311f..4c45591 100644
--- a/src/model/adhesion_functional.hh
+++ b/src/model/adhesion_functional.hh
@@ -1,102 +1,99 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 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 ADHESION_FUNCTIONAL_HH
#define ADHESION_FUNCTIONAL_HH
/* -------------------------------------------------------------------------- */
#include "functional.hh"
#include