diff --git a/src/model/influence.hh b/src/model/influence.hh index 235968d..eb66d98 100644 --- a/src/model/influence.hh +++ b/src/model/influence.hh @@ -1,142 +1,388 @@ /** * @file * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "loop.hh" #include "static_types.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ namespace influence { /// Mother class for influence functions template class Influence { public: Influence(StaticVector domain) : domain(domain) {} protected: StaticVector domain; }; /// Class representing the basic influence functions from Stanley & Kato (1997) template class Basic : public Influence { public: Basic(Real E_hertz, StaticVector domain) : Influence(domain), E_hertz(E_hertz) {} /// Influence function: computes the influence coefficient F CUDA_LAMBDA void operator()(StaticVector&& 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 class Surface : public Influence {}; /// 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 domain) : Influence(domain), E(E), nu(nu) {} CUDA_LAMBDA void operator()(StaticVector&& q_, StaticMatrix&& F) const { q_ *= 2 * M_PI; q_ /= this->domain; - Real q = q_.l2norm(); + 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; - Complex I = Complex(0, 1); + 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 domain) : Influence(domain), E(E), nu(nu) {} CUDA_LAMBDA void operator()(StaticVector&& q_, StaticMatrix&& F) const { q_ *= 2 * M_PI; q_ /= this->domain; - Real q = q_.l2norm(); + 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; - Complex I = Complex(0, 1); + const Complex I = Complex(0, 1); }; + +/// Class for the Kelvin tensor +template +class Kelvin; + +template +class Kelvin<3, derivative_order> : Influence<3> { + constexpr static UInt dim = 3; +public: + Kelvin(Real mu, Real nu, StaticVector domain): + Influence(domain), mu(mu), nu(nu) {} + + /// Compute leading linear term A and leading constant B + template + CUDA_LAMBDA inline void + computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const; + +private: + Real mu, nu; + const Complex I = Complex(0, 1); +}; + +/* -------------------------------------------------------------------------- */ + +/// Kelvin tensor for upper domain +template <> +template <> +inline void +Kelvin<3, 0>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = -bq1 * bq1; + B(0, 0) = 3 * bq1 * bq1 + 4 * bq2 * bq2 - 4 * nu; + + A(0, 1) = A(1, 0) = -bq1 * bq2; + B(0, 1) = B(1, 0) = -bq1 * bq2; + + A(0, 2) = A(2, 0) = -I * bq1; + B(0, 2) = B(2, 0) = 0; + + A(1, 1) = -bq2 * bq2; + B(1, 1) = 4 * bq1 * bq1 + 3 * bq2 * bq2 - 4 * nu; + + A(1, 2) = A(2, 1) = -I * bq2; + B(1, 2) = B(2, 1) = 0; + + A(2, 2) = 1; + B(2, 2) = 3 - 4 * nu; + + A *= (1 / (8 * mu * (1 - nu))); + B *= (1 / (8 * mu * q * (1 - nu))); +} + +/// Kelvin tensor for lower domain +template <> +template <> +inline void +Kelvin<3, 0>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = bq1 * bq1; + B(0, 0) = 3 * bq1 * bq1 + 4 * bq2 * bq2 - 4 * nu; + + A(0, 1) = A(1, 0) = bq1 * bq2; + B(0, 1) = B(1, 0) = -bq1 * bq2; + + A(0, 2) = A(2, 0) = -I * bq1; + B(0, 2) = B(2, 0) = 0; + + A(1, 1) = bq2 * bq2; + B(1, 1) = 4 * bq1 * bq1 + 3 * bq2 * bq2 - 4 * nu; + + A(1, 2) = A(2, 1) = -I * bq2; + B(1, 2) = B(2, 1) = 0; + + A(2, 2) = -1; + B(2, 2) = 3 - 4 * nu; + + A *= (1 / (8 * mu * (1 - nu))); + B *= (1 / (8 * mu * q * (1 - nu))); +} + +/* -------------------------------------------------------------------------- */ + +/// First derivative of Kelvin tensor (upper domain) +template <> +template <> +inline void +Kelvin<3, 1>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = bq1 * bq1; + B(0, 0) = -4 * (1 - nu); + + A(0, 1) = A(1, 0) = bq1 * bq2; + B(0, 1) = B(1, 0) = 0; + + A(0, 2) = A(2, 0) = I * bq1; + B(0, 2) = B(2, 0) = - I * bq1; + + A(1, 1) = bq2 * bq2; + B(1, 1) = -4 * (1 - nu); + + A(1, 2) = A(2, 1) = I * bq2; + B(1, 2) = B(2, 1) = -I * bq2; + + A(2, 2) = -1; + B(2, 2) = 4 * nu - 2; + + A *= (q / (8 * mu * (1 - nu))); + B *= (1 / (8 * mu * (1 - nu))); +} + +/// First derivative of Kelvin tensor (lower domain) +template <> +template <> +inline void +Kelvin<3, 1>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = bq1 * bq1; + B(0, 0) = 4 * (1 - nu); + + A(0, 1) = A(1, 0) = bq1 * bq2; + B(0, 1) = B(1, 0) = 0; + + A(0, 2) = A(2, 0) = - I * bq1; + B(0, 2) = B(2, 0) = - I * bq1; + + A(1, 1) = bq2 * bq2; + B(1, 1) = 4 * (1 - nu); + + A(1, 2) = A(2, 1) = -I * bq2; + B(1, 2) = B(2, 1) = -I * bq2; + + A(2, 2) = -1; + B(2, 2) = -(4 * nu - 2); + + A *= (q / (8 * mu * (1 - nu))); + B *= (1 / (8 * mu * (1 - nu))); +} + +/* -------------------------------------------------------------------------- */ + +/// Regular second derivative of Kelvin tensor (upper domain) +template <> +template <> +inline void +Kelvin<3, 2>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = -bq1 * bq1; + B(0, 0) = bq1 * bq1 - 4. * nu + 4.; + + A(0, 1) = A(1, 0) = -bq1 * bq2; + B(0, 1) = B(1, 0) = bq1 * bq2; + + A(0, 2) = A(2, 0) = -I * bq1; + B(0, 2) = B(2, 0) = 2. * I * bq1; + + A(1, 1) = -bq2 * bq2; + B(1, 1) = bq2 * bq2 - 4. * nu + 4.; + + A(1, 2) = A(2, 1) = -I * bq2; + B(1, 2) = B(2, 1) = 2. * I * bq2; + + A(2, 2) = 1.; + B(2, 2) = 1. - 4. * nu; + + A *= (q * q / (8. * mu * (1. - nu))); + B *= (q / (8. * mu * (1. - nu))); +} + +/// Regular second derivative of Kelvin tensor (lower domain) +template <> +template <> +inline void +Kelvin<3, 2>::computeTerms(StaticVector&& q_, + StaticMatrix&& A, + StaticMatrix&& B) const { + q_ *= 2 * M_PI; + q_ /= StaticVector(this->domain(0)); + + const Real q = q_.l2norm(); + const Real bq1 = q_(0) / q; + const Real bq2 = q_(1) / q; + + A(0, 0) = bq1 * bq1; + B(0, 0) = bq1 * bq1 - 4. * nu + 4.; + + A(0, 1) = A(1, 0) = bq1 * bq2; + B(0, 1) = B(1, 0) = bq1 * bq2; + + A(0, 2) = A(2, 0) = -I * bq1; + B(0, 2) = B(2, 0) = -2. * I * bq1; + + A(1, 1) = bq2 * bq2; + B(1, 1) = bq2 * bq2 - 4. * nu + 4.; + + A(1, 2) = A(2, 1) = -I * bq2; + B(1, 2) = B(2, 1) = -2. * I * bq2; + + A(2, 2) = -1.; + B(2, 2) = 1. - 4. * nu; + + A *= (q * q / (8. * mu * (1. - nu))); + B *= (q / (8. * mu * (1. - nu))); +} + } // namespace influence /* -------------------------------------------------------------------------- */ __END_TAMAAS__