helmholtz_3d_singular_integrals.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2018  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2018  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 NIHU_HELMHOLTZ_3D_SINGULAR_INTEGRALS_HPP_INCLUDED
#define NIHU_HELMHOLTZ_3D_SINGULAR_INTEGRALS_HPP_INCLUDED
 
#include <boost/math/constants/constants.hpp>
 
#include "helmholtz_kernel.hpp"
#include "field_type_helpers.hpp"
#include "guiggiani_1992.hpp"
#include "lib_element.hpp"
#include "laplace_3d_singular_integrals.hpp"
#include "normal_derivative_singular_integrals.hpp"
#include "plane_element_helper.hpp"
#include "../core/match_types.hpp"
#include "../core/singular_integral_shortcut.hpp"
#include "../util/math_functions.hpp"
 
#if NIHU_MEX_DEBUGGING
#include <mex.h>
#endif
 
namespace NiHu
{
 
template <unsigned order>
class helmholtz_3d_SLP_collocation_constant_plane
{
private:
    template <class wavenumber_t>
    static std::complex<double> dynamic_part(double const &r, wavenumber_t const &k)
    {
        using namespace std::complex_literals;
        return -1.i*k * std::exp(-1.i*k*r / 2.0) * sinc(k*r / 2.0);
    }
 
public:
    template <class elem_t, class wavenumber_t>
    static std::complex<double> eval(
        elem_t const &elem,
        typename elem_t::x_t const &x0,
        wavenumber_t const &k)
    {
        using namespace boost::math::double_constants;
        
        typedef regular_quad_store<typename elem_t::domain_t, order> quadr_t;
 
        enum { N = elem_t::domain_t::num_corners };
 
        double r[N], theta[N], alpha[N];
        plane_element_helper(elem, x0, r, theta, alpha);
 
        // integrate static part analytically
        double I_stat = 0.0;
        for (unsigned i = 0; i < N; ++i)
            I_stat += r[i] * std::sin(alpha[i]) *
            std::log(std::tan((alpha[i] + theta[i]) / 2.0) / std::tan(alpha[i] / 2.0)
            );
 
        // integrate dynamic part
        std::complex<double> I_dyn = 0.0;
        for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
        {
            double r = (elem.get_x(it->get_xi()) - x0).norm();
            double jac = elem.get_normal(it->get_xi()).norm();
            I_dyn += dynamic_part(r, k) * it->get_w() * jac;
        }
 
        // assemble result from static and dynamic parts
        return (I_stat + I_dyn) / (4.0 * pi);
    }
};
 
template <unsigned order>
class helmholtz_3d_HSP_collocation_constant_plane
{
private:
    template <class wavenumber_t>
    static std::complex<double> dynamic_part(double const &r, wavenumber_t const &k)
    {
        using namespace std::complex_literals;
        std::complex<double> const ikr(1.i*k*r);
        if (std::abs(r) > 1e-3)
            return (std::exp(-ikr)*(1.0 + ikr) - 1.0 + ikr*ikr / 2.0) / r / r / r;
        else
            return -1.i*k*k*k * (
            1.0 / 3.0 - ikr*(1.0 / 8.0 - ikr*(1.0 / 30.0 - ikr*(1.0 / 144.0 - ikr*(1.0 / 840.0 - ikr / 5760.0))))
            );
    }
 
public:
    template <class elem_t, class wavenumber_t>
    static std::complex<double> eval(
        elem_t const &elem,
        typename elem_t::x_t const &x0,
        wavenumber_t const &k)
    {
        using namespace boost::math::double_constants;
        
        typedef regular_quad_store<typename elem_t::domain_t, order> quadr_t;
        enum { N = elem_t::domain_t::num_corners };
 
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("Integrating Helmholtz HSP over constant plane, N = %d\n", N);
            printed = true;
        }
#endif
 
        double r[N], theta[N], alpha[N];
        plane_element_helper(elem, x0, r, theta, alpha);
 
        // integrate static part
        double IG0 = 0.0, IddG0 = 0.0;
        for (unsigned i = 0; i < N; ++i)
        {
            // integral is zero for highly distorted elements
            if (theta[i] < 1e-3)
                continue;
            IG0 += r[i] * std::sin(alpha[i]) *
                std::log(std::tan((alpha[i] + theta[i]) / 2.0) / tan(alpha[i] / 2.0));
            IddG0 += (std::cos(alpha[i] + theta[i]) - std::cos(alpha[i])) / (r[i] * std::sin(alpha[i]));
        }
 
        // integrate dynamic part
        std::complex<double> I_acc = 0.0;
        for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
        {
            double r = (elem.get_x(it->get_xi()) - x0).norm();
            double jac = elem.get_normal(it->get_xi()).norm();
            I_acc += dynamic_part(r, k) * it->get_w() * jac;
        }
        
        // assemble result from static and dynamic parts
        return (IddG0 + k*k / 2.0 * IG0 + I_acc) / (4.0 * pi);
    }
};
 
template <class WaveNumber, class TestField, class TrialField>
class singular_integral_shortcut<
    helmholtz_3d_SLP_kernel<WaveNumber>, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
        is_collocational<TestField, TrialField>::value &&
        is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_3d_SLP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = helmholtz_3d_SLP_collocation_constant_plane<7>::eval(
            trial_field.get_elem(),
            trial_field.get_elem().get_center(),
            kernel.derived().get_wave_number());
        return result;
    }
};
 
 
template <class WaveNumber, class TestField, class TrialField>
class singular_integral_shortcut<
    helmholtz_3d_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
        is_collocational<TestField, TrialField>::value &&
        is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_3d_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = helmholtz_3d_HSP_collocation_constant_plane<7>::eval(
            trial_field.get_elem(),
            trial_field.get_elem().get_center(),
            kernel.derived().get_wave_number());
        return result;
    }
};
 
template <class WaveNumber, class TestField, class TrialField>
class singular_integral_shortcut<
    helmholtz_3d_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
        is_collocational<TestField, TrialField>::value &&
        !is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_3d_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("Calling Guiggiani now\n");
            printed = true;
        }
#endif
        typedef guiggiani<TrialField, helmholtz_3d_HSP_kernel<WaveNumber>, 5, 9> guiggiani_t;
        auto const &elem = trial_field.get_elem();
        guiggiani_t gui(elem, kernel.derived());
        for (unsigned r = 0; r < TestField::num_dofs; ++r)
        {
            auto const &xi0 = TestField::nset_t::corner_at(r);
            gui.integrate(result.row(r), xi0, elem.get_normal(xi0));
        }
        
        return result;
    }
};
 
} // end of namespace NiHu
 
 
#endif // NIHU_HELMHOLTZ_3D_SINGULAR_INTEGRALS_HPP_INCLUDED