Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/influence.hh
a39991975caafeatures/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(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) {}
/// Compute leading linear term A and leading constant B
template <bool upper>
CUDA_LAMBDA inline void
computeTerms(VectorProxy<Real, dim-1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const;
/// Compute linear integral version of operator
template <Int yj_xi>
CUDA_LAMBDA inline void
linearIntegrate(Real yj, Real xi, Real dl, VectorProxy<Real, dim - 1>&& q,
MatrixProxy<Complex, dim, dim>&& F) 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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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>(VectorProxy<Real, dim - 1>&& q_,
MatrixProxy<Complex, dim, dim>&& A,
MatrixProxy<Complex, dim, dim>&& B) const {
q_ *= 2 * M_PI;
q_ /= VectorProxy<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)));
}
/* -------------------------------------------------------------------------- */
template <UInt dim, UInt derivative, Int yi_xi>
struct KelvinIntegrator;
inline Real kelvin_primitive(Real x, Real alpha, Real beta, Real gamma,
Real lambda) {
return std::exp(alpha * x) / alpha *
(beta * x * x + (gamma - 2 * beta / alpha) * x + lambda +
2 * beta / (alpha * alpha));
}
/// Integrate for yj > xi
template <UInt dim, UInt derivative>
struct KelvinIntegrator<dim, derivative, 1> {
CUDA_LAMBDA inline void integrate(const Kelvin<dim, derivative>& kelvin,
Real yj, Real xi, Real dl,
VectorProxy<Real, dim - 1>&& q,
MatrixProxy<Complex, dim, dim>&& F) {
Complex amem[dim * dim], bmem[dim * dim];
MatrixProxy<Complex, dim, dim> A(*amem), B(*bmem);
const Real dij = yj - xi;
kelvin.computeTerms<true>(q, A, B);
// Integrate [-1, 0] (upper domain)
Real alpha = -q * dl, beta = 0, gamma = 0, lambda = 0;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = A(i, j) * dl;
gamma = A(i, j) * (dl + dij) + B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
// Integrate [0, 1] (also upper domain)
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = -A(i, j) * dl;
gamma = A(i, j) * (dl - dij) - B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
F *= std::exp(dij);
}
};
/* -------------------------------------------------------------------------- */
/// Integrate for yj = xi
template <UInt dim, UInt derivative>
struct KelvinIntegrator<dim, derivative, 0> {
CUDA_LAMBDA inline void integrate(const Kelvin<dim, derivative>& kelvin,
Real yj, Real xi, Real dl,
VectorProxy<Real, dim - 1>&& q,
MatrixProxy<Complex, dim, dim>&& F) {
Complex amem[dim * dim], bmem[dim * dim];
MatrixProxy<Complex, dim, dim> A(*amem), B(*bmem);
const Real dij = yj - xi;
kelvin.computeTerms<false>(q, A, B);
// Integrate [-1, 0] (lower domain)
Real alpha = q * dl, beta = 0, gamma = 0, lambda = 0;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = A(i, j) * dl;
gamma = A(i, j) * (dl + dij) + B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
// Integrate [0, 1] (upper domain)
kelvin.computeTerms<true>(q, A, B);
alpha = -q * dl;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = -A(i, j) * dl;
gamma = A(i, j) * (dl - dij) - B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
F *= std::exp(dij);
}
};
/* -------------------------------------------------------------------------- */
/// Integrate for yj > xi
template <UInt dim, UInt derivative>
struct KelvinIntegrator<dim, derivative, -1> {
CUDA_LAMBDA inline void integrate(const Kelvin<dim, derivative>& kelvin,
Real yj, Real xi, Real dl,
VectorProxy<Real, dim - 1>&& q,
MatrixProxy<Complex, dim, dim>&& F) {
Complex amem[dim * dim], bmem[dim * dim];
MatrixProxy<Complex, dim, dim> A(*amem), B(*bmem);
const Real dij = yj - xi;
kelvin.computeTerms<false>(q, A, B);
// Integrate [-1, 0] (lower domain)
Real alpha = q * dl, beta = 0, gamma = 0, lambda = 0;
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = A(i, j) * dl;
gamma = A(i, j) * (dl + dij) + B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
// Integrate [0, 1] (also lower domain)
for (UInt i = 0; i < dim; ++i) {
for (UInt j = 0; j < dim; ++j) {
beta = -A(i, j) * dl;
gamma = A(i, j) * (dl - dij) - B(i, j);
lambda = A(i, j) * dij + B(i, j);
F(i, j) = kelvin_primitive(0, alpha, beta, gamma, lambda) -
kelvin_primitive(-1, alpha, beta, gamma, lambda);
}
}
F *= std::exp(dij);
}
};
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help