Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/influence.hh
56e1e9fb1c3cfeatures/netcdf-mpi

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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
/// Mother class for influence functions
template <UInt dim>
class Influence {
public:
Influence(StaticVector<const Real, dim> domain) : domain(domain) {}
protected:
StaticVector<const Real, dim> domain;
};
/// Class representing the basic influence functions from Stanley & Kato (1997)
template <UInt dim>
class Basic : public Influence<dim> {
public:
Basic(Real E_hertz, StaticVector<const Real, dim> domain)
: Influence<dim>(domain), E_hertz(E_hertz) {}
/// Influence function: computes the influence coefficient F
CUDA_LAMBDA void operator()(StaticVector<Real, dim>&& q_vec,
Complex& F) const {
q_vec *= 2 * M_PI;
q_vec /= this->domain;
F = 2. / (q_vec.l2norm() * E_hertz);
}
private:
Real E_hertz;
};
template <UInt dim>
class Surface : public Influence<dim> {};
/// Class for the influence with all components of applied tractions (dim = 2)
template <>
class Surface<2> : public Influence<2> {
constexpr static UInt dim = 2;
public:
Surface(Real E, Real nu, StaticVector<const Real, dim> domain)
: Influence<dim>(domain), E(E), nu(nu) {}
CUDA_LAMBDA void
operator()(StaticVector<Real, dim>&& q_,
StaticMatrix<Complex, dim + 1, dim + 1>&& F) const {
q_ *= 2 * M_PI;
q_ /= this->domain;
const Real q = q_.l2norm();
// First line of influence matrix
F(0, 0) =
1. / (q * q) * (-nu * nu * q_(0) * q_(0) - nu * q_(1) * q_(1) + q * q);
F(0, 1) = -q_(0) * q_(1) / (q * q) * (1 + nu) * 2 * nu;
F(0, 2) = I * q_(0) / q * (1 + nu) * (1 - 2 * nu);
// Second line of influence matrix
F(1, 0) = -q_(0) * q_(1) / (q * q) * (1 + nu) * 2 * nu;
F(1, 1) =
1. / (q * q) * (-nu * nu * q_(1) * q_(1) - nu * q_(0) * q_(0) + q * q);
F(1, 2) = I * q_(1) / q * (1 + nu) * (1 - 2 * nu);
// Third line of influence matrix
F(2, 0) = -F(0, 2);
F(2, 1) = -F(1, 2);
F(2, 2) = 2 * (1 - nu * nu);
F *= 1. / (E * q);
}
private:
Real E, nu;
const Complex I = Complex(0, 1);
};
/// Class for the influence with all components of applied tractions (dim = 2)
template <>
class Surface<1> : public Influence<1> {
constexpr static UInt dim = 1;
public:
Surface(Real E, Real nu, StaticVector<const Real, dim> domain)
: Influence<dim>(domain), E(E), nu(nu) {}
CUDA_LAMBDA void
operator()(StaticVector<Real, dim>&& q_,
StaticMatrix<Complex, dim + 1, dim + 1>&& F) const {
q_ *= 2 * M_PI;
q_ /= this->domain;
const Real q = q_.l2norm();
// First line of influence matrix
F(0, 0) = 2. / q * (1 - nu * nu);
F(0, 1) = I / q_(0) * (1 + nu) * (1 - 2 * nu);
// Second line of influence matrix
F(1, 0) = -F(0, 1);
F(1, 1) = F(0, 0);
F *= 1. / E;
}
private:
Real E, nu;
const Complex I = Complex(0, 1);
};
/// Class for the Kelvin tensor
template <UInt dim, UInt derivative_order>
class Kelvin;
template <UInt derivative_order>
class Kelvin<3, derivative_order> : Influence<3> {
constexpr static UInt dim = 3;
public:
Kelvin(Real mu, Real nu, StaticVector<const Real, dim> domain)
: Influence<dim>(domain), mu(mu), nu(nu) {}
private:
/// Compute leading linear term A
template <bool upper>
inline void computeA(StaticVector<const Real, dim - 1>& q_,
StaticMatrix<Complex, dim, dim>& A) const;
/// Compute leading constant term B, which does not depend on y_3
inline void computeB(StaticVector<const Real, dim - 1>& q_,
StaticMatrix<Complex, dim, dim>& B) const;
public:
/// Compute leading terms
template <bool upper>
inline void computeTerms(StaticVector<const Real, dim - 1>& q_,
StaticMatrix<Complex, dim, dim>& A,
StaticMatrix<Complex, dim, dim>& B) const {
computeA<upper>(q_, A);
computeB(q_, B);
}
private:
Real mu, nu;
const Complex I = Complex(0, 1);
};
/* -------------------------------------------------------------------------- */
/* Base Kelvin tensor */
/* -------------------------------------------------------------------------- */
template <>
inline void Kelvin<3, 0>::computeB(StaticVector<const Real, dim - 1>& q_,
StaticMatrix<Complex, dim, dim>& B) const {
const Real q = q_.l2norm();
const Real bq1 = q_(0) / q;
const Real bq2 = q_(1) / q;
for (UInt i = 0; i < dim; ++i)
B(i, i) = 4 * (1 - nu);
B(2, 2) -= 1;
B(0, 0) -= bq1 * bq1;
B(1, 1) -= bq2 * bq2;
B(0, 1) -= bq1 * bq2;
B(1, 0) -= bq1 * bq2;
B(0, 2) = B(1, 2) = B(2, 0) = B(2, 1) = 0;
B *= (1 / (8 * mu * q * (1 - nu)));
}
template <>
template <bool upper>
inline void Kelvin<3, 0>::computeA(StaticVector<const Real, dim - 1>& q_,
StaticMatrix<Complex, dim, dim>& A) const {
const Real q = q_.l2norm();
const Real bq1 = q_(0) / q;
const Real bq2 = q_(1) / q;
constexpr Real sign = (upper) ? 1. : -1;
A(0, 0) = sign * bq1 * bq1;
A(1, 1) = sign * bq2 * bq2;
A(2, 2) = 1;
A(1, 0) = A(0, 1) = sign * bq1 * bq2;
A(2, 0) = A(0, 2) = -I * bq1;
A(2, 1) = A(1, 2) = -I * bq2;
A *= (1 / (8 * mu * (1 - nu)));
}
/* -------------------------------------------------------------------------- */
class KelvinIntegrator {
static inline Complex constantPrimitive(const Real& x, const Real& alpha,
const Complex& beta,
const Complex& gamma) {
return std::exp(alpha * x) / alpha * (beta * x + gamma - beta / alpha);
}
static inline Complex linearPrimitive(const Real& x, const Real& alpha,
const Complex& beta,
const Complex& gamma) {
auto coeff = gamma - 2. * beta / alpha;
return std::exp(alpha * x) / alpha *
(beta * x * x + coeff * x - coeff / alpha);
}
static inline void computePolyCoefficients(Complex& beta, Complex& gamma,
const Real& dl, const Real& dij,
const Complex& A,
const Complex& B) {
beta = A * dl;
gamma = A * dij + B;
}
public:
template <UInt dim>
static inline void
partialIntegration(bool upper, StaticMatrix<Complex, dim, dim>& F,
const Real& alpha, const Real& dij, const Real& dl,
const StaticMatrix<Complex, dim, dim>& A,
const StaticMatrix<Complex, dim, dim>& B) {
Real xsup = (upper)? 1:0;
Real xinf = (upper)? 0:-1;
Real multiplier = (upper)? -1:1;
Complex beta = 0, gamma = 0;
F = 0;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
computePolyCoefficients(beta, gamma, dl, dij, A(i, j), B(i, j));
F(i, j) += constantPrimitive(xsup, alpha, beta, gamma) -
constantPrimitive(xinf, alpha, beta, gamma);
F(i, j) += multiplier * (linearPrimitive(xsup, alpha, beta, gamma) -
linearPrimitive(xinf, alpha, beta, gamma));
F(i, j) *= std::exp(alpha * dij / dl);
}
}
}
template <UInt dim, UInt derivative, Int yj_xi>
static inline void integrate(const Kelvin<dim, derivative>& kelvin, Real yj,
Real xi, Real dl,
StaticVector<const Real, dim - 1>& q,
StaticMatrix<Complex, dim, dim>& F) {
Complex amem[dim * dim];
Complex bmem[dim * dim];
StaticMatrix<Complex, dim, dim> A{*amem};
StaticMatrix<Complex, dim, dim> B{*bmem};
const Real q_norm = q.l2norm();
const Real dij = yj - xi;
constexpr Real sign = (yj_xi > 0) ? 1. : -1;
Real alpha = 0;
// Integrate [-1, 0]
kelvin.template computeTerms<(yj_xi > 0)>(q, A, B);
alpha = -sign * q_norm * dl;
partialIntegration(false, F, alpha, dij, dl, A, B);
// Integrate [0, 1]
if (yj_xi == 0) { // avoid computing the same A & B twice
kelvin.template computeTerms<true>(q, A, B);
alpha *= -1;
}
partialIntegration(true, F, alpha, dij, dl, A, B);
}
};
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help