Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/model/influence.hh
c74466647365maxwell_viscoelastic…

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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INFLUENCE_HH
#define INFLUENCE_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
#include <type_traits>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
namespace {
const Complex I(0, 1); /// < imaginary unit
constexpr inline 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 {{{Complex{0, s * q(0) / q_norm}, Complex{0, s * q(1) / q_norm},
Complex{1, 0}}}};
}
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 {
{{Complex{0, q(0) / q_norm}, Complex{0, q(1) / q_norm}, Complex{0, 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;
constexpr static UInt order = 0;
public:
Kelvin(Real mu, Real nu) : mu(mu), b(4 * (1 - nu)) {}
template <bool upper, bool apply_q_power = false, typename ST>
inline Vector<Complex, dim> applyU0(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const;
template <bool upper, bool apply_q_power = false, 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, bool apply_q_power, typename ST> //< upper 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 q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
d2 *= d2.dot(f);
d2 += b * f;
d2(2) -= f(2);
d2 *= q_power / (2. * mu * b);
return d2;
}
template <bool upper, bool apply_q_power, 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 q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
d *= d.dot(f);
d *= q_power * sign(upper) / (2. * mu * b);
return d;
}
/* -------------------------------------------------------------------------- */
/* Kelvin first derivative */
/* -------------------------------------------------------------------------- */
template <>
class Kelvin<3, 1> : protected Kelvin<3, 0> {
using parent = Kelvin<3, 0>;
protected:
static constexpr UInt dim = parent::dim;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template <bool upper, bool apply_q_power = false,
typename ST> // < apply_q ignored
inline Vector<Complex, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, bool apply_q_power = false,
typename ST> // < apply_q ignored
inline Vector<Complex, dim>
applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
};
template <bool upper, bool apply_q_power, 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);
return res;
}
template <bool upper, bool apply_q_power, 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);
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;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template <bool upper, bool apply_q_power = false, typename ST>
inline Matrix<Complex, dim, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, bool apply_q_power = false, 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, bool apply_q_power, 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);
auto q_power = (apply_q_power) ? q.l2norm() : 1.;
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 *= -1. * q_power; // << -1 for grad_x
return res;
}
template <bool upper, bool apply_q_power, 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 q_power = (apply_q_power) ? q.l2norm() : 1.;
auto res = outer(computeD<upper>(q), tmp);
res *= -1. * -sign(upper) * q_power; // << -1 for grad_x
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}}};
auto res = (f * e3);
res *= -this->b;
res(2) += 2. * f(2, 2);
e3(2) *= -1 / (this->mu * this->b); // << -1 for grad_x
return outer(e3, res);
}
/* -------------------------------------------------------------------------- */
/* Implementation */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt order>
struct 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));
}
template <Int yj_xi, UInt shape_num, 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);
};
namespace detail {
template <UInt shape_num>
struct LinearShape {
using integrator = KelvinIntegrator<1>;
template <bool domain, typename T>
static inline T integrateG0(T& g0_mul, Real alpha, Real beta);
template <bool domain, typename T>
static inline T integrateG1(T& g1_mul, Real alpha, Real beta);
};
template <UInt shape_num>
template <bool domain, typename T>
inline T LinearShape<shape_num>::integrateG0(T& g0_mul, Real alpha,
Real beta) {
constexpr Int mul = (shape_num == 0) ? -1 : 1;
const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
return g0_mul * (integrator::g0<domain>(beta) *
(integrator::int_g0<domain>(bounds, alpha) +
mul * integrator::int_g0_z<domain>(bounds, alpha)));
}
template <UInt shape_num>
template <bool domain, typename T>
inline T LinearShape<shape_num>::integrateG1(T& g1_mul, Real alpha,
Real beta) {
constexpr Int mul = (shape_num == 0) ? -1 : 1;
const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
return g1_mul * (integrator::g1<domain>(beta) *
(integrator::int_g0<domain>(bounds, alpha) +
mul * integrator::int_g0_z<domain>(bounds, alpha)) +
integrator::g0<domain>(beta) *
(integrator::int_g1<domain>(bounds, alpha) +
mul * integrator::int_g1_z<domain>(bounds, alpha)));
}
}
// Integrate against a linear shape function
template <>
template <Int yj_xi, UInt shape_num, 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;
constexpr bool domain = (shape_num == 0) ? yj_xi >= 0 : yj_xi > 0;
// Applying kelvin with the power of q multiplier
auto g0_mul = kelvin.template applyU0<domain, true>(q_, in);
auto g1_mul = kelvin.template applyU1<domain, true>(q_, in);
decltype(g0_mul) res;
res = detail::LinearShape<shape_num>::template integrateG0<domain>(
g0_mul, alpha, beta);
res += detail::LinearShape<shape_num>::template integrateG1<domain>(
g1_mul, alpha, beta);
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;
static constexpr UInt order = 0;
public:
/// Constructor
Boussinesq(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
template <bool apply_q_power = false, typename ST>
inline Vector<Complex, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <bool apply_q_power = false, typename ST>
inline Vector<Complex, dim>
applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
protected:
const Real mu, nu, lambda;
};
template <bool apply_q_power, typename ST>
inline Vector<Complex, Boussinesq<3, 0>::dim>
Boussinesq<3, 0>::applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// This guy has the property Bt(q) = B(-q). Since we know we recieve -q as an
// argument, we can just apply B without transpose
auto dplus = computeD<true>(q), dminus = computeD<false>(q), d2 = computeD2(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
auto res = 2. * t;
res += ((1. - 2. * nu) * dminus.dot(t)) * dplus;
res += d2.dot(t) * d2;
res(2) -= t(2);
res *= q_power / (2. * mu);
return res;
}
template <bool apply_q_power, typename ST>
inline Vector<Complex, Boussinesq<3, 0>::dim>
Boussinesq<3, 0>::applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// This guy does not have the property that Bt(q) = B(-q), so we have to apply
// Bt(q) for real. This makes us correct (in the line below) for the fact that
// we recieve -q as an argument.
auto dplus = computeD<false>(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
auto res = q_power * (dplus.dot(t) / (2. * mu)) * dplus;
return res;
}
/* -------------------------------------------------------------------------- */
template <>
class Boussinesq<3, 1> : protected Boussinesq<3, 0> {
protected:
using parent = Boussinesq<3, 0>;
static constexpr UInt dim = parent::dim;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template <bool apply_q_power = false, typename ST>
inline Matrix<Complex, dim, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <bool apply_q_power = false, typename ST>
inline Matrix<Complex, dim, dim>
applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
};
template <bool apply_q_power, typename ST>
inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
Boussinesq<3, 1>::applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// taking into account K(-q) convolution
auto dbar = computeD<false>(q);
Vector<Real, dim> e3{{{0, 0, 1}}};
auto res = outer(dbar, Boussinesq<3, 0>::applyU0(t, q));
res *= -1.;
res += outer(e3, Boussinesq<3, 0>::applyU1(t, q));
return res;
}
template <bool apply_q_power, typename ST>
inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
Boussinesq<3, 1>::applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// taking into account K(-q) convolution
auto dbar = computeD<false>(q);
auto res = outer(dbar, Boussinesq<3, 0>::applyU1(t, q));
res *= -1;
return res;
}
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif // INFLUENCE_HH
c4science · Help