diff --git a/src/model/boussinesq.cpp b/src/model/boussinesq.cpp
index aa1112e..42849db 100644
--- a/src/model/boussinesq.cpp
+++ b/src/model/boussinesq.cpp
@@ -1,95 +1,97 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "boussinesq.hh"
#include "boussinesq_helper.hh"
#include "influence.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <>
Boussinesq::Boussinesq(Model* model) : parent(model) {
this->initialize(trait::dimension, trait::dimension);
}
template <>
Boussinesq::Boussinesq(Model* model) : parent(model) {
this->initialize(trait::dimension, trait::voigt);
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template
void Boussinesq::apply(GridBase& source,
GridBase& out) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
influence::Boussinesq boussinesq(mu, nu);
influence::ElasticHelper el(mu, nu);
detail::BoussinesqHelper helper;
// the source is a Grid grid normally, so we gotta make sure we
// can cast it into a Grid
GridView view(
source, {}, -1);
parent::transformSource(view);
helper.template apply(this->source_buffer.front(), this->out_buffer,
this->wavevectors,
this->model->getSystemSize().front(), boussinesq);
helper.makeFundamentalModeGreatAgain(this->source_buffer.front(),
this->out_buffer, el);
parent::transformOutput(
- [](auto&& out_buffer, auto layer) { return out_buffer[layer]; }, out);
+ [](auto&& out_buffer, auto layer) ->
+ typename parent::BufferType& { return out_buffer[layer]; },
+ out);
}
/* -------------------------------------------------------------------------- */
template
void Boussinesq::initialize(UInt source_components,
UInt out_components) {
auto hermitian_sizes =
GridHermitian::hermitianDimensions(
this->model->getBoundaryDiscretization());
// Initializing buffers
this->source_buffer.resize(1);
this->out_buffer.resize(this->model->getDiscretization()[0]);
// Resizing source buffer
std::for_each(this->source_buffer.begin(), this->source_buffer.end(),
[&](auto& buffer) {
buffer.setNbComponents(source_components);
buffer.resize(hermitian_sizes);
});
// Resizing output buffer
std::for_each(this->out_buffer.begin(), this->out_buffer.end(),
[&](auto& buffer) {
buffer.setNbComponents(out_components);
buffer.resize(hermitian_sizes);
});
}
} // namespace tamaas
diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index d779937..a950bf8 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,138 +1,141 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "kelvin.hh"
#include "elasto_plastic_model.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
template
Kelvin::Kelvin(Model* model) : VolumePotential(model) {
setIntegrationMethod(integration_method::linear, 1e-12);
}
template
void Kelvin::setIntegrationMethod(integration_method method,
Real cutoff) {
this->method = method;
this->cutoff = cutoff;
if (this->method == integration_method::linear) {
Logger().get(LogLevel::debug)
<< "[Kelvin " << this << "] Setting linear integration method\n";
this->initialize(dtrait::template source_components,
dtrait::template out_components,
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,
dtrait::template out_components, 1);
}
auto max_q = Loop::reduce(
[](auto qv) { return qv.l2norm(); },
range>(
this->wavevectors));
if (this->method == integration_method::linear and
not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
TAMAAS_EXCEPTION("Probable overflow of integral computation (change "
"integration method to integration_method::cutoff)");
}
/* --------------------------------------------------------------------------
*/
/* Operator implementation */
/* --------------------------------------------------------------------------
*/
template
void Kelvin::applyIf(GridBase& source,
GridBase& 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
void Kelvin::linearIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper 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) { return out_buffer[layer]; }, out);
+ [](auto&& out_buffer, auto layer) ->
+ typename parent::BufferType& { return out_buffer[layer]; },
+ out);
}
template
void Kelvin::cutoffIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper helper;
- auto func = [&](auto&& out_buffer, auto layer) {
+ 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;
template class Kelvin;
template class Kelvin;
/* --------------------------------------------------------------------------
*/
} // namespace tamaas
diff --git a/src/model/mindlin.cpp b/src/model/mindlin.cpp
index 105de46..f715624 100644
--- a/src/model/mindlin.cpp
+++ b/src/model/mindlin.cpp
@@ -1,160 +1,163 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "mindlin.hh"
#include "boussinesq_helper.hh"
#include "elasto_plastic_model.hh"
#include "influence.hh"
#include "kelvin_helper.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
Mindlin::Mindlin(Model* model) : parent(model) {
surface_tractions.setNbComponents(trait::components);
surface_tractions.resize(this->source_buffer.front().sizes());
}
/* -------------------------------------------------------------------------- */
template
void Mindlin::applyIf(GridBase& source,
GridBase& out,
filter_t pred) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
influence::Kelvin kelvin_strain(mu, nu);
influence::ElasticHelper elasticity(mu, nu);
detail::SurfaceTractionHelper thelper;
const Real L = this->model->getSystemSize().front();
parent::transformSource(source, pred);
// Reset buffer values
for (auto&& layer : this->out_buffer)
layer = 0;
// computing surface tractions (q power cancelled with
// boussinesq)
surface_tractions = 0;
thelper.template computeSurfaceTractions(
this->source_buffer, surface_tractions, this->wavevectors, L,
kelvin_strain, elasticity);
surface_tractions *= -1;
if (this->method == integration_method::linear)
linearIntegral(out);
else
cutoffIntegral(out);
}
template
void Mindlin::linearIntegral(GridBase& out) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
Real L = this->model->getSystemSize().front();
influence::Kelvin kelvin(mu, nu);
influence::Boussinesq boussinesq(mu, nu);
influence::ElasticHelper elasticity(mu, nu);
detail::KelvinHelper helper;
detail::SurfaceTractionHelper thelper;
detail::BoussinesqHelper bhelper;
// apply kelvin
helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
L, kelvin);
helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
// no need for fundamental correction
// apply boussinesq (q power cancelled with mindlin 2-grad)
bhelper.template apply(surface_tractions, this->out_buffer,
this->wavevectors, L, boussinesq);
// Correcting for fundamental mode
Vector n{{{0, 0, -1}}};
for (UInt i : Loop::range(this->source_buffer.size())) {
typename detail::KelvinTrait::source_t w(
this->source_buffer[i](0));
typename detail::KelvinTrait::out_t u(
this->out_buffer[i](0));
auto t = dense(w) * n;
bhelper.makeFundamentalModeGreatAgain(t, u, elasticity);
}
parent::transformOutput(
- [](auto&& out_buffer, auto layer) { return out_buffer[layer]; }, out);
+ [](auto&& out_buffer, auto layer) ->
+ typename parent::BufferType& { return out_buffer[layer]; },
+ out);
}
template
void Mindlin::cutoffIntegral(GridBase& out) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
Real L = this->model->getSystemSize().front();
influence::Kelvin kelvin(mu, nu);
influence::Boussinesq boussinesq(mu, nu);
influence::ElasticHelper elasticity(mu, nu);
detail::KelvinHelper helper;
detail::SurfaceTractionHelper thelper;
detail::BoussinesqHelper bhelper;
Vector n{{{0, 0, -1}}};
- auto func = [&](auto&& out_buffer, auto layer) {
+ auto func = [&](auto&& out_buffer, auto layer) ->
+ typename parent::BufferType& {
auto&& out = out_buffer.front();
// Reset previous values
out = 0;
// Apply kelvin
helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors, L,
this->cutoff, kelvin);
helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
helper.makeFundamentalGreatAgain(out);
// Apply boussinesq to cancel surface tractions
bhelper.template apply(surface_tractions, out, layer,
this->wavevectors, this->source_buffer.size(),
L, boussinesq);
// Take care of fundamental mode
typename detail::KelvinTrait::source_t w(
this->source_buffer[layer](0));
typename detail::KelvinTrait::out_t u(out(0));
auto t = dense(w) * n;
bhelper.makeFundamentalModeGreatAgain(t, u, elasticity);
return out;
};
parent::transformOutput(func, out);
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Mindlin;
template class Mindlin;
} // namespace tamaas
diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh
index 3462951..ca3758d 100644
--- a/src/model/volume_potential.hh
+++ b/src/model/volume_potential.hh
@@ -1,142 +1,142 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef VOLUME_POTENTIAL_HH
#define VOLUME_POTENTIAL_HH
/* -------------------------------------------------------------------------- */
#include "fft_plan_manager.hh"
#include "grid_hermitian.hh"
#include "grid_view.hh"
#include "integral_operator.hh"
#include "logger.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Trait type for component management
template
struct derivative_traits;
template <>
struct derivative_traits<0> {
template
static constexpr auto source_components = model_type_traits::components;
template
static constexpr auto out_components = model_type_traits::components;
};
template <>
struct derivative_traits<1> {
template
static constexpr auto source_components = model_type_traits::voigt;
template
static constexpr auto out_components = model_type_traits::components;
};
template <>
struct derivative_traits<2> {
template
static constexpr auto source_components = model_type_traits::voigt;
template
static constexpr auto out_components = model_type_traits::voigt;
};
/**
* @brief Volume potential operator class. Applies the operators for computation
* of displacements and strains due to residual/eigen strains
*/
template
class VolumePotential : public IntegralOperator {
using trait = model_type_traits;
protected:
using filter_t = std::function;
public:
VolumePotential(Model* model);
/// Update from model (does nothing)
void updateFromModel() override {}
/// Kind
IntegralOperator::kind getKind() const override {
return IntegralOperator::neumann;
}
/// Type
model_type getType() const override;
/// Apply to all of the source layers
void apply(GridBase& input, GridBase& output) const override {
applyIf(input, output, [](UInt) { return true; });
}
protected:
/// Transform source layer-by-layer
void transformSource(GridBase& in, filter_t pred) const;
/// Transform all source
void transformSource(GridBase& in) const {
transformSource(in, [](auto) { return true; });
}
/// Transform output layer-by-layer
template
void transformOutput(Func func, GridBase& out) const;
/// Initialize fourier buffers
void initialize(UInt source_components, UInt out_components,
UInt source_buffer_size);
protected:
Grid wavevectors;
using BufferType = GridHermitian;
mutable std::vector source_buffer;
mutable std::vector out_buffer;
};
/* -------------------------------------------------------------------------- */
/* Template implementation */
/* -------------------------------------------------------------------------- */
template
template
void VolumePotential::transformOutput(Func func,
GridBase& out) const {
constexpr UInt dim = trait::dimension;
auto& o = dynamic_cast&>(out);
// Transforming output
for (UInt layer : Loop::range(o.sizes().front())) {
Logger().get(LogLevel::debug)
<< "[VolumePotential] Computing layer " << layer << '\n';
auto out_layer = make_view(o, layer);
- auto&& fourier_out_layer = func(out_buffer, layer);
+ auto& fourier_out_layer = func(out_buffer, layer);
FFTPlanManager::get()
.createPlan(out_layer, fourier_out_layer)
.backwardTransform();
}
}
} // namespace tamaas
#endif // VOLUME_POTENTIAL_HH