Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/influence.hh
d3681ba40060features/netcdf-mpi

Last Change
Raw File

influence.hh
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 "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
namespace {
const Complex I(0, 1); /* < imaginary unit */
constexpr Real sign(bool upper) { return (upper) ? 1. : -1.; }
}
template <bool conjugate, UInt dim_q> // dim + 1 to help dim template deduction
inline Vector<Complex, dim_q + 1> computeD(VectorProxy<const Real, dim_q>& q) {
constexpr Real s = -sign(conjugate);
const auto q_norm = q.l2norm();
return {{{s * I * q(0) / q_norm, s * I * q(1) / q_norm, 1}}};
}
template <UInt dim_q> // dim + 1 to help dim template deduction
inline Vector<Complex, dim_q + 1> computeD2(VectorProxy<const Real, dim_q>& q) {
const auto q_norm = q.l2norm();
return {{{I * q(0) / q_norm, I * q(1) / q_norm, 0}}};
}
/// Class for the Kelvin tensor
template <UInt dim, UInt derivative_order>
class Kelvin;
/// See details in notebook and Marc Bonnet's notes
template <>
class Kelvin<3, 0> {
protected:
constexpr static UInt dim = 3;
public:
Kelvin(Real mu, Real nu) : mu(mu), b(4 * (1 - nu)) {}
template <bool upper, typename ST>
inline Vector<Complex, dim> applyU0(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const;
template <bool upper, typename ST>
inline Vector<Complex, dim> applyU1(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const;
protected:
const Real mu, b;
};
/* -------------------------------------------------------------------------- */
/* Base Kelvin tensor */
/* -------------------------------------------------------------------------- */
template <bool upper, typename ST> //< this is ignored
inline Vector<Complex, Kelvin<3, 0>::dim>
Kelvin<3, 0>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const {
auto d2{computeD2(q)};
auto tmp{d2.dot(f) * d2};
tmp += b * f;
tmp(2) -= f(2);
tmp *= 1. / (2. * mu * q.l2norm() * b);
return tmp;
}
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 0>::dim>
Kelvin<3, 0>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const {
auto d{computeD<upper>(q)};
auto tmp{d.dot(f) * d};
tmp *= sign(upper) / (2. * mu * q.l2norm() * b);
return tmp;
}
/* -------------------------------------------------------------------------- */
/* Kelvin first derivative */
/* -------------------------------------------------------------------------- */
template <>
class Kelvin<3, 1> : protected Kelvin<3, 0> {
using parent = Kelvin<3, 0>;
protected:
static constexpr UInt dim = parent::dim;
public:
using parent::parent;
template <bool upper, typename ST>
inline Vector<Complex, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, typename ST>
inline Vector<Complex, dim>
applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
};
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 1>::dim>
Kelvin<3, 1>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp{f * computeD<upper>(q)};
auto res{Kelvin<3, 0>::applyU0<false>(q, tmp)};
res *= -sign(upper);
Vector<Real, dim> e3{{{0, 0, 1}}};
tmp.template mul<false>(f, e3);
res += Kelvin<3, 0>::applyU1<upper>(q, tmp);
res *= q.l2norm();
return res;
}
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 1>::dim>
Kelvin<3, 1>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp{f * computeD<upper>(q)};
auto res{Kelvin<3, 0>::applyU1<upper>(q, tmp)};
res *= -sign(upper) * q.l2norm();
return res;
}
/* -------------------------------------------------------------------------- */
/* Kelvin second derivative */
/* -------------------------------------------------------------------------- */
template <>
class Kelvin<3, 2> : protected Kelvin<3, 1> {
using parent = Kelvin<3, 1>;
static constexpr UInt dim = parent::dim;
public:
using parent::parent;
template <bool upper, typename ST>
inline Matrix<Complex, dim, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, typename ST>
inline Matrix<Complex, dim, dim>
applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyDiscontinuityTerm(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
};
template <bool upper, typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp{Kelvin<3, 1>::applyU0<upper>(q, f)};
tmp *= -sign(upper);
auto res{outer(computeD<upper>(q), tmp)};
Vector<Real, dim> e3{{{0, 0, 1}}};
res += outer(e3, Kelvin<3, 1>::applyU1<upper>(q, f));
res *= q.l2norm();
return res;
}
template <bool upper, typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp{Kelvin<3, 1>::applyU1<upper>(q, f)};
auto res{outer(computeD<upper>(q), tmp)};
res *= -sign(upper) * q.l2norm();
return res;
}
template <typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyDiscontinuityTerm(
VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
Vector<Real, dim> e3{{{0, 0, 1}}};
const auto d2{computeD2(q)};
auto tmp{f * e3};
auto res{Kelvin<3, 0>::applyU0<false>(q, tmp)};
res *= -2 * q.l2norm();
auto scal{tmp.dot(d2)};
auto last_component = tmp(2);
tmp = d2;
tmp *= scal;
tmp(2) += last_component;
tmp *= 1 / (this->mu * this->b);
res += tmp;
return outer(e3, res);
}
/* -------------------------------------------------------------------------- */
/* Implementation */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt order>
class KelvinIntegrator {
static inline Real int_exp(Real a, Real z) { return std::exp(a * z) / a; }
static inline Real int_exp_z(Real a, Real z) {
const Real a_inv = 1. / a;
return std::exp(a * z) * (z - a_inv) * a_inv;
}
static inline Real int_exp_zz(Real a, Real z) {
const Real a_inv = 1. / a;
return std::exp(a * z) * (z * z - 2 * a_inv * z + 2 * a_inv * a_inv) *
a_inv;
}
template <bool upper>
static inline Real g0(Real z) {
return std::exp(-sign(upper) * z);
}
template <bool upper>
static inline Real int_g0(const thrust::pair<Real, Real>& bounds,
Real alpha) {
alpha *= -sign(upper);
return int_exp(alpha, bounds.second) - int_exp(alpha, bounds.first);
}
template <bool upper>
static inline Real int_g0_z(const thrust::pair<Real, Real>& bounds,
Real alpha) {
alpha *= -sign(upper);
return int_exp_z(alpha, bounds.second) - int_exp_z(alpha, bounds.first);
}
template <bool upper>
static inline Real g1(Real z) {
return g0<upper>(z) * z;
}
template <bool upper>
static inline Real int_g1(const thrust::pair<Real, Real>& bounds,
Real alpha) {
return alpha * int_g0_z<upper>(bounds, alpha);
}
template <bool upper>
static inline Real int_g1_z(const thrust::pair<Real, Real>& bounds,
Real alpha) {
return alpha * (int_exp_zz(-alpha * sign(upper), bounds.second) -
int_exp_zz(-alpha * sign(upper), bounds.first));
}
public:
template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
typename KelvinType>
static inline void integrate(OutType&& out, InType&& in, KelvinType&& kelvin,
VectorProxy<const Real, dim_q>& q_, Real dij,
Real dl);
};
// Integrate against a constant shape function
template <>
template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
typename KelvinType>
inline void KelvinIntegrator<0>::integrate(OutType&& out, InType&& in,
KelvinType&& kelvin,
VectorProxy<const Real, dim_q>& q_,
Real dij, Real dl) {
const auto q = q_.l2norm();
const auto alpha = q * dl, beta = q * dij;
// Integrate on [-1, 0]
const auto bounds_low = thrust::make_pair(-0.5, 0);
constexpr bool domain_low = yj_xi > 0;
auto g0_mul{kelvin.template applyU0<domain_low>(q_, in)};
auto g1_mul{kelvin.template applyU1<domain_low>(q_, in)};
decltype(g0_mul) res;
res = g0_mul * (g0<domain_low>(beta) * int_g0<domain_low>(bounds_low, alpha));
res +=
g1_mul * (g1<domain_low>(beta) * int_g0<domain_low>(bounds_low, alpha) +
g0<domain_low>(beta) * int_g1<domain_low>(bounds_low, alpha));
// Integrate on [0, 1]
const auto bounds_up = thrust::make_pair(0, 0.5);
constexpr bool domain_up = yj_xi >= 0;
// Avoid computing
if (yj_xi == 0) {
g0_mul = kelvin.template applyU0<true>(q_, in);
g1_mul = kelvin.template applyU1<true>(q_, in);
}
res += g0_mul * (g0<domain_up>(beta) * int_g0<domain_up>(bounds_up, alpha));
res += g1_mul * (g1<domain_up>(beta) * int_g0<domain_up>(bounds_up, alpha) +
g0<domain_up>(beta) * int_g1<domain_up>(bounds_up, alpha));
res *= dl;
out += res;
}
// Integrate against a linear shape function
template <>
template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
typename KelvinType>
inline void KelvinIntegrator<1>::integrate(OutType&& out, InType&& in,
KelvinType&& kelvin,
VectorProxy<const Real, dim_q>& q_,
Real dij, Real dl) {
const auto q = q_.l2norm();
const auto alpha = q * dl, beta = q * dij;
// Integrate on [-1, 0]
const auto bounds_low = thrust::make_pair(-1, 0);
constexpr bool domain_low = yj_xi > 0;
auto g0_mul{kelvin.template applyU0<domain_low>(q_, in)};
auto g1_mul{kelvin.template applyU1<domain_low>(q_, in)};
decltype(g0_mul) res;
res = g0_mul *
(g0<domain_low>(beta) * (int_g0<domain_low>(bounds_low, alpha) +
int_g0_z<domain_low>(bounds_low, alpha)));
res += g1_mul *
(g1<domain_low>(beta) * (int_g0<domain_low>(bounds_low, alpha) +
int_g0_z<domain_low>(bounds_low, alpha)) +
g0<domain_low>(beta) * (int_g1<domain_low>(bounds_low, alpha) +
int_g1_z<domain_low>(bounds_low, alpha)));
// Integrate on [0, 1]
const auto bounds_up = thrust::make_pair(0, 1);
constexpr bool domain_up = yj_xi >= 0;
// Avoid computing
if (yj_xi == 0) {
g0_mul = kelvin.template applyU0<true>(q_, in);
g1_mul = kelvin.template applyU1<true>(q_, in);
}
res +=
g0_mul * (g0<domain_up>(beta) * (int_g0<domain_up>(bounds_up, alpha) -
int_g0_z<domain_up>(bounds_up, alpha)));
res +=
g1_mul * (g1<domain_up>(beta) * (int_g0<domain_up>(bounds_up, alpha) -
int_g0_z<domain_up>(bounds_up, alpha)) +
g0<domain_up>(beta) * (int_g1<domain_up>(bounds_up, alpha) -
int_g1_z<domain_up>(bounds_up, alpha)));
res *= dl;
out += res;
}
/* -------------------------------------------------------------------------- */
/// Class for the Boussinesq tensor
template <UInt dim, UInt derivative_order>
class Boussinesq;
template <>
class Boussinesq<3, 0> {
protected:
static constexpr UInt dim = 3;
public:
/// Constructor
Boussinesq(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
/// Apply operator
inline void apply(VectorProxy<Complex, dim>& u,
MatrixProxy<Complex, dim, dim>& eps,
VectorProxy<const Real, dim - 1>& q, Real y);
Matrix<Complex, dim, dim> computeSigma(MatrixProxy<Complex, dim, dim>& gradu) {
auto trace{gradu.trace()};
Matrix<Complex, dim, dim> sigma;
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
sigma(i, j) =
(i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
return sigma;
}
protected:
const Real mu, nu, lambda;
};
inline void Boussinesq<3, 0>::apply(VectorProxy<Complex, dim>& u,
MatrixProxy<Complex, dim, dim>& gradu,
VectorProxy<const Real, dim - 1>& q,
Real y) {
Vector<Real, dim> e3{{{0, 0, -1}}}; // WHY???!!!
auto q_norm{q.l2norm()};
auto traction{computeSigma(gradu) * e3};
auto d{computeD<false>(q)}, dbar{computeD<true>(q)}, d2{computeD2(q)};
auto scal{dbar.dot(traction)};
auto res{(scal * q_norm * y) * dbar};
res += 2. * traction;
res += ((1 - 2 * nu) * scal) * d;
res += d2.dot(traction) * d2;
res(2) -= traction(2);
u += (std::exp(-q_norm * y) / (2*mu*q_norm)) * res;
}
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help