Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/model/influence.hh
77d78e14930fmaxwell_viscoelastic…

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) {}
/// Compute leading linear term A and leading constant B
template <bool upper>
CUDA_LAMBDA inline void
computeTerms(StaticVector<Real, dim-1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& 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<true>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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<false>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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<true>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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<false>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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<true>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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<false>(StaticVector<Real, dim - 1>&& q_,
StaticMatrix<Complex, dim, dim>&& A,
StaticMatrix<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= StaticVector<const Real, dim - 1>(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__
c4science · Help