diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index 996b41ba..525607d0 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,88 +1,89 @@
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 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 "surface_statistics.hh"
#include "statistics.hh"
#include "wrap.hh"
#include <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace wrap {
/* -------------------------------------------------------------------------- */
template <UInt dim>
void wrapStatistics(py::module& mod) {
auto name = makeDimensionName("Statistics", dim);
py::class_<Statistics<dim>>(mod, name.c_str())
.def_static("computePowerSpectrum",
&Statistics<dim>::computePowerSpectrum)
.def_static("computeAutocorrelation",
&Statistics<dim>::computeAutocorrelation)
.def_static("computeMoments", &Statistics<dim>::computeMoments)
.def_static("computeSpectralRMSSlope",
- &Statistics<dim>::computeSpectralRMSSlope);
+ &Statistics<dim>::computeSpectralRMSSlope)
+ .def_static("computeRMSHeights", &Statistics<dim>::computeRMSHeights);
}
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod) {
// Exposing SurfaceStatistics (legacy)
py::class_<SurfaceStatistics>(mod, "SurfaceStatistics")
.def_static("computeMaximum", &SurfaceStatistics::computeMaximum)
.def_static("computeMinimum", &SurfaceStatistics::computeMinimum)
.def_static("computeSpectralRMSSlope",
&SurfaceStatistics::computeSpectralRMSSlope)
.def_static("computeRMSSlope", &SurfaceStatistics::computeRMSSlope)
.def_static("computeMoments", &SurfaceStatistics::computeMoments)
.def_static("computeSkewness", &SurfaceStatistics::computeSkewness)
.def_static("computeKurtosis", &SurfaceStatistics::computeKurtosis)
.def_static("computeSpectralMeanCurvature",
&SurfaceStatistics::computeSpectralMeanCurvature)
.def_static("computeSpectralStdev",
&SurfaceStatistics::computeSpectralStdev)
.def_static("computeAnalyticFractalMoment",
&SurfaceStatistics::computeAnalyticFractalMoment)
.def_static("computePerimeter", &SurfaceStatistics::computePerimeter)
.def_static("computeContactArea", &SurfaceStatistics::computeContactArea)
.def_static("computeContactAreaRatio",
&SurfaceStatistics::computeContactAreaRatio)
.def_static("computeSpectralDistribution",
&SurfaceStatistics::computeSpectralDistribution)
.def_static("computeSum", &SurfaceStatistics::computeSum)
.def_static("computeAutocorrelation",
&SurfaceStatistics::computeAutocorrelation)
.def_static("computePowerSpectrum",
&SurfaceStatistics::computePowerSpectrum);
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 7b622629..5975865a 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,170 +1,184 @@
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 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 "wrap.hh"
#include "isopowerlaw.hh"
#include "surface_generator_filter.hh"
#include "surface_generator_filter_fft.hh"
+#include "surface_generator_random_phase.hh"
+#include "wrap.hh"
#include <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template <UInt dim>
class PyFilter : public Filter<dim> {
public:
using Filter<dim>::Filter;
// Overriding pure virtual functions
void
computeFilter(GridHermitian<Real, dim>& filter_coefficients) const override {
PYBIND11_OVERLOAD_PURE(void, Filter<dim>, computeFilter,
filter_coefficients);
}
};
template <UInt dim>
void wrapFilter(py::module& mod) {
auto name = makeDimensionName("Filter", dim);
py::class_<Filter<dim>, PyFilter<dim>>(mod, name.c_str())
- .def(py::init<>())
- .def("computeFilter", &Filter<dim>::computeFilter);
+ .def(py::init<>())
+ .def("computeFilter", &Filter<dim>::computeFilter);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void wrapIsopowerlaw(py::module& mod) {
std::string name = makeDimensionName("Isopowerlaw", dim);
py::class_<Isopowerlaw<dim>, Filter<dim>>(mod, name.c_str())
.def(py::init<>())
.def_property("q0", &Isopowerlaw<dim>::getQ0, &Isopowerlaw<dim>::setQ0)
.def_property("q1", &Isopowerlaw<dim>::getQ1, &Isopowerlaw<dim>::setQ1)
.def_property("q2", &Isopowerlaw<dim>::getQ2, &Isopowerlaw<dim>::setQ2)
.def_property("hurst", &Isopowerlaw<dim>::getHurst,
&Isopowerlaw<dim>::setHurst)
.def("rmsHeights", &Isopowerlaw<dim>::rmsHeights)
.def("moments", &Isopowerlaw<dim>::moments)
.def("alpha", &Isopowerlaw<dim>::alpha)
.def("rmsSlopes", &Isopowerlaw<dim>::rmsSlopes)
// legacy wrapper code
.def("getQ0", &Isopowerlaw<dim>::getQ0,
py::return_value_policy::reference)
.def("getQ1", &Isopowerlaw<dim>::getQ1,
py::return_value_policy::reference)
.def("getQ2", &Isopowerlaw<dim>::getQ2,
py::return_value_policy::reference)
.def("setQ0", &Isopowerlaw<dim>::setQ0,
py::return_value_policy::reference)
.def("setQ1", &Isopowerlaw<dim>::setQ1,
py::return_value_policy::reference)
.def("setQ2", &Isopowerlaw<dim>::setQ2,
py::return_value_policy::reference)
.def("setHurst", &Isopowerlaw<dim>::setHurst)
.def("getHurst", &Isopowerlaw<dim>::getHurst,
py::return_value_policy::reference);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
-void wrapSurfaceGeneratorFilter(py::module& mod) {
- std::string name = makeDimensionName("SurfaceGeneratorFilter", dim);
- py::class_<SurfaceGeneratorFilter<dim>>(mod, name.c_str())
- .def(py::init<>())
- .def("buildSurface", &SurfaceGeneratorFilter<dim>::buildSurface,
- "Build rough surface")
- .def("setFilter", &SurfaceGeneratorFilter<dim>::setFilter,
- "Set PSD filter", "filter"_a)
+void wrapSurfaceGenerators(py::module& mod) {
+ std::string generator_name = makeDimensionName("SurfaceGenerator", dim);
+ py::class_<SurfaceGenerator<dim>>(mod, generator_name.c_str())
+ .def("buildSurface", &SurfaceGenerator<dim>::buildSurface)
.def("setSizes",
py::overload_cast<const std::array<UInt, dim>&>(
- &SurfaceGeneratorFilter<dim>::setSizes),
+ &SurfaceGenerator<dim>::setSizes),
"surface_sizes"_a)
- .def_property("random_seed", &SurfaceGeneratorFilter<dim>::getRandomSeed,
- &SurfaceGeneratorFilter<dim>::setRandomSeed)
+ .def_property("random_seed", &SurfaceGenerator<dim>::getRandomSeed,
+ &SurfaceGenerator<dim>::setRandomSeed);
+
+ std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim);
+ py::class_<SurfaceGeneratorFilter<dim>, SurfaceGenerator<dim>>(
+ mod, filter_name.c_str())
+ .def(py::init<>())
+ .def("setFilter", &SurfaceGeneratorFilter<dim>::setFilter,
+ "Set PSD filter", "filter"_a)
+ .def("setSpectrum", &SurfaceGeneratorFilter<dim>::setFilter,
+ "Set PSD spectrum", "filter"_a)
// legacy wrapper code
.def("setRandomSeed", &SurfaceGeneratorFilter<dim>::setRandomSeed,
"[[Deprecated: use property]]");
+
+ std::string random_name =
+ makeDimensionName("SurfaceGeneratorRandomPhase", dim);
+ py::class_<SurfaceGeneratorRandomPhase<dim>, SurfaceGenerator<dim>>(
+ mod, random_name.c_str())
+ .def(py::init<>())
+ .def("setSpectrum", &SurfaceGeneratorRandomPhase<dim>::setSpectrum,
+ "Set PSD Spectrum", "spectrum"_a);
}
/* -------------------------------------------------------------------------- */
// Legacy wrap
void wrapSurfaceGeneratorFilterFFT(py::module& mod) {
py::class_<SurfaceGeneratorFilterFFT>(mod, "SurfaceGeneratorFilterFFT")
.def(py::init<>())
.def("getQ0",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<UInt>(&f.getQ0());
})
.def("getQ1",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<UInt>(&f.getQ1());
})
.def("getQ2",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<UInt>(&f.getQ2());
})
.def("getHurst",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<Real>(&f.getHurst());
})
.def("analyticalRMS", &SurfaceGeneratorFilterFFT::analyticalRMS)
.def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface)
.def("getGridSize",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<UInt>(&f.getGridSize());
})
.def("getRandomSeed",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<long>(&f.getRandomSeed());
})
.def("getRMS",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer<Real>(&f.getRMS());
})
.def("Init", &SurfaceGeneratorFilterFFT::Init);
}
/* -------------------------------------------------------------------------- */
void wrapSurface(py::module& mod) {
wrapFilter<1>(mod);
wrapFilter<2>(mod);
wrapIsopowerlaw<1>(mod);
wrapIsopowerlaw<2>(mod);
- wrapSurfaceGeneratorFilter<1>(mod);
- wrapSurfaceGeneratorFilter<2>(mod);
+ wrapSurfaceGenerators<1>(mod);
+ wrapSurfaceGenerators<2>(mod);
// legacy wrap
wrapSurfaceGeneratorFilterFFT(mod);
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/SConscript b/src/SConscript
index b005dcfb..7e220943 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,131 +1,132 @@
import os
Import('main_env')
def prepend(path, list):
return [os.path.join(path, x) for x in list]
env = main_env.Clone()
# Core
core_list = """
fft_plan_manager.cpp
fftransform.cpp
fftransform_fftw.cpp
grid.cpp
grid_hermitian.cpp
statistics.cpp
surface.cpp
tamaas.cpp
legacy_types.cpp
loop.cpp
""".split()
core_list = prepend('core', core_list)
core_list.append('tamaas_info.cpp')
# Lib roughcontact
generator_list = """
surface_generator.cpp
surface_generator_filter.cpp
surface_generator_filter_fft.cpp
+surface_generator_random_phase.cpp
isopowerlaw.cpp
""".split()
generator_list = prepend('surface', generator_list)
# Lib PERCOLATION
percolation_list = """
flood_fill.cpp
""".split()
percolation_list = prepend('percolation', percolation_list)
# BEM PERCOLATION
bem_list = """
bem_kato.cpp
bem_polonski.cpp
bem_gigi.cpp
bem_gigipol.cpp
bem_penalty.cpp
bem_uzawa.cpp
bem_fft_base.cpp
bem_functional.cpp
bem_meta_functional.cpp
elastic_energy_functional.cpp
exponential_adhesion_functional.cpp
squared_exponential_adhesion_functional.cpp
maugis_adhesion_functional.cpp
complimentary_term_functional.cpp
bem_grid.cpp
bem_grid_polonski.cpp
bem_grid_kato.cpp
bem_grid_teboulle.cpp
bem_grid_condat.cpp
""".split()
bem_list = prepend('bem', bem_list)
# Model
model_list = """
model.cpp
model_factory.cpp
model_type.cpp
model_template.cpp
be_engine.cpp
westergaard.cpp
elastic_functional.cpp
meta_functional.cpp
adhesion_functional.cpp
volume_potential.cpp
kelvin.cpp
mindlin.cpp
boussinesq.cpp
elasto_plastic/isotropic_hardening.cpp
elasto_plastic/elasto_plastic_model.cpp
elasto_plastic/residual.cpp
""".split()
model_list = prepend('model', model_list)
# Solvers
solvers_list = """
contact_solver.cpp
polonsky_keer_rey.cpp
kato.cpp
beck_teboulle.cpp
condat.cpp
polonsky_keer_tan.cpp
ep_solver.cpp
epic.cpp
""".split()
solvers_list = prepend('solvers', solvers_list)
# GPU API
gpu_list = """
fftransform_cufft.cpp
""".split()
gpu_list = prepend('gpu', gpu_list)
# Assembling total list
rough_contact_list = \
core_list + generator_list + percolation_list + model_list + solvers_list
# For some reason collapse OMP loops don't compile on intel
if env['CXX'] != 'icpc':
rough_contact_list += bem_list
# Adding GPU if needed
if env['backend'] == 'cuda':
rough_contact_list += gpu_list
# Generating dependency files
# env.AppendUnique(CXXFLAGS=['-MMD'])
# for file_name in rough_contact_list:
# obj = file_name.replace('.cpp', '.os')
# dep_file_name = file_name.replace('.cpp', '.d')
# env.SideEffect(dep_file_name, obj)
# env.ParseDepends(dep_file_name)
# Adding extra warnings for Tamaas base lib
env.AppendUnique(CXXFLAGS=['-Wextra'])
libTamaas = env.SharedLibrary('Tamaas', rough_contact_list)
Export('libTamaas')
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 79a5583d..61d5df62 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,159 +1,162 @@
/**
* @file
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-2017 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 "statistics.hh"
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
+template <UInt dim>
+Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
+ return std::sqrt(surface.var());
+}
+
template <UInt dim>
Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
auto h_size = GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
auto wavevectors =
FFTransform<Real, dim>::template computeFrequencies<true>(h_size);
auto psd = computePowerSpectrum(surface);
- Real rms_slope_mean = 0;
-
- rms_slope_mean = Loop::reduce<operation::plus>(
+ Real rms_slope_mean = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, dim> q, Complex & psd_val) {
return q.l2squared() * psd_val.real();
},
range<VectorProxy<Real, dim>>(wavevectors), psd);
- rms_slope_mean /= wavevectors.getNbPoints(); // average
+ //rms_slope_mean /= wavevectors.getNbPoints(); // average
rms_slope_mean *= 2; // because we have a hermitian surface
- return rms_slope_mean;
+ return std::sqrt(rms_slope_mean);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
GridHermitian<Real, dim>
Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
auto h_size = GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
GridHermitian<Real, dim> psd(h_size, surface.getNbComponents());
FFTPlanManager::get().createPlan(surface, psd).forwardTransform();
Real factor = surface.getNbPoints();
factor = 1. / (factor * factor);
// Squaring the fourier transform of surface and normalizing
Loop::loop(
[factor] CUDA_LAMBDA(Complex & c) {
c *= conj(c);
c *= factor;
},
psd);
FFTPlanManager::get().destroyPlan(surface, psd);
return psd;
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Grid<Real, dim>
Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
auto psd = computePowerSpectrum(surface);
FFTPlanManager::get().createPlan(acf, psd).backwardTransform();
acf *= acf.getNbPoints();
FFTPlanManager::get().destroyPlan(acf, psd);
return acf;
}
/* -------------------------------------------------------------------------- */
namespace {
template <UInt dim>
class moment_helper {
public:
moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
CUDA_LAMBDA Complex operator()(VectorProxy<const Real, dim> q,
const Complex& phi) const {
Real mul = 1;
for (UInt i = 0; i < dim; ++i)
mul *= std::pow(q(i), exponent[i]);
return mul * phi;
}
private:
std::array<UInt, dim> exponent;
};
} // namespace
template <>
std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
constexpr UInt dim = 1;
std::vector<Real> moments(3);
const auto psd = computePowerSpectrum(surface);
const auto wavevectors =
FFTransform<Real, dim>::template computeFrequencies<true>(psd.sizes());
moments[0] = psd.sum().real();
moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
range<PVector>(wavevectors), psd)
.real();
moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
range<PVector>(wavevectors), psd)
.real();
return moments;
}
template <>
std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
constexpr UInt dim = 2;
std::vector<Real> moments(3);
const auto psd = computePowerSpectrum(surface);
const auto wavevectors =
FFTransform<Real, dim>::template computeFrequencies<true>(psd.sizes());
moments[0] = psd.sum().real();
auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
range<PVector>(wavevectors), psd)
.real();
auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
range<PVector>(wavevectors), psd)
.real();
moments[1] = 0.5 * (m02 + m20);
auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
range<PVector>(wavevectors), psd)
.real();
auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
range<PVector>(wavevectors), psd)
.real();
auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
range<PVector>(wavevectors), psd)
.real();
moments[2] = (3 * m22 + m40 + m04) / 3.;
return moments;
}
template struct Statistics<1>;
template struct Statistics<2>;
__END_TAMAAS__
diff --git a/src/core/statistics.hh b/src/core/statistics.hh
index 3b98ad3b..7fcd6a67 100644
--- a/src/core/statistics.hh
+++ b/src/core/statistics.hh
@@ -1,57 +1,59 @@
/**
* @file
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-2017 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __STATISTICS_HH__
#define __STATISTICS_HH__
/* -------------------------------------------------------------------------- */
#include "fft_plan_manager.hh"
#include "grid.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
/**
* @brief Suitcase class for all statistics related functions
*/
template <UInt dim>
struct Statistics {
+ /// Compute hrms
+ static Real computeRMSHeights(Grid<Real, dim>& surface);
/// Compute hrms' in fourier space
static Real computeSpectralRMSSlope(Grid<Real, dim>& surface);
/// Compute PSD of surface
static GridHermitian<Real, dim>
computePowerSpectrum(Grid<Real, dim>& surface);
/// Compute autocorrelation
static Grid<Real, dim> computeAutocorrelation(Grid<Real, dim>& surface);
/// Compute spectral moments
static std::vector<Real> computeMoments(Grid<Real, dim>& surface);
using PVector = VectorProxy<const Real, dim>;
};
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif // __STATISTICS_HH__
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index 9044ab1e..31634f51 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,158 +1,159 @@
/**
* @mainpage Welcome to Tamaas !
*
* @section Introduction
* Tamaas is a spectral-boundary-element based contact library. It is made with
* love to be fast and friendly!
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Lucas Frérot <lucas.frerot@epfl.ch>
* @author Valentine Rey <valentine.rey@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-2017 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __TAMAAS_HH__
#define __TAMAAS_HH__
/* -------------------------------------------------------------------------- */
// Standard includes
#include <exception>
#include <iostream>
#include <memory>
#include <string>
/* -------------------------------------------------------------------------- */
// Special thrust includes
#include <thrust/complex.h>
#include <thrust/random.h>
#ifdef USE_CUDA
#include "unified_allocator.hh"
#else
#include "fftw_allocator.hh"
#endif
/// @section Namespace macros
#define __BEGIN_TAMAAS__ namespace tamaas {
#define __END_TAMAAS__ }
/// @section Cuda specific definitions
#define CUDA_LAMBDA __device__ __host__
#ifdef USE_CUDA
#define DEFAULT_ALLOCATOR UnifiedAllocator<T>
#else
#define DEFAULT_ALLOCATOR FFTWAllocator<T>
#endif
/* -------------------------------------------------------------------------- */
/// @section Standard definitions pulled from C++14
namespace std {
#if __cplusplus < 201402L
/// exchange
template <class T, class U = T>
T exchange(T& obj, U&& new_value) {
T old_value = move(obj);
obj = forward<U>(new_value);
return old_value;
}
/// make_unique
template <typename T, typename... Args>
unique_ptr<T> make_unique(Args&&... args) {
return unique_ptr<T>(new T(forward<Args>(args)...));
}
#endif
}
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
/// @section Common types definitions
using Real = double; ///< default floating point type
using UInt = unsigned int; ///< default unsigned integer type
using Int = int; ///< default signed integer type
template <typename T>
using complex = thrust::complex<T>; ///< template complex wrapper
using Complex = complex<Real>; ///< default floating point complex type
static constexpr Real zero_threshold = 1e-14;
/// @section Defining random toolbox
using ::thrust::random::normal_distribution;
+using ::thrust::random::uniform_real_distribution;
using random_engine = ::thrust::random::default_random_engine;
/* -------------------------------------------------------------------------- */
/// initialize tamaas (0 threads => let OMP_NUM_THREADS decide)
void initialize(UInt num_threads = 0);
/// cleanup tamaas
void finalize();
/* -------------------------------------------------------------------------- */
/// Generic exception class
class Exception : public std::exception {
public:
/// Constructor
Exception(const std::string& mesg) : msg(mesg) {}
virtual const char* what() const noexcept { return msg.c_str(); }
virtual ~Exception() = default;
private:
std::string msg; ///< message of exception
};
/* -------------------------------------------------------------------------- */
/// Enumeration of reduction operations
enum class operation { plus, times, min, max };
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
/// @section Convenience macros
#define TAMAAS_EXCEPTION(mesg) \
{ \
std::stringstream sstr; \
sstr << __FILE__ << ":" << __LINE__ << ":FATAL: " << mesg << std::endl; \
std::cerr.flush(); \
throw ::tamaas::Exception(sstr.str()); \
}
#define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
#if defined(TAMAAS_DEBUG)
#define TAMAAS_ASSERT(cond, reason) \
do { \
if (!(cond)) { \
TAMAAS_EXCEPTION(#cond " assert failed: " << reason); \
} \
} while (0)
#define TAMAAS_DEBUG_EXCEPTION(reason) TAMAAS_EXCEPTION(reason)
#else
#define TAMAAS_ASSERT(cond, reason)
#define TAMAAS_DEBUG_EXCEPTION(reason)
#endif
#define TAMAAS_ACCESSOR(var, type, name) \
type& get##name() { return var; } \
void set##name(const type& new_var) { var = new_var; }
/* -------------------------------------------------------------------------- */
#endif // TAMAAS_HH
diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp
index ee43a104..55e2358b 100644
--- a/src/surface/isopowerlaw.cpp
+++ b/src/surface/isopowerlaw.cpp
@@ -1,110 +1,110 @@
/**
* @file
*
* @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 {
auto wavevectors = FFTransform<Real, dim>::template computeFrequencies<true>(
filter_coefficients.sizes());
Loop::loop( // loop here cannot capture this
[this] CUDA_LAMBDA(Complex & coeff, VectorProxy<Real, dim> q) {
coeff = (*this)(q);
},
filter_coefficients, range<VectorProxy<Real, dim>>(wavevectors));
}
/* -------------------------------------------------------------------------- */
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;
-
+ moments.reserve(3);
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");
}
/* -------------------------------------------------------------------------- */
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(2 * m[1]);
}
template class Isopowerlaw<1>;
template class Isopowerlaw<2>;
__END_TAMAAS__
diff --git a/src/surface/surface_generator_filter.hh b/src/surface/surface_generator_filter.hh
index 884202d6..304be513 100644
--- a/src/surface/surface_generator_filter.hh
+++ b/src/surface/surface_generator_filter.hh
@@ -1,81 +1,81 @@
/**
* @file
*
* @author Guillaume Anciaux <guillaume.anciaux@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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef SURFACE_GENERATOR_FILTER_H
#define SURFACE_GENERATOR_FILTER_H
/* -------------------------------------------------------------------------- */
#include "filter.hh"
#include "surface_generator.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <UInt dim>
class SurfaceGeneratorFilter : public SurfaceGenerator<dim> {
public:
/// Default constructor
SurfaceGeneratorFilter() = default;
/// Destructor
virtual ~SurfaceGeneratorFilter();
public:
/// Build surface with Hu & Tonder algorithm
Grid<Real, dim>& buildSurface() override;
- /// Set filter object (need namespace for swig)
- void setFilter(::tamaas::Filter<dim>* new_filter);
+ /// Set filter object
+ void setFilter(Filter<dim>* new_filter);
protected:
/// Apply filter coefficients on white noise
void applyFilterOnSource();
/// Generate white noise with given distribution
template <typename T>
void generateWhiteNoise();
protected:
Filter<dim>* filter = nullptr; ///< not owned
GridHermitian<Real, dim> filter_coefficients;
Grid<Real, dim> white_noise;
};
/* -------------------------------------------------------------------------- */
/* Template implementations */
/* -------------------------------------------------------------------------- */
template <UInt dim>
template <typename T>
void SurfaceGeneratorFilter<dim>::generateWhiteNoise() {
random_engine gen(this->random_seed);
T distribution;
for (auto& noise : white_noise)
noise = distribution(gen);
}
__END_TAMAAS__
#endif
diff --git a/src/surface/surface_generator_random_phase.cpp b/src/surface/surface_generator_random_phase.cpp
new file mode 100644
index 00000000..1dcdd0f9
--- /dev/null
+++ b/src/surface/surface_generator_random_phase.cpp
@@ -0,0 +1,76 @@
+/**
+ * @file
+ *
+ * @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 "surface_generator_random_phase.hh"
+#include "fft_plan_manager.hh"
+#include "fftransform.hh"
+#include <iostream>
+/* -------------------------------------------------------------------------- */
+namespace tamaas {
+/* -------------------------------------------------------------------------- */
+
+template <UInt dim>
+Grid<Real, dim>& SurfaceGeneratorRandomPhase<dim>::buildSurface() {
+ if (this->grid.dataSize() == 0)
+ TAMAAS_EXCEPTION("the size of the grid is zero, did you call setSizes() ?");
+ if (this->filter == nullptr)
+ TAMAAS_EXCEPTION("filter is null, did you call setFilter() ?");
+
+ // Resizing arrays
+ auto hermitian_dim =
+ GridHermitian<Real, dim>::hermitianDimensions(this->grid.sizes());
+ this->white_noise.resize(hermitian_dim);
+ this->surface_coefficients.resize(hermitian_dim);
+
+ /// Generating white noise with uniform distribution for phase
+ this->generateWhiteNoise<uniform_real_distribution<Real>>();
+
+ /// Applying filters
+ filter->computeFilter(this->surface_coefficients);
+
+ Loop::loop([](Complex& coeff,
+ Real& phase) { coeff *= thrust::polar(1., 2 * M_PI * phase); },
+ this->surface_coefficients, this->white_noise);
+
+ FFTPlanManager::get()
+ .createPlan(this->grid, this->surface_coefficients)
+ .backwardTransform();
+ FFTPlanManager::get().destroyPlan(this->grid, this->surface_coefficients);
+
+ /// Normalizing surface for hrms independent of number of points
+ this->grid *= this->grid.dataSize();
+
+ return this->grid;
+}
+
+/* -------------------------------------------------------------------------- */
+template class SurfaceGeneratorRandomPhase<1>;
+template class SurfaceGeneratorRandomPhase<2>;
+/* -------------------------------------------------------------------------- */
+
+} // namespace tamaas
diff --git a/src/surface/surface_generator_filter.hh b/src/surface/surface_generator_random_phase.hh
similarity index 73%
copy from src/surface/surface_generator_filter.hh
copy to src/surface/surface_generator_random_phase.hh
index 884202d6..f4e86376 100644
--- a/src/surface/surface_generator_filter.hh
+++ b/src/surface/surface_generator_random_phase.hh
@@ -1,81 +1,75 @@
/**
* @file
*
- * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
+ * @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/>.
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef SURFACE_GENERATOR_FILTER_H
-#define SURFACE_GENERATOR_FILTER_H
+#ifndef SURFACE_GENERATOR_RANDOM_PHASE_H
+#define SURFACE_GENERATOR_RANDOM_PHASE_H
/* -------------------------------------------------------------------------- */
#include "filter.hh"
#include "surface_generator.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
+namespace tamaas {
+/* -------------------------------------------------------------------------- */
template <UInt dim>
-class SurfaceGeneratorFilter : public SurfaceGenerator<dim> {
-public:
- /// Default constructor
- SurfaceGeneratorFilter() = default;
- /// Destructor
- virtual ~SurfaceGeneratorFilter();
-
+class SurfaceGeneratorRandomPhase : public SurfaceGenerator<dim> {
public:
- /// Build surface with Hu & Tonder algorithm
+ /// Build surface with uniform random phase
Grid<Real, dim>& buildSurface() override;
- /// Set filter object (need namespace for swig)
- void setFilter(::tamaas::Filter<dim>* new_filter);
+ /// Set filter object
+ void setSpectrum(Filter<dim>* new_filter) {
+ filter = new_filter;
+ }
protected:
- /// Apply filter coefficients on white noise
- void applyFilterOnSource();
/// Generate white noise with given distribution
template <typename T>
void generateWhiteNoise();
protected:
Filter<dim>* filter = nullptr; ///< not owned
- GridHermitian<Real, dim> filter_coefficients;
+ GridHermitian<Real, dim> surface_coefficients;
Grid<Real, dim> white_noise;
};
/* -------------------------------------------------------------------------- */
/* Template implementations */
/* -------------------------------------------------------------------------- */
template <UInt dim>
template <typename T>
-void SurfaceGeneratorFilter<dim>::generateWhiteNoise() {
+void SurfaceGeneratorRandomPhase<dim>::generateWhiteNoise() {
random_engine gen(this->random_seed);
T distribution;
for (auto& noise : white_noise)
noise = distribution(gen);
}
-__END_TAMAAS__
+} // namespace tamaas
#endif