diff --git a/examples/rough_contact.py b/examples/rough_contact.py
index 1374ec4d..10f145ba 100644
--- a/examples/rough_contact.py
+++ b/examples/rough_contact.py
@@ -1,74 +1,74 @@
#!/usr/bin/env python3
# @file
# @section LICENSE
#
# Copyright (©) 2016-19 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/>.
import tamaas as tm
-import sys
import matplotlib.pyplot as plt
# Initialize threads and fftw
tm.initialize()
+tm.set_log_level(tm.LogLevel.info) # Show progression of solver
# Surface size
n = 512
# Surface generator
sg = tm.SurfaceGeneratorFilter2D()
sg.setSizes([n, n])
sg.random_seed = 0
# Spectrum
spectrum = tm.Isopowerlaw2D()
sg.setFilter(spectrum)
# Parameters
spectrum.q0 = 16
spectrum.q1 = 16
spectrum.q2 = 64
spectrum.hurst = 0.8
# Generating surface
surface = sg.buildSurface()
#surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
surface /= n
#print(spectrum.rmsSlopes())
#print(tm.SurfaceStatistics.computeRMSSlope(surface))
plt.imshow(surface)
# Creating model
model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n])
# Solver
solver = tm.PolonskyKeerRey(model, surface, 1e-12,
tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.pressure)
# Solve for target pressure
p_target = 0.1
solver.solve(p_target)
plt.figure()
plt.imshow(model.getTraction())
print(model.getTraction().mean())
# Cleanup threads
tm.finalize()
plt.show()
diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index 57b0930d..4f344c64 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,103 +1,112 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "logger.hh"
#include "statistics.hh"
#include "surface_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,
py::return_value_policy::move)
.def_static("computeAutocorrelation",
&Statistics<dim>::computeAutocorrelation,
py::return_value_policy::move)
.def_static("computeMoments", &Statistics<dim>::computeMoments)
.def_static("computeSpectralRMSSlope",
&Statistics<dim>::computeSpectralRMSSlope)
.def_static("computeRMSHeights", &Statistics<dim>::computeRMSHeights);
}
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod) {
+ // Exposing logging facility
+ py::enum_<LogLevel>(mod, "LogLevel")
+ .value("debug", LogLevel::debug)
+ .value("info", LogLevel::info)
+ .value("warning", LogLevel::warning)
+ .value("error", LogLevel::error);
+ mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); });
+
// 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,
py::return_value_policy::copy)
.def_static("computePowerSpectrum",
&SurfaceStatistics::computePowerSpectrum,
py::return_value_policy::copy);
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
mod.def("to_voigt",
[](const Grid<Real, 3>& field) {
if (field.getNbComponents() == 9) {
Grid<Real, 3> voigt(field.sizes(), 6);
Loop::loop([](auto in, auto out) { out.symmetrize(in); },
range<MatrixProxy<const Real, 3, 3>>(field),
range<SymMatrixProxy<Real, 3>>(voigt));
return voigt;
} else
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/SConscript b/src/SConscript
index 5b784cf1..570e3625 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,160 +1,161 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
# Copyright (©) 2016-19 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/>.
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
eigenvalues.cpp
+logger.cpp
""".split()
core_list = prepend('core', core_list)
if not main_env['strip_info']:
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
regularized_powerlaw.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
integration/element.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/surface/surface_generator_filter_fft.cpp b/src/core/logger.cpp
similarity index 50%
copy from src/surface/surface_generator_filter_fft.cpp
copy to src/core/logger.cpp
index 4cc35974..3b7fb32d 100644
--- a/src/surface/surface_generator_filter_fft.cpp
+++ b/src/core/logger.cpp
@@ -1,55 +1,62 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "logger.hh"
+#include <boost/preprocessor/seq.hpp>
+#include <boost/preprocessor/stringize.hpp>
+#include <cstdlib>
/* -------------------------------------------------------------------------- */
-#include "surface_generator_filter_fft.hh"
-#include "isopowerlaw.hh"
-#include "surface_statistics.hh"
+namespace tamaas {
/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-extern template class Grid<Real, 2>;
-
-SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT() {
- pl_filter = std::make_shared<Isopowerlaw<2>>();
- generator.setFilter(pl_filter);
+LogLevel Logger::current_level = LogLevel::info;
+
+std::ostream& operator<<(std::ostream& o, const LogLevel& val) {
+ switch (val) {
+#define PRINT_LOGLEVEL(r, data, type) \
+ case data::type: \
+ o << BOOST_PP_STRINGIZE(type); \
+ break;
+ BOOST_PP_SEQ_FOR_EACH(PRINT_LOGLEVEL, LogLevel,
+ (debug)(info)(warning)(error));
+#undef PRINT_LOGLEVEL
+ }
+ return o;
}
-SurfaceGeneratorFilterFFT::~SurfaceGeneratorFilterFFT() {}
-
-Real SurfaceGeneratorFilterFFT::constrainRMS() {
- std::cerr << "Waring /!\\ :: Statistical bias in using the generated RMS as "
- "normalizing factor"
- << std::endl;
- Real hrms_grid = SurfaceStatistics::computeStdev(generator.grid);
- generator.grid *= rms / hrms_grid;
- return rms;
+Logger::~Logger() {
+ if (wish_level >= current_level) {
+ fprintf(stderr, "%s", stream.str().c_str());
+ fflush(stderr);
+ }
}
-Grid<Real, 2>& SurfaceGeneratorFilterFFT::buildSurface() {
- generator.setSizes({{grid_size, grid_size}});
- return generator.buildSurface();
+std::ostream& Logger::get(LogLevel level) {
+ wish_level = level;
+ // std::string name;
+ // stream << level << ": ";
+ return stream;
}
-__END_TAMAAS__
+void Logger::setLevel(LogLevel level) { current_level = level; }
+
+/* -------------------------------------------------------------------------- */
+} // namespace tamaas
diff --git a/src/surface/surface_generator_filter_fft.cpp b/src/core/logger.hh
similarity index 57%
copy from src/surface/surface_generator_filter_fft.cpp
copy to src/core/logger.hh
index 4cc35974..2594d98e 100644
--- a/src/surface/surface_generator_filter_fft.cpp
+++ b/src/core/logger.hh
@@ -1,55 +1,56 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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/>.
*
*/
/* -------------------------------------------------------------------------- */
-
+#ifndef LOGGER_HH
+#define LOGGER_HH
+/* -------------------------------------------------------------------------- */
+#include "tamaas.hh"
+#include <sstream>
/* -------------------------------------------------------------------------- */
-#include "surface_generator_filter_fft.hh"
-#include "isopowerlaw.hh"
-#include "surface_statistics.hh"
+namespace tamaas {
/* -------------------------------------------------------------------------- */
-__BEGIN_TAMAAS__
+/// Log levels enumeration
+enum class LogLevel { debug = 0, info = 1, warning = 2, error = 3 };
-extern template class Grid<Real, 2>;
+/// Logging class
+class Logger {
+public:
+ /// Writing to stderr
+ ~Logger();
-SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT() {
- pl_filter = std::make_shared<Isopowerlaw<2>>();
- generator.setFilter(pl_filter);
-}
+ /// Get stream
+ std::ostream& get(LogLevel level = LogLevel::debug);
-SurfaceGeneratorFilterFFT::~SurfaceGeneratorFilterFFT() {}
+ /// Set acceptable logging level
+ static void setLevel(LogLevel level);
-Real SurfaceGeneratorFilterFFT::constrainRMS() {
- std::cerr << "Waring /!\\ :: Statistical bias in using the generated RMS as "
- "normalizing factor"
- << std::endl;
- Real hrms_grid = SurfaceStatistics::computeStdev(generator.grid);
- generator.grid *= rms / hrms_grid;
- return rms;
-}
+private:
+ static LogLevel current_level;
-Grid<Real, 2>& SurfaceGeneratorFilterFFT::buildSurface() {
- generator.setSizes({{grid_size, grid_size}});
- return generator.buildSurface();
-}
+ std::ostringstream stream;
+ LogLevel wish_level = LogLevel::debug;
+};
+/* -------------------------------------------------------------------------- */
+} // namespace tamaas
-__END_TAMAAS__
+#endif // LOGGER_HH
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index 2e1652d5..4c191288 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,175 +1,174 @@
/**
* @mainpage Welcome to Tamaas !
*
* @section Introduction
* Tamaas is a spectral-integral-equation 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>
* @author Son Pham-Ba <son.phamba@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-19 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#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
} // namespace std
/* -------------------------------------------------------------------------- */
__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;
namespace detail {
template <bool acc, template <typename> class Trait, typename Head,
typename... Tail>
struct fold_trait_tail_rec
: std::integral_constant<bool,
fold_trait_tail_rec<acc and Trait<Head>::value,
Trait, Tail...>::value> {};
template <bool acc, template <typename> class Trait, typename Head>
struct fold_trait_tail_rec<acc, Trait, Head>
: std::integral_constant<bool, acc and Trait<Head>::value> {};
} // namespace detail
template <template <typename> class Trait, typename... T>
struct fold_trait : detail::fold_trait_tail_rec<true, Trait, T...> {};
/* -------------------------------------------------------------------------- */
/// 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(); \
+ sstr << __FILE__ << ":" << __LINE__ << ":FATAL: " << mesg << '\n'; \
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/model/be_engine.cpp b/src/model/be_engine.cpp
index b0f3f1e0..a77d89cf 100644
--- a/src/model/be_engine.cpp
+++ b/src/model/be_engine.cpp
@@ -1,82 +1,85 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "be_engine.hh"
+#include "logger.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
void solve(IntegralOperator::kind kind,
const std::map<IntegralOperator::kind, IntegralOperator*>& operators,
GridBase<Real>& input, GridBase<Real>& output) try {
operators.at(kind)->apply(input, output);
} catch (std::out_of_range& e) {
- std::cerr << "Operator (" << "Westergaard::" << kind << ") not registered\n";
+ Logger().get(LogLevel::warning)
+ << "Operator ("
+ << "Westergaard::" << kind << ") not registered\n";
throw e;
}
template <model_type mtype, IntegralOperator::kind otype>
void registerWestergaardOperator(
std::map<IntegralOperator::kind, IntegralOperator*>& operators,
Model& model) {
std::stringstream sstr;
sstr << "Westergaard::" << otype;
if (operators.find(otype) == operators.end())
operators[otype] =
model.template registerIntegralOperator<Westergaard<mtype, otype>>(
sstr.str());
}
template <model_type type>
void BEEngineTmpl<type>::solveNeumann(GridBase<Real>& neumann,
GridBase<Real>& dirichlet) const {
solve(IntegralOperator::neumann, this->operators, neumann, dirichlet);
}
template <model_type type>
void BEEngineTmpl<type>::solveDirichlet(GridBase<Real>& dirichlet,
GridBase<Real>& neumann) const {
solve(IntegralOperator::dirichlet, this->operators, dirichlet, neumann);
}
template <model_type type>
void BEEngineTmpl<type>::registerNeumann() {
registerWestergaardOperator<type, IntegralOperator::neumann>(this->operators,
*this->model);
}
/// Register dirichlet operator
template <model_type type>
void BEEngineTmpl<type>::registerDirichlet() {
registerWestergaardOperator<type, IntegralOperator::dirichlet>(
this->operators, *this->model);
}
template class BEEngineTmpl<model_type::basic_1d>;
template class BEEngineTmpl<model_type::basic_2d>;
template class BEEngineTmpl<model_type::surface_1d>;
template class BEEngineTmpl<model_type::surface_2d>;
template class BEEngineTmpl<model_type::volume_1d>;
template class BEEngineTmpl<model_type::volume_2d>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/model/model.cpp b/src/model/model.cpp
index 49fd58bd..68110286 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,151 +1,154 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "model.hh"
#include "be_engine.hh"
+#include "logger.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
Model::Model(const std::vector<Real>& system_size,
const std::vector<UInt>& discretization)
: system_size(system_size), discretization(discretization) {}
/* -------------------------------------------------------------------------- */
Model::~Model() = default;
/* -------------------------------------------------------------------------- */
void Model::setElasticity(Real E, Real nu) {
this->E = E;
this->nu = nu;
updateOperators();
}
/* -------------------------------------------------------------------------- */
Real Model::getHertzModulus() const { return E / (1 - nu * nu); }
/* -------------------------------------------------------------------------- */
GridBase<Real>& Model::getTraction() { return getField("traction"); }
const GridBase<Real>& Model::getTraction() const {
return getField("traction");
}
/* -------------------------------------------------------------------------- */
GridBase<Real>& Model::getDisplacement() { return getField("displacement"); }
const GridBase<Real>& Model::getDisplacement() const {
return getField("displacement");
}
/* -------------------------------------------------------------------------- */
const std::vector<Real>& Model::getSystemSize() const { return system_size; }
/* -------------------------------------------------------------------------- */
const std::vector<UInt>& Model::getDiscretization() const {
return discretization;
}
/* -------------------------------------------------------------------------- */
void Model::solveNeumann() {
engine->registerNeumann();
engine->solveNeumann(getField("traction"), getField("displacement"));
}
void Model::solveDirichlet() {
engine->registerDirichlet();
engine->solveDirichlet(getField("displacement"), getField("traction"));
}
/* -------------------------------------------------------------------------- */
IntegralOperator* Model::getIntegralOperator(const std::string& name) {
return operators[name].get();
}
/* -------------------------------------------------------------------------- */
void Model::updateOperators() {
for (auto& op : operators)
op.second->updateFromModel();
}
/* -------------------------------------------------------------------------- */
void Model::addDumper(std::shared_ptr<ModelDumper> dumper) {
this->dumpers.push_back(dumper);
}
void Model::dump() const {
for (auto& dumper : dumpers)
if (dumper)
dumper->dump(*this);
}
/* -------------------------------------------------------------------------- */
void Model::registerField(const std::string& name,
std::shared_ptr<GridBase<Real>> field) {
fields[name] = std::move(field);
}
const GridBase<Real>& Model::getField(const std::string& name) const try {
return *fields.at(name);
} catch (std::out_of_range& e) {
- std::cerr << "Field " << name << " not registered in model\n";
+ Logger().get(LogLevel::warning)
+ << "Field " << name << " not registered in model\n";
throw e;
}
GridBase<Real>& Model::getField(const std::string& name) try {
return *fields.at(name);
} catch (std::out_of_range& e) {
- std::cerr << "Field " << name << " not registered in model\n";
+ Logger().get(LogLevel::warning)
+ << "Field " << name << " not registered in model\n";
throw e;
}
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const Model& _this) {
o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
<< ", nu = " << _this.getPoissonRatio() << ")\n";
// Printing domain size
o << " - domain = [";
std::for_each(_this.getSystemSize().begin(), _this.getSystemSize().end() - 1,
[&o](const Real& x) { o << x << ", "; });
o << _this.getSystemSize().back() << "]\n";
// Printing discretization
o << " - discretization = [";
std::for_each(_this.getDiscretization().begin(),
_this.getDiscretization().end() - 1,
[&o](const UInt& x) { o << x << ", "; });
o << _this.getDiscretization().back() << "]";
return o;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/solvers/beck_teboulle.cpp b/src/solvers/beck_teboulle.cpp
index 903ada04..ce0072e4 100644
--- a/src/solvers/beck_teboulle.cpp
+++ b/src/solvers/beck_teboulle.cpp
@@ -1,107 +1,108 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "beck_teboulle.hh"
+#include "logger.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BeckTeboulle::BeckTeboulle(Model& model, const GridBase<Real>& surface,
Real tolerance, Real mu)
: Kato(model, surface, tolerance, mu) {}
/* -------------------------------------------------------------------------- */
Real BeckTeboulle::solve(GridBase<Real>& g0) {
if (g0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
- "Target mean gap does not have the right number of components");
+ "Target mean gap does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveTmpl<model_type::surface_1d>(g0);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(g0);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real BeckTeboulle::solveTmpl(GridBase<Real>& g0) {
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
*pressure = 0;
GridBase<Real> z(*pressure);
GridBase<Real> pressure_old(*pressure);
Real t0 = 1;
Real t1 = 1;
Real m = 0;
Real lipschitz = engine.getNeumannNorm();
do {
engine.solveNeumann(*pressure, *gap);
addUniform<comp>(*gap, g0);
*gap -= *surfaceComp;
z = *gap;
z /= -lipschitz;
*pressure += z;
enforcePressureCoulomb<comp>();
z = *pressure;
t1 = 1 + std::sqrt(4 * t0 * t0 + 1) / 2;
m = 1 + (t0 - 1) / t1;
z -= pressure_old;
z *= m;
z += pressure_old;
t0 = t1;
pressure_old = *pressure;
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
__END_TAMAAS__
-/* -------------------------------------------------------------------------- */
\ No newline at end of file
+/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/condat.cpp b/src/solvers/condat.cpp
index 582d6cb3..c8015bbd 100644
--- a/src/solvers/condat.cpp
+++ b/src/solvers/condat.cpp
@@ -1,149 +1,152 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "condat.hh"
+#include "logger.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
-Condat::Condat(Model& model, const GridBase<Real>& surface,
- Real tolerance, Real mu)
+Condat::Condat(Model& model, const GridBase<Real>& surface, Real tolerance,
+ Real mu)
: Kato(model, surface, tolerance, mu) {
pressure_old =
- allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
- model.getTraction().getNbComponents());
+ allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
+ model.getTraction().getNbComponents());
}
/* -------------------------------------------------------------------------- */
Real Condat::solve(GridBase<Real>& p0, Real grad_step) {
if (p0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
- "Target mean pressure does not have the right number of components");
+ "Target mean pressure does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveTmpl<model_type::surface_1d>(p0, grad_step);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0, grad_step);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Condat::solveTmpl(GridBase<Real>& p0, Real grad_step) {
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
pressure->uniformSetComponents(p0);
*pressure_old = 0;
- Grid<Real, 1> q({1}, comp); q = 0;
+ Grid<Real, 1> q({1}, comp);
+ q = 0;
Real lipschitz = engine.getNeumannNorm();
Real sigma = 2 * grad_step / lipschitz;
Vector<Real, comp> h_mean = computeMean<comp>(*surfaceComp);
*surfaceComp -= h_mean;
do {
updateGap<comp>(sigma, grad_step, q);
*pressure -= *gap;
enforcePressureCoulomb<comp>();
updateLagrange<comp>(q, p0);
*pressure_old = *pressure;
cost = computeCost();
printState(n, cost, cost);
- } while ((cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations);
+ } while ((cost > this->tolerance || cost == 0.0) &&
+ n++ < this->max_iterations);
*surfaceComp += h_mean;
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
// template <UInt comp>
// void Condat::updateGap(Real sigma, Real grad_step, GridBase<Real>& q) {
// engine.solveNeumann(*pressure, *gap);
// VectorProxy<Real, comp> _q = q(0);
// Vector<Real, comp> qt = _q;
// qt *= (1 - grad_step);
// Loop::stridedLoop(
// [sigma, qt] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
// VectorProxy<Real, comp>&& h) {
// g -= h;
// g *= sigma;
// g += qt;
// },
// *gap, *surfaceComp);
// }
template <UInt comp>
void Condat::updateGap(Real sigma, Real grad_step, GridBase<Real>& q) {
computeGradient<comp>();
VectorProxy<Real, comp> _q(q(0));
Vector<Real, comp> qt = _q;
qt *= (1 - grad_step);
Loop::loop(
[sigma, qt] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
g *= sigma;
g += qt;
},
range<VectorProxy<Real, comp>>(*gap));
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Condat::updateLagrange(GridBase<Real>& q, GridBase<Real>& p0) {
*pressure_old *= -0.5;
*pressure_old += *pressure;
Vector<Real, comp> corr = computeMean<comp>(*pressure_old);
corr *= 2;
q += corr;
q -= p0;
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/contact_solver.cpp b/src/solvers/contact_solver.cpp
index 31bb4fa5..85e9efd0 100644
--- a/src/solvers/contact_solver.cpp
+++ b/src/solvers/contact_solver.cpp
@@ -1,70 +1,70 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "contact_solver.hh"
#include "fft_plan_manager.hh"
+#include "logger.hh"
#include <iomanip>
#include <iostream>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
ContactSolver::ContactSolver(Model& model, const GridBase<Real>& surface,
Real tolerance)
- : model(model), surface(), functional(),
- tolerance(tolerance), surface_stddev(std::sqrt(surface.var())) {
+ : model(model), surface(), functional(), tolerance(tolerance),
+ surface_stddev(std::sqrt(surface.var())) {
auto discretization = model.getBoundaryDiscretization();
auto n_points = model.getTraction().getNbPoints();
if (n_points != surface.dataSize())
TAMAAS_EXCEPTION("Model size and surface size do not match!");
this->surface.wrap(surface);
_gap = allocateGrid<true, Real>(model.getType(), discretization);
-
}
/* -------------------------------------------------------------------------- */
ContactSolver::~ContactSolver() {
FFTPlanManager::get().clean(); /// TODO find a smarter way to clean up
}
/* -------------------------------------------------------------------------- */
void ContactSolver::printState(UInt iter, Real cost_f, Real error) {
// UInt precision = std::cout.precision();
if (iter % dump_frequency)
return;
- std::cout << std::setw(5) << iter << " ";
- std::cout << std::setw(15) << std::scientific << cost_f << " ";
- std::cout << std::setw(15) << error << std::endl;
- std::cout << std::fixed;
+ Logger().get(LogLevel::info)
+ << std::setw(5) << iter << " " << std::setw(15) << std::scientific
+ << cost_f << " " << std::setw(15) << error << std::endl
+ << std::fixed;
}
/* -------------------------------------------------------------------------- */
void ContactSolver::addFunctionalTerm(functional::Functional* func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/epic.cpp b/src/solvers/epic.cpp
index dd594ac9..ecbf0994 100644
--- a/src/solvers/epic.cpp
+++ b/src/solvers/epic.cpp
@@ -1,95 +1,96 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "epic.hh"
+#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
EPICSolver::EPICSolver(ContactSolver& csolver, EPSolver& epsolver,
Real tolerance, Real relaxation)
: csolver(csolver), epsolver(epsolver), tolerance(tolerance),
relaxation(relaxation) {
auto& model = csolver.getModel();
surface.wrap(csolver.getSurface());
pressure.resize(model.getTraction().getNbPoints());
#define SET_VIEWS(_, __, type) \
case type: { \
setViews<type>(); \
break; \
}
switch (model.getType()) {
BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
}
#undef SET_VIEWS
}
/* -------------------------------------------------------------------------- */
void EPICSolver::solve(std::vector<Real> load) {
UInt n = 0, nmax = 1000;
Real error = 0.;
const GridBase<Real> initial_surface(surface);
Real normalizing_factor = std::sqrt(initial_surface.var());
GridBase<Real> previous_residual(*residual_disp);
GridBase<Real> relaxation_direction(*residual_disp);
previous_residual = 0;
do {
csolver.solve(load); // solve contact problem
*pressure_inc -= pressure; // compute pressure increment
epsolver.solve(); // compute residual displacement
relaxation_direction = *residual_disp;
relaxation_direction -= previous_residual;
relaxation_direction *= relaxation;
error = computeError(*residual_disp, previous_residual, normalizing_factor);
previous_residual += relaxation_direction;
surface = initial_surface;
surface -= previous_residual; // update surface
std::cout << "error = " << std::scientific << error << std::fixed
<< std::endl;
} while (error > tolerance && n++ < nmax);
// computing full pressure
pressure += *pressure_inc;
// setting the model pressure to full converged pressure
*pressure_inc = pressure;
// updating plastic state
epsolver.updateState();
}
/* -------------------------------------------------------------------------- */
Real EPICSolver::computeError(const GridBase<Real>& current,
const GridBase<Real>& prev, Real factor) const {
GridBase<Real> error(current);
error -= prev;
return std::sqrt(error.dot(error)) / factor;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp
index 8a0ecf52..62d36129 100644
--- a/src/solvers/kato.cpp
+++ b/src/solvers/kato.cpp
@@ -1,511 +1,502 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "kato.hh"
#include "elastic_functional.hh"
-#include "loop.hh"
#include "fft_plan_manager.hh"
+#include "logger.hh"
+#include "loop.hh"
#include <iomanip>
-#include <iterator>
#include <iostream>
+#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
-Kato::Kato(Model& model, const GridBase<Real>& surface, Real tolerance,
- Real mu)
- : ContactSolver(model, surface, tolerance),
- engine(model.getBEEngine()), mu(mu) {
+Kato::Kato(Model& model, const GridBase<Real>& surface, Real tolerance, Real mu)
+ : ContactSolver(model, surface, tolerance), engine(model.getBEEngine()),
+ mu(mu) {
if (model.getType() != model_type::surface_1d &&
model.getType() != model_type::surface_2d) {
- TAMAAS_EXCEPTION(
- "Model type is not compatible with Kato solver");
+ TAMAAS_EXCEPTION("Model type is not compatible with Kato solver");
}
- gap = this->_gap.get(); // locally allocated
+ gap = this->_gap.get(); // locally allocated
pressure = &model.getTraction();
N = pressure->getNbPoints();
if (model.getType() == model_type::surface_1d) {
initSurfaceWithComponents<model_type::surface_1d>();
} else {
initSurfaceWithComponents<model_type::surface_2d>();
}
engine.registerNeumann();
}
/* -------------------------------------------------------------------------- */
Real Kato::solve(GridBase<Real>& p0, UInt proj_iter) {
if (p0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
- "Target mean pressure does not have the right number of components");
+ "Target mean pressure does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveTmpl<model_type::surface_1d>(p0, proj_iter);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0, proj_iter);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveTmpl(GridBase<Real>& p0, UInt proj_iter) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
pressure->uniformSetComponents(p0);
do {
computeGradient<comp>();
*pressure -= *gap;
enforcePressureConstraints<comp>(p0, proj_iter);
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
Real Kato::solveRelaxed(GridBase<Real>& g0) {
if (g0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
- "Target mean gap does not have the right number of components");
+ "Target mean gap does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveRelaxedTmpl<model_type::surface_1d>(g0);
break;
case model_type::surface_2d:
cost = solveRelaxedTmpl<model_type::surface_2d>(g0);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveRelaxedTmpl(GridBase<Real>& g0) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
*pressure = 0;
do {
engine.solveNeumann(*pressure, *gap);
addUniform<comp>(*gap, g0);
*gap -= *surfaceComp;
*pressure -= *gap;
enforcePressureCoulomb<comp>();
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
Real Kato::solveRegularized(GridBase<Real>& p0, Real r) {
if (p0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
- "Target mean pressure does not have the right number of components");
+ "Target mean pressure does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveRegularizedTmpl<model_type::surface_1d>(p0, r);
break;
case model_type::surface_2d:
cost = solveRegularizedTmpl<model_type::surface_2d>(p0, r);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveRegularizedTmpl(GridBase<Real>& p0, Real r) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
pressure->uniformSetComponents(p0);
do {
// enforcePressureMean<comp>(p0);
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
// Impose zero tangential displacement in non-sliding zone
UInt count_static = 0;
Vector<Real, comp> g_static = Loop::reduce<operation::plus>(
[&] CUDA_LAMBDA(VectorProxy<Real, comp> g,
VectorProxy<Real, comp> p) -> Vector<Real, comp> {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact
count_static++;
return g; // to compute mean of g_T
} else {
return 0;
}
},
range<VectorProxy<Real, comp>>(*gap),
range<VectorProxy<Real, comp>>(*pressure));
g_static /= count_static != 0 ? count_static : 1;
g_static(comp - 1) = 0;
Loop::loop(
[this, r, g_static] CUDA_LAMBDA(VectorProxy<Real, comp> p,
VectorProxy<Real, comp> g) {
// Add frictional term to gradient of functional
g -= g_static;
// Vector<Real, comp> _g = g; // copy
VectorProxy<Real, comp - 1> g_T(g(0));
VectorProxy<Real, 1> g_N(g(comp - 1));
Real g_T_norm = g_T.l2norm();
// g_N += mu * regularize(g_T_norm, r) * g_T_norm;
g_N += mu * g_T_norm;
// Update pressure with gradient
// _g *= 0.1;
p -= g;
// Truncate negative normal pressure
VectorProxy<Real, comp - 1> p_T(p(0));
VectorProxy<Real, 1> p_N(p(comp - 1));
if (p_N(0) < 0)
p_N = 0;
// Set tangential pressure
p_T = g_T;
if (g_T_norm != 0)
p_T *= -mu * p_N(0) * regularize(g_T_norm, r) / g_T_norm;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*gap));
// enforcePressureMean<comp>(p0);
// enforcePressureCoulomb<comp>();
enforcePressureConstraints<comp>(p0, 50);
cost = computeCost();
printState(n, cost, cost);
} while (std::abs(cost) > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Kato::initSurfaceWithComponents() {
constexpr UInt comp = model_type_traits<type>::components;
- surfaceComp =
- allocateGrid<true, Real>(type, model.getDiscretization(), comp);
+ surfaceComp = allocateGrid<true, Real>(type, model.getDiscretization(), comp);
*surfaceComp = 0;
- Loop::loop(
- [] CUDA_LAMBDA(Real& s, VectorProxy<Real, comp> sc) {
- sc(comp - 1) = s;
- },
- surface, range<VectorProxy<Real, comp>>(*surfaceComp));
+ Loop::loop([] CUDA_LAMBDA(Real & s,
+ VectorProxy<Real, comp> sc) { sc(comp - 1) = s; },
+ surface, range<VectorProxy<Real, comp>>(*surfaceComp));
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureMean(GridBase<Real>& p0) {
Vector<Real, comp> corr = computeMean<comp>(*pressure);
VectorProxy<Real, comp> _p0(p0(0));
corr -= _p0;
*pressure -= corr;
}
/* -------------------------------------------------------------------------- */
// template <UInt comp>
// void Kato::enforcePressureCoulomb() {
// Loop::stridedLoop(
// [this] CUDA_LAMBDA(VectorProxy<Real, comp>&& p) {
// VectorProxy<Real, comp - 1> p_T(p(0));
// Real p_N = p(comp - 1);
// Real p_T_sqrd= p_T.l2squared();
// // Projection normale au cône de friction
// bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N);
// bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu);
// if (cond2) {
// p_T = 0;
// p(comp - 1) = 0;
// } else if (!cond1) {
// Real p_T_norm = std::sqrt(p_T_sqrd);
// Real k = (p_N + mu * p_T_norm) / (1 + mu * mu);
// p_T *= k * mu / p_T_norm;
// p(comp - 1) = k;
// }
// },
// *pressure);
// }
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureTresca() {
Loop::loop(
[this] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
if (p_N < 0) {
p_T = 0;
p(comp - 1) = 0;
} else if (p_T_norm > mu) { // TODO: replace name of variable mu
p_T *= mu / p_T_norm;
// p(comp - 1) += p_T_norm - mu;
}
},
range<VectorProxy<Real, comp>>(*pressure));
}
template void Kato::enforcePressureTresca<2>();
template void Kato::enforcePressureTresca<3>();
/* -------------------------------------------------------------------------- */
/**
* Compute mean of the field taking each component separately.
*/
// template <UInt comp>
// Vector<Real, comp> Kato::computeMean(GridBase<Real>& field) {
// Vector<Real, comp> mean = Loop::stridedReduce<operation::plus>(
// [] CUDA_LAMBDA(VectorProxy<Real, comp>&& f) -> Vector<Real, comp> {
// return f;
// },
// field);
// mean /= N;
// return mean;
// }
/* -------------------------------------------------------------------------- */
// template <UInt comp>
// void Kato::addUniform(GridBase<Real>& field, GridBase<Real>& vec) {
// VectorProxy<Real, comp> _vec(vec(0));
// field += _vec;
// }
/* -------------------------------------------------------------------------- */
Real Kato::computeCost(bool use_tresca) {
UInt N = pressure->getNbPoints();
Grid<Real, 1> lambda({N}, 1);
Grid<Real, 1> eta({N}, 1);
Grid<Real, 1> p_N({N}, 1);
Grid<Real, 1> p_C({N}, 1);
switch (model.getType()) {
case model_type::surface_1d:
if (!use_tresca) {
computeValuesForCost<model_type::surface_1d>(lambda, eta, p_N, p_C);
} else {
computeValuesForCostTresca<model_type::surface_1d>(lambda, eta, p_N, p_C);
}
break;
case model_type::surface_2d:
if (!use_tresca) {
computeValuesForCost<model_type::surface_2d>(lambda, eta, p_N, p_C);
} else {
computeValuesForCostTresca<model_type::surface_2d>(lambda, eta, p_N, p_C);
}
break;
default:
break;
}
return p_N.dot(lambda) + p_C.dot(eta);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Kato::computeValuesForCost(GridBase<Real>& lambda, GridBase<Real>& eta,
- GridBase<Real>& p_N, GridBase<Real>& p_C) {
+ GridBase<Real>& p_N, GridBase<Real>& p_C) {
constexpr UInt comp = model_type_traits<type>::components;
Real g_N_min = Loop::reduce<operation::min>(
- [] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
- return g(comp - 1);
- },
- range<VectorProxy<Real, comp>>(*gap));
+ [] CUDA_LAMBDA(VectorProxy<Real, comp> g) { return g(comp - 1); },
+ range<VectorProxy<Real, comp>>(*gap));
Loop::loop(
- [this, g_N_min] CUDA_LAMBDA(
- VectorProxy<Real, comp> p, VectorProxy<Real, comp> g, Real & lambda_,
- Real & eta_, Real & p_N_, Real & p_C_) {
+ [this, g_N_min] CUDA_LAMBDA(VectorProxy<Real, comp> p,
+ VectorProxy<Real, comp> g, Real & lambda_,
+ Real & eta_, Real & p_N_, Real & p_C_) {
VectorProxy<Real, comp - 1> g_T(g(0));
Real g_N = g(comp - 1);
Real g_T_norm = g_T.l2norm();
lambda_ = g_N - g_N_min;
eta_ = g_T_norm;
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
p_N_ = p_N;
p_C_ = (p_N > 0) ? mu * p_N - p_T_norm : 0;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*gap), lambda, eta, p_N, p_C);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
-void Kato::computeValuesForCostTresca(GridBase<Real>& lambda, GridBase<Real>& eta,
- GridBase<Real>& p_N, GridBase<Real>& p_C) {
+void Kato::computeValuesForCostTresca(GridBase<Real>& lambda,
+ GridBase<Real>& eta, GridBase<Real>& p_N,
+ GridBase<Real>& p_C) {
constexpr UInt comp = model_type_traits<type>::components;
Real g_N_min = Loop::reduce<operation::min>(
- [] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
- return g(comp - 1);
- },
- range<VectorProxy<Real, comp>>(*gap));
+ [] CUDA_LAMBDA(VectorProxy<Real, comp> g) { return g(comp - 1); },
+ range<VectorProxy<Real, comp>>(*gap));
Real p_C_min = Loop::reduce<operation::min>(
- [this] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
- VectorProxy<Real, comp - 1> p_T(p(0));
- Real p_T_norm = p_T.l2norm();
- return mu - p_T_norm;
- },
- range<VectorProxy<Real, comp>>(*pressure));
+ [this] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
+ VectorProxy<Real, comp - 1> p_T(p(0));
+ Real p_T_norm = p_T.l2norm();
+ return mu - p_T_norm;
+ },
+ range<VectorProxy<Real, comp>>(*pressure));
Loop::loop(
[this, g_N_min, p_C_min] CUDA_LAMBDA(
VectorProxy<Real, comp> p, VectorProxy<Real, comp> g, Real & lambda_,
Real & eta_, Real & p_N_, Real & p_C_) {
VectorProxy<Real, comp - 1> g_T(g(0));
Real g_N = g(comp - 1);
Real g_T_norm = g_T.l2norm();
lambda_ = g_N - g_N_min;
eta_ = g_T_norm;
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
p_N_ = p_N;
p_C_ = (p_N > 0) ? mu - p_T_norm - p_C_min : 0;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*gap), lambda, eta, p_N, p_C);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::computeFinalGap() {
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
Real g_N_min = Loop::reduce<operation::min>(
- [] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
- return g(comp - 1);
- },
- range<VectorProxy<Real, comp>>(*gap));
+ [] CUDA_LAMBDA(VectorProxy<Real, comp> g) { return g(comp - 1); },
+ range<VectorProxy<Real, comp>>(*gap));
Grid<Real, 1> g_shift({comp}, 1);
g_shift = 0;
g_shift(comp - 1) = -g_N_min;
*gap += *surfaceComp;
addUniform<comp>(*gap, g_shift);
model.getDisplacement() = *gap;
}
/* -------------------------------------------------------------------------- */
Real Kato::regularize(Real x, Real r) {
Real xr = x / r;
return xr / (1 + std::abs(xr));
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp
index 1a3d582e..eed2ea75 100644
--- a/src/solvers/polonsky_keer_rey.cpp
+++ b/src/solvers/polonsky_keer_rey.cpp
@@ -1,292 +1,293 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "polonsky_keer_rey.hh"
#include "elastic_functional.hh"
#include "fft_plan_manager.hh"
+#include "logger.hh"
#include "loop.hh"
#include "model_type.hh"
#include <boost/preprocessor/seq.hpp>
#include <iomanip>
#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase<Real>& surface,
Real tolerance, type variable_type,
type constraint_type)
: ContactSolver(model, surface, tolerance), variable_type(variable_type),
constraint_type(constraint_type) {
#define SET_VIEWS(_, __, type) \
case type: { \
setViews<type>(); \
break; \
}
switch (model.getType()) {
BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
}
#undef SET_VIEWS
search_direction = allocateGrid<true, Real>(
operation_type, model.getBoundaryDiscretization());
projected_search_direction = allocateGrid<true, Real>(
operation_type, model.getBoundaryDiscretization());
switch (variable_type) {
case gap: {
model.getBEEngine().registerDirichlet();
primal = gap_view.get();
dual = pressure_view.get();
this->functional.addFunctionalTerm(
new functional::ElasticFunctionalGap(*integral_op, this->surface),
true);
break;
}
case pressure: {
model.getBEEngine().registerNeumann();
primal = pressure_view.get();
dual = gap_view.get();
this->functional.addFunctionalTerm(
new functional::ElasticFunctionalPressure(*integral_op, this->surface),
true);
break;
}
}
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerRey::solve(std::vector<Real> target_v) {
const Real target = target_v.back();
Real G = 0, Gold = 1, error = 0, error_norm = 1;
UInt n = 0;
bool delta = false;
*search_direction = 0;
//*dual = 0;
// Printing column headers
- std::cout << std::setw(5) << "Iter"
- << " " << std::setw(15) << "Cost_f"
- << " " << std::setw(15) << "Error" << '\n'
- << std::fixed;
+ Logger().get(LogLevel::info) << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
// Setting uniform value if constraint
if (constraint_type == variable_type && std::abs(primal->sum()) <= 0)
*primal = target;
else if (constraint_type == variable_type)
*primal *= target / primal->mean();
else if (constraint_type != variable_type)
*primal = this->surface_stddev;
do {
// Computing functional gradient
functional.computeGradF(*primal, *dual);
Real dbar = meanOnUnsaturated(*dual);
// Enforcing dual constraint via gradient
if (constraint_type != variable_type) {
*dual += 2 * target + dbar;
} else {
// Centering dual on primal > 0
*dual -= dbar;
}
// Computing gradient norm
G = computeSquaredNorm(*dual);
// Updating search direction (conjugate gradient)
updateSearchDirection(delta * G / Gold);
Gold = G;
// Computing critical step
Real tau = computeCriticalStep(target);
// Update primal variable
delta = updatePrimal(tau);
// We should scale to constraint
if (constraint_type == variable_type)
*primal *= target / primal->mean();
error = computeError() / error_norm;
Real cost_f = functional.computeF(*primal, *dual);
printState(n, cost_f, error);
} while (error > this->tolerance && n++ < this->max_iterations);
// Final update of dual variable
functional.computeGradF(*primal, *dual);
Real dual_min = dual->min();
*dual -= dual_min;
// Primal is pressure: need to make sure gap is positive
if (variable_type == pressure) {
*displacement_view = *dual;
*displacement_view += this->surface;
}
// Primal is gap: need to make sure dual is positive (pressure + adhesion)
else {
*displacement_view = *primal;
*displacement_view += this->surface;
integral_op->apply(*displacement_view, *pressure_view);
*pressure_view -= dual_min;
}
return error;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$
*/
Real PolonskyKeerRey::meanOnUnsaturated(const GridBase<Real>& field) const {
// Sum on unsaturated
Real sum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
*primal, field);
// Number of unsaturated points
UInt n_unsat = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal);
return sum / n_unsat;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
*/
Real PolonskyKeerRey::computeSquaredNorm(const GridBase<Real>& field) const {
return Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; },
*primal, field);
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i'
* t_i} } \f$
*/
Real PolonskyKeerRey::computeCriticalStep(Real target) {
integral_op->apply(*search_direction, *projected_search_direction);
Real rbar = meanOnUnsaturated(*projected_search_direction);
if (constraint_type == variable_type)
*projected_search_direction -= rbar;
else {
*projected_search_direction += 2 * target + rbar;
}
Real num = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(const Real& p, const Real& q, const Real& t) {
return (p > 0) ? q * t : 0;
},
*primal, *dual, *search_direction);
Real denum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(const Real& p, const Real& r, const Real& t) {
return (p > 0) ? r * t : 0;
},
*primal, *projected_search_direction, *search_direction);
return num / denum;
}
/* -------------------------------------------------------------------------- */
/**
* Update steps:
* 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$
* 2. Truncate all \f$p\f$ negative
* 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i =
* p_i - \tau q_i\f$
*/
bool PolonskyKeerRey::updatePrimal(Real step) {
UInt na_num = Loop::reduce<operation::plus>(
[step] CUDA_LAMBDA(Real & p, const Real& q, const Real& t) -> UInt {
p -= step * t; // Updating primal
if (p < 0)
p = 0; // Truncating neg values
if (p == 0 && q < 0) { // If non-admissible state
p -= step * q;
return 1;
} else
return 0;
},
*primal, *dual, *search_direction);
return na_num == 0;
}
/* -------------------------------------------------------------------------- */
/**
* Error is based on \f$ \sum{p_i q_i} \f$
*/
Real PolonskyKeerRey::computeError() {
/// Making sure dual respects constraint
*dual -= dual->min();
Real norm = 1;
Real error = primal->dot(*dual);
if (std::isnan(error))
TAMAAS_EXCEPTION("Encountered NaN in complementarity error: this may be "
"caused by a contact area of a single node.");
if (variable_type == pressure)
norm = std::abs(primal->sum() * this->surface_stddev);
else
norm = std::abs(dual->sum() * this->surface_stddev);
norm *= primal->getNbPoints();
return std::abs(error) / norm;
}
/* -------------------------------------------------------------------------- */
/**
* Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t}
* \f$
*/
void PolonskyKeerRey::updateSearchDirection(Real factor) {
Loop::loop(
[factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
t = (p > 0) ? q + factor * t : 0;
},
*primal, *dual, *search_direction);
}
/* -------------------------------------------------------------------------- */
void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/surface/surface_generator_filter_fft.cpp b/src/surface/surface_generator_filter_fft.cpp
index 4cc35974..ddac0fb9 100644
--- a/src/surface/surface_generator_filter_fft.cpp
+++ b/src/surface/surface_generator_filter_fft.cpp
@@ -1,55 +1,57 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "surface_generator_filter_fft.hh"
#include "isopowerlaw.hh"
+#include "logger.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
extern template class Grid<Real, 2>;
SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT() {
pl_filter = std::make_shared<Isopowerlaw<2>>();
generator.setFilter(pl_filter);
}
SurfaceGeneratorFilterFFT::~SurfaceGeneratorFilterFFT() {}
Real SurfaceGeneratorFilterFFT::constrainRMS() {
- std::cerr << "Waring /!\\ :: Statistical bias in using the generated RMS as "
- "normalizing factor"
- << std::endl;
+ Logger().get(LogLevel::warning)
+ << "Statistical bias in using the generated RMS as "
+ "normalizing factor"
+ << std::endl;
Real hrms_grid = SurfaceStatistics::computeStdev(generator.grid);
generator.grid *= rms / hrms_grid;
return rms;
}
Grid<Real, 2>& SurfaceGeneratorFilterFFT::buildSurface() {
generator.setSizes({{grid_size, grid_size}});
return generator.buildSurface();
}
__END_TAMAAS__