Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/isopowerlaw.cpp
f05be88c05d7archive/legacy_mpi

Last Change
Raw File

isopowerlaw.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include "fftransform.hh"
#include <map>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <UInt dim>
void Isopowerlaw<dim>::computeFilter(GridHermitian<Real, dim> & filter_coefficients) const {
const Real l[dim] = {1.};
Grid<Real, dim> wavevectors(filter_coefficients.sizes(), l, dim);
FFTransform<Real, dim>::template computeFrequencies<true>(wavevectors);
auto end = wavevectors.end(dim);
#pragma omp parallel for
for (auto it = wavevectors.begin(dim) ; it < end ; ++it) {
VectorProxy<Real> q(*it, dim);
UInt i = wavevectors.getNbPoints() - (end - it);
filter_coefficients(i) = (*this)(q);
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::rmsHeights() const {
return std::sqrt(M_PI * ( (hurst + 1)/hurst * q1*q1
- 1./hurst * std::pow(q1, 2*(hurst+1)) * std::pow(q2, -2*hurst)
- q0*q0 ));
}
/* -------------------------------------------------------------------------- */
/**
* 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 xi = q0 / q1;
Real zeta = q2 / q1;
std::vector<Real> moments;
using namespace std;
for (int q : {0, 2, 4}) {
Real m = T[q] * pow(q1, q-2*hurst) * ((1-pow(xi, q+2))/(q+2) +
(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");
return std::vector<Real>();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::alpha() const {
std::vector<Real> m = moments();
return m[0] * m[2] / (m[1]*m[1]);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Isopowerlaw<dim>::rmsSlopes() const {
std::vector<Real> m = moments();
return std::sqrt(M_PI * m[1] / 2);
}
template class Isopowerlaw<1>;
template class Isopowerlaw<2>;
__END_TAMAAS__
c4science · Help