Page MenuHomec4science
Files F66182994

isopowerlaw.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Jun 8, 20:44

isopowerlaw.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* 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]);
}
/* -------------------------------------------------------------------------- */
template <>
Real Isopowerlaw<2>::radialPSDMoment(Real q) const {
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));
Real exp = q - 2 * hurst;
Real zeta_pow = [](Real zeta, Real exp) {
if (exp == 0.) // handle edge case when H = 0.5 and q = 1
return std::log(zeta);
return (std::pow(zeta, exp) - 1) / exp;
}(zeta, exp);
return A * std::pow(q1, exp) *
((1 - std::pow(xi, q + 2)) / (q + 2) + zeta_pow);
}
template <>
Real Isopowerlaw<1>::radialPSDMoment(Real /*q*/) const {
throw not_implemented_error{
TAMAAS_MSG("Moments have not been implemented for 1D surfaces")};
}
/* -------------------------------------------------------------------------- */
/// \cond DO_NOT_DOCUMENT
/*
* 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;
std::vector<Real> moments;
moments.reserve(3);
for (UInt q : {0, 2, 4}) {
moments.push_back(radialPSDMoment(q) * T[q]);
}
return moments;
}
template <>
std::vector<Real> Isopowerlaw<1>::moments() const {
throw not_implemented_error{
TAMAAS_MSG("Moments have not been implemented for 1D surfaces")};
}
/// \endcond
/* -------------------------------------------------------------------------- */
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 <UInt dim>
Real Isopowerlaw<dim>::elasticEnergy() const {
auto integral_2d = 2 * M_PI * radialPSDMoment(1);
auto q_converted = 2 * M_PI * integral_2d;
return q_converted / 4;
}
template class Isopowerlaw<1>;
template class Isopowerlaw<2>;
} // namespace tamaas

Event Timeline

c4science · Help