stokes_kernel.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2014  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2014  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 stokes_KERNEL_HPP_INCLUDED
#define stokes_KERNEL_HPP_INCLUDED
 
#include <boost/math/constants/constants.hpp>
 
#include "../core/global_definitions.hpp"
#include "../core/kernel.hpp"
#include "../core/gaussian_quadrature.hpp"
#include "location_normal.hpp"
#include "guiggiani_1992.hpp"
 
 
namespace NiHu
{
 
class stokes_kernel
{
public:
    stokes_kernel(double mu)
        : m_mu(mu)
    {
    }
 
 
 
    double get_viscosity() const
    {
        return m_mu;
    }
 
private:
    double m_mu;
};
 
 
 
class stokes_2d_U_kernel;
 
template <>
struct kernel_traits<stokes_2d_U_kernel>
{
    typedef location_input_2d test_input_t;
    typedef location_input_2d trial_input_t;
    typedef Eigen::Matrix<double, 2, 2> result_t;
    typedef gauss_family_tag quadrature_family_t;
    static bool const is_symmetric = true;
 
    typedef asymptotic::log<1> far_field_behaviour_t;
 
    static bool const is_singular = true;
};
 
 
template <>
struct singular_kernel_traits<stokes_2d_U_kernel>
{
    typedef asymptotic::log<1> singularity_type_t;
    static unsigned const singular_quadrature_order = 7;
    typedef stokes_2d_U_kernel singular_core_t;
};
 
 
class stokes_2d_U_kernel
    : public kernel_base<stokes_2d_U_kernel>
    , public stokes_kernel
{
public:
    typedef location_input_2d::x_t x_t;
 
    stokes_2d_U_kernel(double mu)
        : stokes_kernel(mu)
    {
    }
 
    result_t operator()(x_t const &x, x_t const &y) const
    {
        using boost::math::double_constants::pi;
        double mu = get_viscosity();
        x_t rvec = y - x;   // distance vector
        double r = rvec.norm();             // scalar distance
        x_t gradr = rvec.normalized();      // gradient of the distance vector
        return ( - 1.* result_t::Identity() * (2.*std::log(r) +1.) + 2.* (gradr * gradr.transpose()) ) / (8.*pi*mu);
    }
 
    result_t operator()(
        location_input_2d const &x,
        location_input_2d const &y) const
    {
        return (*this)(x.get_x(), y.get_x());
    }
 
};
 
 
class stokes_2d_T_kernel;
 
template <>
struct kernel_traits<stokes_2d_T_kernel>
{
    typedef location_input_2d test_input_t;
    typedef location_normal_input_2d trial_input_t;
    typedef Eigen::Matrix<double, 2, 2> result_t;
    typedef gauss_family_tag quadrature_family_t;
    static bool const is_symmetric = true;
 
    typedef asymptotic::inverse<1> far_field_behaviour_t;
 
    static bool const is_singular = true;
};
 
 
template <>
struct singular_kernel_traits<stokes_2d_T_kernel>
{
    typedef asymptotic::inverse<1> singularity_type_t;
    static unsigned const singular_quadrature_order = 7;
    typedef stokes_2d_T_kernel singular_core_t;
};
 
class stokes_2d_T_kernel
    : public kernel_base<stokes_2d_T_kernel>
    , public stokes_kernel
{
public:
    typedef location_input_2d::x_t x_t;
 
    stokes_2d_T_kernel(double mu)
        : stokes_kernel(mu)
    {
    }
 
    result_t operator()(
        location_input_2d const &x,
        location_normal_input_2d const &y) const
    {
        using boost::math::double_constants::pi;
 
        x_t rvec = y.get_x() - x.get_x();
        double r = rvec.norm();
        x_t gradr = rvec.normalized();
        x_t const &n = y.get_unit_normal();
        double rdny = gradr.dot(n);
        return
          -1.0 *( rdny *( gradr*gradr.transpose()) ) / (pi*r);
    }
};
 
 
class stokes_3d_U_kernel;
 
template <>
struct kernel_traits<stokes_3d_U_kernel>
{
    typedef location_input_3d test_input_t;
    typedef location_input_3d trial_input_t;
    typedef Eigen::Matrix<double, 3, 3> result_t;
    typedef gauss_family_tag quadrature_family_t;
    static bool const is_symmetric = true;
    typedef asymptotic::inverse<1> far_field_behaviour_t;
    static bool const is_singular = true;
};
 
 
template <>
struct singular_kernel_traits<stokes_3d_U_kernel>
{
    typedef asymptotic::inverse<1> singularity_type_t;
    static unsigned const singular_quadrature_order = 7;
    typedef stokes_3d_U_kernel singular_core_t;
};
 
class stokes_3d_U_kernel
    : public kernel_base<stokes_3d_U_kernel>
    , public stokes_kernel
{
public:
    typedef location_input_3d::x_t x_t;
 
    stokes_3d_U_kernel(double mu)
        : stokes_kernel(mu)
    {
    }
 
    result_t operator()(
        x_t const& x,
        x_t const& y) const
    {
        using boost::math::double_constants::pi;
 
        double mu = get_viscosity();
        x_t rvec = y - x;
        double r = rvec.norm();
        x_t gradr = rvec.normalized();
        return ( result_t::Identity() + (gradr * gradr.transpose())) / (8. * pi * r * mu);
    }
 
    result_t operator()(
        location_input_3d const &x,
        location_input_3d const &y) const
    {
        return (*this)(x.get_x(), y.get_x());
    }
};
 
class stokes_3d_T_kernel;
 
template <>
struct kernel_traits<stokes_3d_T_kernel>
{
    typedef location_input_3d test_input_t;
    typedef location_normal_input_3d trial_input_t;
    typedef Eigen::Matrix<double, 3, 3> result_t;
    typedef gauss_family_tag quadrature_family_t;
    static bool const is_symmetric = false;
    typedef asymptotic::inverse<2> far_field_behaviour_t;
    static bool const is_singular = true;
};
 
template <>
struct singular_kernel_traits<stokes_3d_T_kernel>
{
    typedef asymptotic::inverse<2> singularity_type_t;
    static unsigned const singular_quadrature_order = 7;
    typedef stokes_3d_T_kernel singular_core_t;
};
 
class stokes_3d_T_kernel
    : public kernel_base<stokes_3d_T_kernel>
    , public stokes_kernel
{
public:
    typedef location_input_3d::x_t x_t;
 
    stokes_3d_T_kernel(double mu)
        : stokes_kernel(mu)
    {
    }
 
    result_t operator()(
        x_t const &x,
        x_t const &y,
        x_t const &ny) const
    {
        using boost::math::double_constants::pi;
 
        auto rvec = y - x;
        auto r = rvec.norm();
        auto gradr = rvec.normalized();
        auto rdny = gradr.dot(ny);
        return (-rdny * ( 3. * (gradr * gradr.transpose())) ) / (4. * pi * r * r);
    }
 
    result_t operator()(
        location_input_3d const &x,
        location_normal_input_3d const &y) const
    {
        return (*this)(x.get_x(), y.get_x(), y.get_unit_normal());
    }
};
 
template <>
class polar_laurent_coeffs<stokes_3d_T_kernel>
{
public:
    template <class guiggiani>
    static void eval(guiggiani &obj)
    {
        // in this case the traction kernel is only weakly singular
        obj.set_laurent_coeff(_m1(), guiggiani::laurent_coeff_t::Zero());
    }
};
 
}
 
#endif // STOKES_KERNEL_HPP_INCLUDED