diff --git a/src/model/influence.hh b/src/model/influence.hh index 524d1b5..365eeee 100644 --- a/src/model/influence.hh +++ b/src/model/influence.hh @@ -1,292 +1,295 @@ /** * @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; 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 domain) : Influence(domain), E(E), nu(nu) {} CUDA_LAMBDA void operator()(StaticVector&& q_, StaticMatrix&& 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 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) {} private: /// Compute leading linear term A template inline void computeA(StaticVector& q_, StaticMatrix& A) const; /// Compute leading constant term B, which does not depend on y_3 inline void computeB(StaticVector& q_, StaticMatrix& B) const; public: /// Compute leading terms template inline void computeTerms(StaticVector& q_, StaticMatrix& A, StaticMatrix& B) const { computeA(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& q_, StaticMatrix& 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(1, 0) = 0; + B(0, 2) = B(1, 2) = B(2, 0) = B(2, 1) = 0; B *= (1 / (8 * mu * q * (1 - nu))); } template <> template inline void Kelvin<3, 0>::computeA(StaticVector& q_, StaticMatrix& 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 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 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.; + 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& kelvin, Real yj, Real xi, Real dl, StaticVector& q, StaticMatrix& F) { Complex amem[dim * dim]; Complex bmem[dim * dim]; StaticMatrix A{*amem}; StaticMatrix 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 = q_norm * dl; + alpha = -sign * q_norm * dl; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { - computePolyCoefficients(beta, gamma, lambda, dl, dij, A(i, j), + computePolyCoefficients(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 + if (yj_xi == 0) { // avoid computing the same A & B twice kelvin.template computeTerms(q, A, B); - alpha = q_norm * dl; + alpha *= -1; + } + for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { - computePolyCoefficients(beta, gamma, lambda, dl, dij, A(i, j), + computePolyCoefficients(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__ diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index e0164dd..13ce316 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,70 +1,76 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 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 . # ----------------------------------------------------------------------------- import tamaas as tm import numpy as np import matplotlib.pyplot as plt from pyevtk.hl import gridToVTK def dumpGrid(name, coords, field, func=lambda x: x): field_vtk = [np.zeros([coords[i].shape[0] for i in range(3)]) for j in range(3)] for i in range(3): field_vtk[i][:, :, :] = func(field[:, :, :, i]) gridToVTK(name, *coords, pointData={ 'X':field_vtk[0], 'Y':field_vtk[1], 'Z':field_vtk[2], }) -N = 32 +N = 27 -domain = [1.] * 3 + +domain = np.array([1.] * 3) + +omega = 2 * np.pi * np.array([1, 1]) / domain[:2] discretization = [N] * 3 +coords = [np.linspace(0, domain[i], discretization[i]) for i in range(3)] +x = coords[1] +y = coords[2] +x, y = np.meshgrid(x, y) model = tm.ModelFactory.createModel(tm.model_type.volume_2d, domain, discretization) displacement = model.getDisplacement() engine = tm._tamaas._test_features.Kelvin(model) source = np.zeros_like(displacement) -source[N//2, 0, 0, 2] = 1 +source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y) print("Starting computation") engine.applyVolumeForcePotential(source, displacement) print("End computation") -coords = [np.linspace(0, domain[i], discretization[i]) for i in range(3)] dumpGrid('kelvin', coords, displacement, lambda x: np.einsum('ijk->jki', x)) dumpGrid('source', coords, source, lambda x: np.einsum('ijk->jki', x)) plt.imshow(displacement[:, 4, :, 2]) plt.colorbar() #plt.show()