Homec4science
Diffusion tamaas (ci-fix-static-checks) Browse

tamaas/src/model/westergaard.cpp
04c23013552dci-fix-static-checks

Last Change
Raw File

westergaard.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* 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 "grid_hermitian.hh"
#include "grid_view.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "model_type.hh"
#include "static_types.hh"
#include "westergaard.hh"
#include <functional>
#include <numeric>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
Westergaard<mtype, otype>::Westergaard(Model* model)
: IntegralOperator(model), buffer(), engine(FFTEngine::makeEngine()) {
// Copy sizes
auto bdisc = model->getBoundaryDiscretization();
auto boundary_hermitian_sizes =
GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
bdisc);
constexpr UInt nb_components = trait::components;
buffer.setNbComponents(nb_components);
buffer.resize(boundary_hermitian_sizes);
influence = allocateGrid<mtype, true, Real, GridHermitian>(
boundary_hermitian_sizes, nb_components * nb_components);
(*this)["influence"] = influence;
initInfluence();
}
/* -------------------------------------------------------------------------- */
namespace detail {
template <UInt m, UInt j>
struct boundary_fft_helper {
template <typename Buffer, typename Out>
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
auto view = make_view(out, 0);
e.backward(view, buffer);
}
};
template <UInt m>
struct boundary_fft_helper<m, m> {
template <typename Buffer, typename Out>
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
e.backward(out, buffer);
}
};
} // namespace detail
template <model_type mtype, IntegralOperator::kind otype>
template <typename Functor>
void Westergaard<mtype, otype>::fourierApply(Functor func, GridBase<Real>& in,
GridBase<Real>& out) const try {
auto& i = dynamic_cast<Grid<Real, bdim>&>(in);
auto& full_out = dynamic_cast<Grid<Real, dim>&>(out);
engine->forward(i, buffer);
/// Applying influence
func(buffer, *influence);
/// Backward fourier transform on boundary only
detail::boundary_fft_helper<bdim, dim>::backwardTransform(*engine, buffer,
full_out);
} catch (const std::bad_cast& e) {
throw model_type_error{
TAMAAS_MSG("Neumann and dirichlet types do not match model type")};
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
template <typename Functor>
void Westergaard<mtype, otype>::initFromFunctor(Functor func) {
// Compute wavevectors for influence
auto wavevectors = FFTEngine::template computeFrequencies<Real, bdim, true>(
influence->sizes());
// Get boundary physical size
auto system_size = this->model->getBoundarySystemSize();
VectorProxy<const Real, bdim> domain(system_size[0]);
// Normalize wavevectors
wavevectors *= 2 * M_PI;
wavevectors /= domain;
// Apply functor
Loop::loop(func, range<VectorProxy<Real, bdim>>(wavevectors),
range<MatrixProxy<Complex, comp, comp>>(*influence));
if (mpi::rank() == 0) {
// Set fundamental mode to zero
MatrixProxy<Complex, comp, comp> mat((*influence)(0));
mat = 0;
}
}
/* -------------------------------------------------------------------------- */
/// \cond DO_NOT_DOCUMENT
#define NEUMANN_BASIC(type) \
template <> \
void Westergaard<type, IntegralOperator::neumann>::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> q, \
MatrixProxy<Complex, comp, comp> k) { \
k(0, 0) = 2. / (E_star * q.l2norm()); \
}; \
initFromFunctor(basic); \
}
#define DIRICHLET_BASIC(type) \
template <> \
void Westergaard<type, IntegralOperator::dirichlet>::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> q, \
MatrixProxy<Complex, comp, comp> k) { \
k(0, 0) = E_star * q.l2norm() / 2; \
}; \
initFromFunctor(basic); \
}
NEUMANN_BASIC(model_type::basic_1d)
NEUMANN_BASIC(model_type::basic_2d)
DIRICHLET_BASIC(model_type::basic_1d)
DIRICHLET_BASIC(model_type::basic_2d)
#undef NEUMANN_BASIC
#undef DIRICHLET_BASIC
/* -------------------------------------------------------------------------- */
template <>
void Westergaard<model_type::surface_1d,
IntegralOperator::neumann>::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> q_,
MatrixProxy<Complex, comp, comp> F) {
Real q = q_(0);
q *= 1. / q_.l2norm();
F(0, 0) = 2 * (1 + nu) * (1 - nu * q * q);
F(1, 1) = 2 * (1 - nu * nu);
F(0, 1) = I * q * (1 + nu) * (1 - 2 * nu);
F(1, 0) = -F(0, 1);
F *= 1. / (E * q_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <>
void Westergaard<model_type::surface_2d,
IntegralOperator::neumann>::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> q_,
MatrixProxy<Complex, comp, comp> F) {
Vector<Real, bdim> q(q_);
q *= 1 / q_.l2norm();
F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0));
F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1));
F(2, 2) = 2 * (1 - nu * nu);
F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu);
F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu);
F(1, 2) = I * q(1) * (1 + nu) * (1. - 2. * nu);
F(2, 0) = -F(0, 2);
F(2, 1) = -F(1, 2);
F *= 1. / (E * q_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
void Westergaard<mtype, otype>::initInfluence() {}
/// \endcond
/* ------------------------------------------------------------------------ */
template <model_type mtype, IntegralOperator::kind otype>
void Westergaard<mtype, otype>::apply(GridBase<Real>& input,
GridBase<Real>& output) const {
auto apply = [](decltype(buffer)& buffer,
const decltype(*influence)& influence) {
Loop::loop(
[] CUDA_LAMBDA(VectorProxy<Complex, comp> i,
MatrixProxy<const Complex, comp, comp> F) { i = F * i; },
range<VectorProxy<Complex, comp>>(buffer),
range<MatrixProxy<const Complex, comp, comp>>(influence));
};
fourierApply(apply, input, output);
}
/* -------------------------------------------------------------------------- */
template <model_type type, IntegralOperator::kind otype>
Real Westergaard<type, otype>::getOperatorNorm() {
constexpr UInt comp = model_type_traits<type>::components;
Real _norm = 0.0;
_norm = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(MatrixProxy<Complex, comp, comp> m) {
Real n = thrust::norm(m.l2squared());
return std::isnan(n) ? 0 : n;
},
range<MatrixProxy<Complex, comp, comp>>(*influence));
auto size = model->getSystemSize();
auto disc = model->getDiscretization();
auto dim = model_type_traits<type>::dimension;
/// TODO: why?
switch (dim) {
case 1:
_norm /= (size[0] / disc[0]) * (size[0] / disc[0]);
break;
default:
for (UInt i = 0; i < dim; i++) {
_norm /= size[i] / disc[i];
}
}
return std::sqrt(_norm);
}
/* -------------------------------------------------------------------------- */
template <model_type type, IntegralOperator::kind otype>
std::pair<UInt, UInt> Westergaard<type, otype>::matvecShape() const {
auto n = model->getBoundaryDiscretization();
auto N = std::accumulate(n.begin(), n.end(), comp, std::multiplies<>());
return std::make_pair(N, N);
}
/* -------------------------------------------------------------------------- */
template <model_type type, IntegralOperator::kind otype>
GridBase<Real> Westergaard<type, otype>::matvec(GridBase<Real>& X) const {
auto n = model->getBoundaryDiscretization();
// Creating correct input
Grid<Real, bdim> x{n, comp, X.view()};
// Creating correct output
if (bdim != dim)
n.insert(n.begin(), 1);
Grid<Real, dim> y(n, comp);
apply(x, y);
return y;
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Westergaard<model_type::basic_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::basic_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::basic_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::basic_2d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::surface_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::surface_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::surface_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::surface_2d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::volume_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::volume_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::volume_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::volume_2d, IntegralOperator::dirichlet>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
/* -------------------------------------------------------------------------- */
c4science · Help