Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65751902
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
Wed, Jun 5, 23:44
Size
16 KB
Mime Type
text/x-c++
Expires
Fri, Jun 7, 23:44 (2 d)
Engine
blob
Format
Raw Data
Handle
18121797
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"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
namespace {
const Complex I(0, 1); /// < imaginary unit
constexpr inline Real sign(bool upper) { return (upper) ? 1. : -1.; }
}
template <UInt dim>
struct ElasticHelper {
ElasticHelper(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
template <typename DT, typename ST>
Matrix<DT, dim, dim>
operator()(const StaticMatrix<DT, ST, dim, dim>& gradu) const {
auto trace{gradu.trace()};
Matrix<DT, 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;
}
const Real mu, nu, lambda;
};
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 *= -1. * q.l2norm(); // << -1 for grad_x
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 *= -1. * -sign(upper) * q.l2norm(); // << -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}}};
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;
e3(2) *= -1; // << -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 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)) {}
template <typename ST>
inline Vector<Complex, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <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 <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 {
// anticipating K(-q) convolution (no change for d2 cause d2xd2 not sensitive)
auto d{computeD<true>(q)}, dbar{computeD<false>(q)}, d2{computeD2(q)};
auto res{2. * t};
res += ((1. - 2. * nu) * dbar.dot(t)) * d;
res += d2.dot(t) * d2;
res(2) -= t(2);
res *= (1. / (2. * mu * q.l2norm()));
return res;
}
template <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 {
// anticipating K(-q) convolution (no change for d2 cause d2xd2 not sensitive)
auto dbar{computeD<false>(q)};
auto res{(dbar.dot(t) / (2. * mu * q.l2norm())) * dbar};
return res;
}
/* -------------------------------------------------------------------------- */
template <>
class Boussinesq<3, 1> : protected Boussinesq<3, 0> {
protected:
using parent = Boussinesq<3, 0>;
static const UInt dim = parent::dim;
public:
using parent::parent;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
};
template <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));
res *= q.l2norm();
return res;
}
template <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 *= -q.l2norm();
return res;
}
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif // INFLUENCE_HH
Event Timeline
Log In to Comment