Homec4science
Diffusion tamaas (feature/mpi_fft_engine) Browse

tamaas/src/model/boussinesq_helper.hh
02240a9f1952feature/mpi_fft_engine

Last Change
Raw File

boussinesq_helper.hh
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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BOUSSINESQ_HELPER_HH
#define BOUSSINESQ_HELPER_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "grid_hermitian.hh"
#include "influence.hh"
#include "integration/accumulator.hh"
#include "kelvin_helper.hh"
#include "model.hh"
#include "model_type.hh"
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace detail {
template <model_type type, typename boussinesq_t>
struct BoussinesqHelper {
using trait = model_type_traits<type>;
static constexpr UInt dim = trait::dimension;
static constexpr UInt bdim = trait::boundary_dimension;
using BufferType = GridHermitian<Real, bdim>;
using source_t = typename KelvinTrait<boussinesq_t>::source_t;
using out_t = typename KelvinTrait<boussinesq_t>::out_t;
template <bool apply_q_power>
void apply(BufferType& tractions, std::vector<BufferType>& out,
const Grid<Real, bdim>& wavevectors, Real domain_size,
const boussinesq_t& boussinesq) {
for (UInt l : Loop::range(out.size()))
apply<apply_q_power>(tractions, out[l], l, wavevectors, out.size(),
domain_size, boussinesq);
}
template <bool apply_q_power>
void apply(BufferType& tractions, BufferType& out, UInt layer,
const Grid<Real, bdim>& wavevectors, UInt discretization,
Real domain_size, const boussinesq_t& boussinesq) {
accumulator.makeUniformMesh(discretization, domain_size);
const auto xl = accumulator.nodePositions()[layer];
Loop::loop(
[xl, &boussinesq](auto qv, auto u, auto t) {
const auto q = qv.l2norm();
const auto e = std::exp(-q * xl);
// g0 part
auto tmp = boussinesq.template applyU0<apply_q_power>(t, qv);
tmp *= e;
u += tmp;
// g1 part
tmp = boussinesq.template applyU1<apply_q_power>(t, qv);
tmp *= q * xl * e;
u += tmp;
},
range<VectorProxy<const Real, bdim>>(wavevectors).headless(),
range<out_t>(out).headless(), range<source_t>(tractions).headless());
}
// add constant term if necessary
void makeFundamentalModeGreatAgain(BufferType& /*tractions*/,
std::vector<BufferType>& /*out*/,
influence::ElasticHelper<dim>& /*el*/) {}
// add constant term to a single layer
template <typename ST>
void makeFundamentalModeGreatAgain(StaticVector<Complex, ST, dim>& /*t*/,
out_t& /*u*/,
influence::ElasticHelper<dim>& /*el*/) {}
protected:
Accumulator<type, source_t> accumulator; ///< really only here for mesh
};
template <>
template <typename ST>
inline void
BoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
makeFundamentalModeGreatAgain(StaticVector<Complex, ST, dim>& t, out_t& u,
influence::ElasticHelper<dim>& el) {
if (mpi::rank() != 0)
return;
Matrix<Complex, 3, 3> correction({0, 0, 0, 0, 0, 0, -t(0) / el.mu,
-t(1) / el.mu,
-t(2) / (el.lambda + 2 * el.mu)});
u.symmetrize(correction);
}
template <>
inline void
BoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
makeFundamentalModeGreatAgain(BufferType& tractions,
std::vector<BufferType>& out,
influence::ElasticHelper<dim>& el) {
if (mpi::rank() != 0)
return;
VectorProxy<Complex, dim> t(&tractions(0));
for (auto&& grid : out) {
out_t u(&grid(0));
makeFundamentalModeGreatAgain(t, u, el);
}
}
} // namespace detail
} // namespace tamaas
#endif // BOUSSINESQ_HELPER_HH
/* -------------------------------------------------------------------------- */
c4science · Help