Page MenuHomec4science
Files F60720068

influence.hh
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 04:36

influence.hh
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INFLUENCE_HH
#define INFLUENCE_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
#include <type_traits>
/* -------------------------------------------------------------------------- */
namespace 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>
__device__ __host__ 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
__device__ __host__ 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}}}};
}
/* -------------------------------------------------------------------------- */
/// Functor to apply Hooke's tensor
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>
__device__ __host__ Matrix<std::remove_cv_t<DT>, dim, dim>
operator()(const StaticMatrix<DT, ST, dim, dim>& gradu) const {
auto trace = gradu.trace();
Matrix<std::remove_cv_t<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;
}
template <typename DT, typename ST>
__device__ __host__ SymMatrix<std::remove_cv_t<DT>, dim>
operator()(const StaticSymMatrix<DT, ST, dim>& eps) const {
SymMatrix<std::remove_cv_t<DT>, dim> sigma;
auto trace = eps.trace();
for (UInt i = 0; i < dim; ++i)
sigma(i) = lambda * trace + 2 * mu * eps(i);
for (UInt i = dim; i < voigt_size<dim>::value; ++i)
sigma(i) = 2 * mu * eps(i);
return sigma;
}
template <typename DT, typename ST>
__device__ __host__ Matrix<std::remove_cv_t<DT>, dim, dim>
inverse(const StaticMatrix<DT, ST, dim, dim>& sigma) const {
Matrix<std::remove_cv_t<DT>, dim, dim> epsilon;
auto trace = sigma.trace();
auto c1 = 1. / (2 * mu);
auto c2 = -lambda / (3 * lambda + 2 * mu);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) = c1 * (sigma(i, j) + c2 * trace * (i == j));
}
const Real mu, nu, lambda;
};
/* -------------------------------------------------------------------------- */
/// 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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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> __device__ __host__
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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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>
__device__ __host__ 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
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // INFLUENCE_HH

Event Timeline

c4science · Help