Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/model/influence.hh
dc7ef4e3f247features/declarative

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.; }
} // namespace
template <bool conjugate, UInt dim_q> // dim + 1 to help dim template deduction
inline Vector<Complex, dim_q + 1>
computeD(const 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(const VectorProxy<const Real, dim_q>& q) {
const auto q_norm = (q.l2norm() < 1e-16) ? 1 : 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;
/**
* @brief Kelvin base tensor
* See arXiv:1811.11558 for details.
*/
template <>
class Kelvin<3, 0> {
protected:
static constexpr UInt dim = 3;
static constexpr 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(const VectorProxy<const Real, dim - 1>& q,
const StaticVector<Complex, ST, dim>& f) const {
// upper is ignored here
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 = false, typename ST>
inline Vector<Complex, dim> applyU1(const VectorProxy<const Real, dim - 1>& q,
const 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;
}
protected:
const Real mu, b;
};
/* -------------------------------------------------------------------------- */
/// 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>
inline Vector<Complex, dim>
applyU0(const VectorProxy<const Real, dim - 1>& q,
const StaticMatrix<Complex, ST, dim, dim>& f) const {
// apply_q_power is ignored here
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 = false, typename ST>
inline Vector<Complex, dim>
applyU1(const VectorProxy<const Real, dim - 1>& q,
const StaticMatrix<Complex, ST, dim, dim>& f) const {
// apply_q_power is ignored here
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(const VectorProxy<const Real, dim - 1>& q,
const 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 = false, typename ST>
inline Matrix<Complex, dim, dim>
applyU1(const VectorProxy<const Real, dim - 1>& q,
const 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, dim, dim>
applyDiscontinuityTerm(const 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);
}
};
/* -------------------------------------------------------------------------- */
/// 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(const StaticVector<Complex, ST, dim>& t,
const 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 = false, typename ST>
inline Vector<Complex, dim>
applyU1(const StaticVector<Complex, ST, dim>& t,
const 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;
}
protected:
const Real mu, nu, lambda;
};
/* -------------------------------------------------------------------------- */
/// Boussinesq first gradient
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(const StaticVector<Complex, ST, dim>& t,
const 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 = false, typename ST>
inline Matrix<Complex, dim, dim>
applyU1(const StaticVector<Complex, ST, dim>& t,
const 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