elastodynamics_kernel.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2015  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2015  Peter Rucz <rucz@hit.bme.hu>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
#ifndef ELASTODYNAMICS_KERNEL_HPP_INCLUDED
#define ELASTODYNAMICS_KERNEL_HPP_INCLUDED
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/bessel.hpp>
 
#include "../core/global_definitions.hpp"
#include "../core/kernel.hpp"
#include "../core/gaussian_quadrature.hpp"
#include "location_normal.hpp"
#include "elastostatics_kernel.hpp"
 
namespace NiHu
{
 
class elastodynamics_data
{
public:
    elastodynamics_data(double nu, double rho, double mu, double omega)
        : m_nu(nu), m_rho(rho), m_mu(mu), m_omega(omega) {}
    double get_poisson_ratio() const { return m_nu; }
    double get_mass_density() const { return m_rho; }
    double get_shear_modulus() const { return m_mu; }
    double get_frequency() const { return m_omega; }
 
private:
    double m_nu;    
    double m_rho;   
    double m_mu;    
    double m_omega; 
};
 
class elastodynamics_3d_U_kernel;
 
template <>
struct kernel_traits<elastodynamics_3d_U_kernel>
{
    using test_input_t = location_input_3d;
    using trial_input_t = location_input_3d ;
    using result_t = Eigen::Matrix<std::complex<double>, 3, 3>;
    using quadrature_family_t = gauss_family_tag;
    static bool const is_symmetric = true;
    using far_field_behaviour_t = asymptotic::inverse<1>;
    static bool const is_singular = true;
};
 
template <>
struct singular_kernel_traits<elastodynamics_3d_U_kernel>
{
    using singularity_type_t = asymptotic::inverse<1>;
    static unsigned const singular_quadrature_order = 7;
    using singular_core_t = elastostatics_3d_U_kernel;
};
 
class elastodynamics_3d_U_kernel
    : public kernel_base<elastodynamics_3d_U_kernel>
    , public elastodynamics_data
{
public:
    elastodynamics_3d_U_kernel(double nu, double rho, double mu, double omega)
        : elastodynamics_data(nu, rho, mu, omega)
    {
    }
 
    result_t operator()(
        location_input_3d const &x,
        location_input_3d const &y) const
    {
        using namespace std::complex_literals;
        using namespace boost::math::double_constants;
    
        // compute distance and its gradient
        auto rvec = y.get_x() - x.get_x();
        auto r = rvec.norm();
        auto gradr = rvec.normalized();
 
        // material properties and Helmholtz numbers
        auto nu = get_poisson_ratio();
        auto a2 = (1.-2.*nu)/2./(1.-nu);
        auto a = std::sqrt(a2);
        auto mu = get_shear_modulus();
        auto rho = get_mass_density();
        auto cS = std::sqrt(mu/rho);
 
        auto om = get_frequency();
        auto kSr = om / cS * r;
        auto kPr = a * kSr;
 
        std::complex<double> psi =
            std::exp(-1.i*kPr) * a2 * (1.i/kPr + 1./(kPr*kPr)) +
            std::exp(-1.i*kSr) * (1. - 1.i/kSr - 1./(kSr*kSr));
 
        std::complex<double> chi =
            std::exp(-1.i*kPr) * a2 * (1. - 3.i/kPr - 3./(kPr*kPr)) -
            std::exp(-1.i*kSr) * (1. - 3.i/kSr - 3./(kSr*kSr));
 
        // matrix valued result
        return ( psi * result_t::Identity() + chi * (gradr * gradr.transpose()) ) / (4.*pi*mu*r);
    }
}; // end of class elastodynamics_3d_U_kernel
 
 
 
class elastodynamics_3d_T_kernel;
 
template <>
struct kernel_traits<elastodynamics_3d_T_kernel>
{
    using test_input_t = location_input_3d;
    using trial_input_t = location_normal_input_3d;
    using result_t = Eigen::Matrix<std::complex<double>, 3, 3>;
    using quadrature_family_t = gauss_family_tag;
    static bool const is_symmetric = false;
    using far_field_behaviour_t = asymptotic::inverse<2>;
    static bool const is_singular = true;
};
 
template <>
struct singular_kernel_traits<elastodynamics_3d_T_kernel>
{
    using singularity_type_t = asymptotic::inverse<2>;
    static unsigned const singular_quadrature_order = 7;
    using singular_core_t = elastostatics_3d_T_kernel;
};
 
class elastodynamics_3d_T_kernel
    : public kernel_base<elastodynamics_3d_T_kernel>
    , public elastodynamics_data
{
public:
    elastodynamics_3d_T_kernel(double nu, double rho, double mu, double omega)
        : elastodynamics_data(nu, rho, mu, omega)
    {
    }
        
    result_t operator()(
        location_input_3d const &x,
        location_normal_input_3d const &y) const
    {
        using namespace std::complex_literals;
        using namespace boost::math::double_constants;
 
        // compute distance and its gradient
        auto rvec = y.get_x() - x.get_x();
        auto r = rvec.norm();
        auto const &n = y.get_unit_normal();
        auto gradr = rvec.normalized();
        double rdn = gradr.dot(n);
 
        // material properties and Helmholtz numbers
        auto nu = get_poisson_ratio();
        auto a2 = (1.-2.*nu)/2./(1.-nu);
        auto a = std::sqrt(a2);
        auto mu = get_shear_modulus();
        auto lambdapermu = (1./a2 - 2.);
        auto rho = get_mass_density();
        auto cS = std::sqrt(mu/rho);
 
        auto om = get_frequency();
        auto kSr = om / cS * r;
        auto kPr = a * kSr;
 
        std::complex<double> psi =
            std::exp(-1.i*kPr) * a2 * (     1.i/kPr + 1./(kPr*kPr)) +
            std::exp(-1.i*kSr) *      (1. - 1.i/kSr - 1./(kSr*kSr));
 
        std::complex<double> chi =
            std::exp(-1.i*kPr) * a2 * (1. - 3.i/kPr - 3./(kPr*kPr)) -
            std::exp(-1.i*kSr) *      (1. - 3.i/kSr - 3./(kSr*kSr));
 
        std::complex<double> dpsi =
            std::exp(-1.i*kPr) * a2 * (1.         - 2.i/kPr - 2./(kPr*kPr)) -
            std::exp(-1.i*kSr) *      (1. + 1.i*kSr - 2.i/kSr - 2./(kSr*kSr));
 
        std::complex<double> dchi =
            std::exp(-1.i*kSr) *      (3. + 1.i*kSr - 6.i/kSr - 6./(kSr*kSr)) -
            std::exp(-1.i*kPr) * a2 * (3. + 1.i*kPr - 6.i/kPr - 6./(kPr*kPr));
 
        std::complex<double> A = dpsi - psi;
        std::complex<double> B = dchi - 3.*chi;
        std::complex<double> C = chi;
 
        // matrix valued result
        return (
            rdn * (result_t::Identity() * (A+C) + 2.*B*(gradr * gradr.transpose()))
            +
            (2.*C + lambdapermu*(A+B+4.*C)) * (gradr * n.transpose())
            +
            (A+C) * (n * gradr.transpose())
        ) / (4.*pi*r*r);
    }
}; // end of class elastodynamics_3d_T_kernel
 
 
 
class elastodynamics_2d_U_kernel;
 
template <>
struct kernel_traits<elastodynamics_2d_U_kernel>
{
    using test_input_t = location_input_2d;
    using trial_input_t = location_input_2d ;
    using result_t = Eigen::Matrix<std::complex<double>, 2, 2>;
    using quadrature_family_t = gauss_family_tag;
    static bool const is_symmetric = true;
    using far_field_behaviour_t = asymptotic::inverse<1>;   // \todo check this
    static bool const is_singular = true;
};
 
template <>
struct singular_kernel_traits<elastodynamics_2d_U_kernel>
{
    using singularity_type_t = asymptotic::log<1>;
    static unsigned const singular_quadrature_order = 7;
    using singular_core_t = elastostatics_2d_U_kernel;
};
 
class elastodynamics_2d_U_kernel
    : public kernel_base<elastodynamics_2d_U_kernel>
    , public elastodynamics_data
{
public:
    elastodynamics_2d_U_kernel(double nu, double rho, double mu, double omega)
        : elastodynamics_data(nu, rho, mu, omega)
    {
    }
 
    result_t operator()(
        location_input_2d const &x,
        location_input_2d const &y) const
    {
        using namespace std::complex_literals;
        using namespace boost::math::double_constants;
        using boost::math::cyl_hankel_2;
 
        // compute distance and its gradient
        auto rvec = y.get_x() - x.get_x();
        auto r = rvec.norm();
        auto gradr = rvec.normalized();
 
        // material properties and Helmholtz numbers
        auto nu = get_poisson_ratio();
        auto a2 = (1.-2.*nu)/2./(1.-nu);
        auto a = std::sqrt(a2);
        auto mu = get_shear_modulus();
        auto rho = get_mass_density();
        auto cS = std::sqrt(mu/rho);
 
        auto om = get_frequency();
        auto kSr = om / cS * r;
        auto kPr = a * kSr;
 
        std::complex<double> psi = .25i * ((cyl_hankel_2(1, kSr)/kSr - a2 * cyl_hankel_2(1, kPr) / kPr) - cyl_hankel_2(0, kSr));
 
        std::complex<double> chi = .25i * (a2 *cyl_hankel_2(2, kPr) - cyl_hankel_2(2, kSr));
 
        // matrix valued result
        return ( psi * result_t::Identity() + chi * (gradr * gradr.transpose()) ) / mu;
    }
}; // end of class elastodynamics_2d_U_kernel
 
} // end of namespace NiHu
 
#endif // ELASTODYNAMICS_KERNEL_HPP_INCLUDED