Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/surface/isopowerlaw.cpp
9a360c12e700features/netcdf-mpi

Last Change
Raw File

isopowerlaw.cpp
View Options

/*
* 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 <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include <map>
/* -------------------------------------------------------------------------- */
namespace tamaas {
template <UInt dim>
void Isopowerlaw<dim>::computeFilter(
GridHermitian<Real, dim>& filter_coefficients) const {
Filter<dim>::computeFilter(
[this] CUDA_LAMBDA(Complex & coeff, VectorProxy<Real, dim> q) {
coeff = (*this)(q);
},
filter_coefficients);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::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<Real> Isopowerlaw<2>::moments() const {
std::map<UInt, Real> T;
T[0] = 2 * M_PI;
T[2] = M_PI;
T[4] = 3 * M_PI / 4.;
Real q1 = static_cast<Real>(this->q1); // shadowing this->q1
Real xi = q0 / q1;
Real zeta = q2 / q1;
Real A = std::pow(q1, 2 * (hurst + 1));
std::vector<Real> 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<Real> Isopowerlaw<1>::moments() const {
TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces");
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::alpha() const {
const auto m = moments();
return m[0] * m[2] / (m[1] * m[1]);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::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
c4science · Help