Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/kelvin.cpp
4281604d9e17features/netcdf-mpi

Last Change
Raw File

kelvin.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 "kelvin.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
Kelvin<type, derivative>::Kelvin(Model* model) : VolumePotential<type>(model) {
setIntegrationMethod(integration_method::linear, 1e-12);
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::setIntegrationMethod(integration_method method,
Real cutoff) {
this->method = method;
this->cutoff = cutoff;
Logger logger;
if (this->method == integration_method::linear) {
logger.get(LogLevel::debug)
<< "[Kelvin " << this << "] Setting linear integration method\n";
this->initialize(dtrait::template source_components<type>,
dtrait::template out_components<type>,
this->model->getDiscretization()[0]);
}
else {
logger.get(LogLevel::debug)
<< "[Kelvin " << this << "] Setting cutoff integration method (cutoff "
<< this->cutoff << ")\n";
this->initialize(dtrait::template source_components<type>,
dtrait::template out_components<type>, 1);
}
auto max_q = Loop::reduce<operation::max>(
[](auto qv) { return qv.l2norm(); },
range<VectorProxy<const Real, trait::boundary_dimension>>(
this->wavevectors));
if (this->method == integration_method::linear and
not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
logger.get(LogLevel::warning)
<< "[Kelvin " << this
<< "] Probable overflow of integral computation (consider "
"changing integration method to integration_method::cutoff or "
"compiling with real_type='long double')\n";
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::applyIf(GridBase<Real>& source,
GridBase<Real>& out,
filter_t pred) const {
KelvinInfluence kelvin(this->model->getShearModulus(),
this->model->getPoissonRatio());
parent::transformSource(source, pred);
// Reset buffer values
for (auto&& layer : this->out_buffer)
layer = 0;
if (method == integration_method::linear) {
linearIntegral(out, kelvin);
} else {
cutoffIntegral(out, kelvin);
}
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::linearIntegral(GridBase<Real>& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper<type, KelvinInfluence> helper;
helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
this->model->getSystemSize().front(), kelvin);
helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
helper.makeFundamentalGreatAgain(this->out_buffer);
parent::transformOutput(
[](auto&& out_buffer, auto layer) ->
typename parent::BufferType& { return out_buffer[layer]; },
out);
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::cutoffIntegral(GridBase<Real>& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper<type, KelvinInfluence> helper;
auto func = [&](auto&& out_buffer, auto layer) ->
typename parent::BufferType& {
auto&& out = out_buffer.front();
out = 0;
helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors,
this->model->getSystemSize().front(), cutoff, kelvin);
helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
helper.makeFundamentalGreatAgain(out);
return out;
};
parent::transformOutput(func, out);
}
/* --------------------------------------------------------------------------
*/
/* Template instanciation */
/* --------------------------------------------------------------------------
*/
template class Kelvin<model_type::volume_2d, 0>;
template class Kelvin<model_type::volume_2d, 1>;
template class Kelvin<model_type::volume_2d, 2>;
/* --------------------------------------------------------------------------
*/
} // namespace tamaas
c4science · Help