Page MenuHomec4science
Files F65854822

influence.hh
No OneTemporary
Actions

File Metadata

Created
Thu, Jun 6, 16:22

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)));
}
/* -------------------------------------------------------------------------- */
template <UInt dim, UInt derivative, Int yj_xi>
class KelvinIntegrator {
static inline Complex primitive(const Real& x, const Real& alpha,
const Complex& beta, const Complex& gamma,
const Complex& lambda) {
return std::exp(alpha * x) / alpha *
(beta * x * x + (gamma - 2. * beta / alpha) * x + lambda +
2. * beta / (alpha * alpha));
}
template <bool upper>
static inline void computePolyCoefficients(Complex& beta, Complex& gamma,
Complex& lambda, const Real& dl,
const Real& dij, const Complex& A,
const Complex& B) {
constexpr Real sign = (upper) ? -1. : 1.;
beta = sign * A * dl;
gamma = A * (dl + sign * dij) + sign * B;
lambda = A * dij + B;
}
public:
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;
Complex beta = 0, gamma = 0, lambda = 0;
// Integrate [-1, 0]
kelvin.template computeTerms<(yj_xi > 0)>(q, A, B);
alpha = -sign * q_norm * dl;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
computePolyCoefficients<false>(beta, gamma, lambda, dl, dij, A(i, j),
B(i, j));
F(i, j) = primitive(0, alpha, beta, gamma, lambda) -
primitive(-1, alpha, beta, gamma, lambda);
}
}
// Integrate [0, 1]
if (yj_xi == 0) { // avoid computing the same A & B twice
kelvin.template computeTerms<true>(q, A, B);
alpha *= -1;
}
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
computePolyCoefficients<true>(beta, gamma, lambda, dl, dij, A(i, j),
B(i, j));
F(i, j) = primitive(1, alpha, beta, gamma, lambda) -
primitive(0, alpha, beta, gamma, lambda);
}
}
}
};
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__

Event Timeline

c4science · Help