Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/boussinesq.cpp
e4452faac64ffeatures/netcdf-mpi

Last Change
Raw File

boussinesq.cpp
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boussinesq.hh"
#include "influence.hh"
#include "kelvin_helper.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace detail {
template <model_type type, typename T>
class BoussinesqHelper {
using boussinesq_t = T;
using trait = model_type_traits<type>;
static constexpr UInt dim = trait::dimension;
static constexpr UInt bdim = trait::boundary_dimension;
public:
BoussinesqHelper(const T& boussinesq, Real y_3, Real cutoff)
: boussinesq(boussinesq), y_3(y_3), cutoff(cutoff) {}
inline void applyIntegral(GridHermitian<Real, bdim>& out,
GridHermitian<Real, bdim>& source,
const Grid<Real, bdim>& wavevectors) const {
Loop::stridedLoop(
[&](typename KelvinTrait<T>::out_t&& out_local,
VectorProxy<Complex, dim>&& traction,
VectorProxy<const Real, bdim>&& q) {
// Cutoff
if (-q.l2norm() * std::abs(y_3) < std::log(cutoff))
return;
out_local +=
boussinesq.applyU0(traction, q) *
influence::KelvinIntegrator<0>::g0<true>(q.l2norm() * y_3);
out_local +=
(boussinesq.applyU1(traction, q) *
influence::KelvinIntegrator<0>::g1<true>(q.l2norm() * y_3));
},
out, source, wavevectors);
}
inline void applyConstantTerm(typename KelvinTrait<T>::out_t& out,
typename KelvinTrait<T>::source_t& source,
const influence::ElasticHelper<dim>& el) const {
out = 0;
}
protected:
const T& boussinesq;
const Real y_3, cutoff;
};
template <>
inline void
BoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
applyConstantTerm(
typename KelvinTrait<influence::Boussinesq<3, 1>>::out_t& out,
typename KelvinTrait<influence::Boussinesq<3, 1>>::source_t& source,
const influence::ElasticHelper<dim>& el) const {
out = 0;
out(2, 0) = -source(0) / el.mu;
out(2, 1) = -source(1) / el.mu;
out(2, 2) = -source(2) / (el.lambda + 2 * el.mu);
}
} // namespace detail
template <>
Boussinesq<model_type::volume_2d, 0>::Boussinesq(Model* model) : parent(model) {
this->initialize(trait::dimension, trait::dimension);
}
template <>
Boussinesq<model_type::volume_2d, 1>::Boussinesq(Model* model) : parent(model) {
this->initialize(trait::dimension, trait::dimension * trait::dimension);
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
void Boussinesq<type, derivative>::apply(GridBase<Real>& source,
GridBase<Real>& out) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
influence::Boussinesq<trait::dimension, derivative> boussinesq(mu, nu);
auto apply = [&](UInt i, decltype(this->source_buffers)& source_buffers,
decltype(this->disp_buffer)& out_buffer) {
const Real L = this->model->getSystemSize().front();
const UInt N = this->model->getDiscretization().front();
const Real dl = L / (N - 1);
out_buffer = 0;
const Real xi = i * dl;
detail::BoussinesqHelper<type, decltype(boussinesq)> helper(boussinesq, xi,
1e-2);
helper.applyIntegral(out_buffer, source_buffers.front(), this->wavevectors);
typename detail::KelvinTrait<decltype(boussinesq)>::out_t out_fundamental(
out_buffer(0));
typename detail::KelvinTrait<decltype(boussinesq)>::source_t in_fundamental(
source_buffers.front()(0));
// Uniform shift in case of gradient computation
helper.applyConstantTerm(out_fundamental, in_fundamental,
influence::ElasticHelper<trait::dimension>(mu, nu));
};
this->fourierApply(apply, source, out);
}
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
void Boussinesq<type, derivative>::initialize(UInt source_components,
UInt out_components) {
// Copy horizontal sizes
std::array<UInt, trait::boundary_dimension> sizes;
std::copy(this->model->getDiscretization().begin() + 1,
this->model->getDiscretization().end(), sizes.begin());
auto hermitian_sizes =
GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
sizes);
this->disp_buffer.setNbComponents(out_components);
this->disp_buffer.resize(hermitian_sizes);
this->source_buffers.resize(1);
auto& source = this->source_buffers.front();
source.setNbComponents(source_components);
source.resize(hermitian_sizes);
}
__END_TAMAAS__
c4science · Help