Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/westergaard.cpp
966dd0718da2features/netcdf-mpi

Last Change
Raw File

westergaard.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 "fft_plan_manager.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "grid_view.hh"
#include "static_types.hh"
#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::type otype>
Westergaard<mtype, otype>::Westergaard(Model* model) : IntegralOperator(model) {
// Copy sizes
std::array<UInt, trait::dimension> sizes;
std::copy(model->getDiscretization().begin(),
model->getDiscretization().end(), sizes.begin());
// Copy boundary sizes
auto sizes_it = sizes.begin();
std::array<UInt, trait::boundary_dimension> boundary_hermitian_sizes;
// Ignoring first dimension if model is volumetric
if (trait::dimension > trait::boundary_dimension)
++sizes_it;
std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin());
boundary_hermitian_sizes =
GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
boundary_hermitian_sizes);
constexpr UInt nb_components = trait::components;
buffer.setNbComponents(nb_components);
buffer.resize(boundary_hermitian_sizes);
influence.setNbComponents(nb_components * nb_components);
influence.resize(boundary_hermitian_sizes);
initInfluence();
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::type otype>
template <typename Functor>
void Westergaard<mtype, otype>::fourierApply(Functor func, GridBase<Real>& in,
GridBase<Real>& out) const try {
Grid<Real, bdim>& i = dynamic_cast<decltype(i)>(in);
Grid<Real, dim>& full_out = dynamic_cast<decltype(full_out)>(out);
FFTPlanManager::get().createPlan(i, buffer).forwardTransform();
/// Applying influence
func(buffer, influence);
/// Backward fourier transform on boundary only
if constexpr (bdim == dim)
FFTPlanManager::get().createPlan(full_out, buffer).backwardTransform();
else {
auto view = make_view(full_out, 0);
FFTPlanManager::get()
.createPlan(view, buffer)
.backwardTransform();
}
} catch (const std::bad_cast& e) {
TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::type otype>
template <typename Functor>
void Westergaard<mtype, otype>::initFromFunctor(Functor func) {
// Compute wavevectors for influence
auto wavevectors = FFTransform<Real, bdim>::template computeFrequencies<true>(
influence.sizes());
// Get boundary physical size
VectorProxy<const Real, bdim> domain(
this->model->getSystemSize()[(bdim == dim) ? 0 : 1]);
// Normalize wavevectors
wavevectors *= 2 * M_PI;
wavevectors /= domain;
// Apply functor
Loop::stridedLoop(func, wavevectors, influence);
// Set fundamental mode to zero
MatrixProxy<Complex, trait::components, trait::components> mat(influence(0));
mat = 0;
}
/* -------------------------------------------------------------------------- */
#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, \
Complex & k) { k = 2. / 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, \
Complex & k) { k = 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_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) {
Real q_norm = q.l2norm();
q /= q_norm;
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_norm);
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::type otype>
void Westergaard<mtype, otype>::initInfluence() {
TAMAAS_EXCEPTION("the requested operator has not been implemented");
}
/* ------------------------------------------------------------------------ */
template <model_type mtype, IntegralOperator::type otype>
void Westergaard<mtype, otype>::apply(GridBase<Real>& input,
GridBase<Real>& output) const {
auto apply = [](decltype(buffer)& buffer,
const decltype(influence)& influence) {
Loop::stridedLoop([] CUDA_LAMBDA(VectorProxy<Complex, comp> && i,
MatrixProxy<const Complex, comp, comp> &&
F) { i = F * i; },
buffer, influence);
};
fourierApply(apply, input, output);
}
/* -------------------------------------------------------------------------- */
/* 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>;
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help