Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62512897
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
Mon, May 13, 17:39
Size
9 KB
Mime Type
text/x-c++
Expires
Wed, May 15, 17:39 (2 d)
Engine
blob
Format
Raw Data
Handle
17657020
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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
/// Mother class for influence functions
template <UInt dim>
class Influence {
public:
Influence(VectorProxy<const Real, dim> domain) : domain(domain) {}
protected:
VectorProxy<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, VectorProxy<const Real, dim> domain)
: Influence<dim>(domain), E_hertz(E_hertz) {}
/// Influence function: computes the influence coefficient F
CUDA_LAMBDA void operator()(VectorProxy<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, VectorProxy<const Real, dim> domain)
: Influence<dim>(domain), E(E), nu(nu) {}
CUDA_LAMBDA void
operator()(VectorProxy<Real, dim>&& q_,
MatrixProxy<Complex, dim + 1, dim + 1>&& F) const {
q_ *= 2 * M_PI;
q_ /= this->domain;
const Real q = q_.l2norm();
q_ /= q;
// Diagonal
F(0, 0) = 2 * (1 + nu) * (1 - nu * q_(0) * q_(0));
F(1, 1) = 2 * (1 + nu) * (1 - nu * q_(1) * q_(1));
F(2, 2) = 2 * (1 - nu * nu);
F(0, 1) = F(1, 0) = -q_(0) * q_(1) * 2 * nu * (1 + nu);
F(0, 2) = I * q_(0) * (1 + nu) * (1 - 2 * nu);
F(1, 2) = I * q_(1) * (1 + nu) * (1 - 2 * nu);
F(2, 0) = -F(0, 2);
F(2, 1) = -F(1, 2);
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, VectorProxy<const Real, dim> domain)
: Influence<dim>(domain), E(E), nu(nu) {}
CUDA_LAMBDA void
operator()(VectorProxy<Real, dim>&& q_,
MatrixProxy<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, VectorProxy<const Real, dim> domain)
: Influence<dim>(domain), mu(mu), nu(nu) {}
private:
/// Compute leading linear term A
template <bool upper>
inline void computeA(VectorProxy<const Real, dim - 1>& q_,
Matrix<Complex, dim, dim>& A) const;
/// Compute leading constant term B, which does not depend on y_3
inline void computeB(VectorProxy<const Real, dim - 1>& q_,
Matrix<Complex, dim, dim>& B) const;
public:
/// Compute leading terms
template <bool upper>
inline void computeTerms(VectorProxy<const Real, dim - 1>& q_,
Matrix<Complex, dim, dim>& A,
Matrix<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(VectorProxy<const Real, dim - 1>& q_,
Matrix<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) = 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(VectorProxy<const Real, dim - 1>& q_,
Matrix<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) = sign;
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 <bool upper, UInt dim>
static inline void partialIntegration(Matrix<Complex, dim, dim>& F,
const Real& alpha, const Real& dij,
const Real& dl,
const Matrix<Complex, dim, dim>& A,
const Matrix<Complex, dim, dim>& B) {
Real xsup = (upper) ? 1 : 0;
Real xinf = (upper) ? 0 : -1;
Real multiplier = (upper) ? -1 : 1;
Real e_t = std::exp(alpha * dij / dl);
Complex beta = 0, gamma = 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) += e_t * (constantPrimitive(xsup, alpha, beta, gamma) -
constantPrimitive(xinf, alpha, beta, gamma));
F(i, j) +=
e_t * multiplier * (linearPrimitive(xsup, alpha, beta, gamma) -
linearPrimitive(xinf, alpha, beta, gamma));
}
}
}
template <UInt dim, UInt derivative, Int yj_xi>
static inline void integrate(const Kelvin<dim, derivative>& kelvin, Real yj,
Real xi, Real dl,
VectorProxy<const Real, dim - 1>& q,
Matrix<Complex, dim, dim>& F) {
F = 0;
Matrix<Complex, dim, dim> A, B;
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__
Event Timeline
Log In to Comment