Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/westergaard.cpp
383b32a7e7befeatures/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 "grid_view.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "static_types.hh"
#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
Westergaard<mtype, otype>::Westergaard(Model* model) : IntegralOperator(model) {
// Copy sizes
std::array<UInt, trait::boundary_dimension> sizes;
auto bdisc = model->getBoundaryDiscretization();
// Ignoring first dimension if model is volumetric
std::copy(bdisc.begin(), bdisc.end(), sizes.begin());
auto boundary_hermitian_sizes =
GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
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();
}
/* -------------------------------------------------------------------------- */
namespace detail {
template <UInt m, UInt j>
struct boundary_fft_helper {
template <typename Buffer, typename Out>
static void backwardTransform(Buffer&& buffer, Out&& out) {
auto view = make_view(out, 0);
FFTPlanManager::get().createPlan(view, buffer).backwardTransform();
}
};
template <UInt m>
struct boundary_fft_helper<m, m> {
template <typename Buffer, typename Out>
static void backwardTransform(Buffer&& buffer, Out&& out) {
FFTPlanManager::get().createPlan(out, buffer).backwardTransform();
}
};
}
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 {
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
detail::boundary_fft_helper<bdim, dim>::backwardTransform(buffer, full_out);
} catch (const std::bad_cast& e) {
TAMAAS_EXCEPTION("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 = FFTransform<Real, bdim>::template computeFrequencies<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));
// Set fundamental mode to zero
MatrixProxy<Complex, comp, comp> 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, \
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() {}
/* ------------------------------------------------------------------------ */
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 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