Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65879888
influence.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 6, 20:17
Size
16 KB
Mime Type
text/x-c++
Expires
Sat, Jun 8, 20:17 (2 d)
Engine
blob
Format
Raw Data
Handle
18140641
Attached To
rTAMAAS tamaas
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
Event Timeline
Log In to Comment