helmholtz_2d_singular_double_integrals.hpp

// 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_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED
#define NIHU_HELMHOLTZ_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED
 
#include <boost/math/constants/constants.hpp>
 
#include "helmholtz_kernel.hpp"
#include "field_type_helpers.hpp"
#include "lib_element.hpp"
#include "laplace_2d_singular_double_integrals.hpp"
#include "normal_derivative_singular_integrals.hpp"
#include "../core/match_types.hpp"
#include "../core/singular_integral_shortcut.hpp"
#include "../util/math_functions.hpp"
 
namespace NiHu
{
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_SLP_galerkin_face_general
{
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef typename test_field_t::nset_t test_shape_t;
    typedef typename trial_field_t::nset_t trial_shape_t;
 
    static size_t const nTest = test_shape_t::num_nodes;
    static size_t const nTrial = trial_shape_t::num_nodes;
 
    typedef Eigen::Matrix<std::complex<double>, nTest, nTrial> result_t;
 
    typedef typename trial_field_t::elem_t elem_t;
    typedef typename elem_t::domain_t domain_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename elem_t::x_t x_t;
 
    typedef regular_quad_store<domain_t, order> quadrature_t;
 
public:
    template <class wave_number_t>
    static result_t eval(elem_t const &elem, wave_number_t const &k)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_SLP_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
        typedef helmholtz_2d_SLP_kernel<wave_number_t> kernel_t;
        kernel_t kernel(k);
 
        result_t I = result_t::Zero();
 
        // loop over u
        for (auto itu = quadrature_t::quadrature.begin(); itu != quadrature_t::quadrature.end(); ++itu) {
            xi_t u = itu->get_xi();
            double wu = itu->get_w();
 
            result_t NNt0 = test_shape_t::template eval_shape<0>(u) *
                trial_shape_t::template eval_shape<0>(u).transpose();
            double J0 = elem.get_dx(u).norm();
 
            // loop over v
            for (auto itv = quadrature_t::quadrature.begin(); itv != quadrature_t::quadrature.end(); ++itv) {
                double v = (itv->get_xi()(0) + 1.) / 2.;
                double wv = itv->get_w() / 2.;
 
                xi_t xi, eta;
 
                // loop over domains (0) upper triangle, (1) lower triangle 
                for (int d = 0; d < 2; ++d) {
                    if (d == 0) {
                        xi(0) = u(0) + (+1. - u(0)) * v;
                        eta(0) = u(0) + (-1. - u(0)) * v;
                    }
                    else {
                        xi(0) = u(0) + (-1. - u(0)) * v;
                        eta(0) = u(0) + (+1. - u(0)) * v;
                    }
 
                    std::complex<double> G = kernel(elem.get_x(xi), elem.get_x(eta));
                    double Jxi = elem.get_dx(xi).norm();
                    double Jeta = elem.get_dx(eta).norm();
                    auto Nxi = test_shape_t::template eval_shape<0>(xi);
                    auto Neta = trial_shape_t::template eval_shape<0>(eta);
 
                    // evaluate integrand
                    result_t F = G * Jxi * Jeta * (Nxi * Neta.transpose());
 
                    // evaluate integrand's singular part
                    double G0 = -std::log(2. * v * J0) / two_pi;
                    result_t F0 = G0 * J0 * J0 * NNt0;
 
                    // integrate difference numerically
                    I += (F - F0) * (2. * (1. - v)) * wu * wv;
                } // loop over domains
            } // loop over v
 
            // compensate with analytical inner integral of the regular part
            I += 2. * (1.5 - std::log(2. * J0)) / two_pi * J0 * J0 * NNt0 * wu;
        } // loop over u
 
        return I;
    } // function eval
}; // class helmholtz_2d_SLP_galerkin_face_general
 
 
template <unsigned expansion_length>
class helmholtz_2d_SLP_galerkin_face_constant_line
{
public:
    template <class wavenumber_t>
    static std::complex<double> eval(double R, wavenumber_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_SLP_galerkin_face_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
        using namespace std::complex_literals;
    
        wavenumber_t kR = k * R; 
        wavenumber_t logkR = std::log(kR);
        
        auto I = 1. - 2.i/pi * (logkR - 1.5 + euler);
        
        wavenumber_t q = -kR * kR;
        wavenumber_t pow = 1.;
        double Cn = 0.;
        
        for (unsigned n = 1; n <= expansion_length; ++n) {
            unsigned d = (n+1) * (2*n+1);
            double Fn = 1./d;
            wavenumber_t Gn = logkR/d - (4*n + 3)/2./(d*d);
            pow *= q/(n*n);
            Cn += 1./n;
            I += (Fn-(2.i/pi)*(Gn+(euler - Cn)*Fn)) * pow;
        }
        
        return I * (R*R) * -1.0i;
    } // end of function eval
}; // end of class helmholtz_2d_SLP_galerkin_face_constant_line
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_DLP_galerkin_face_general
{
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef typename test_field_t::nset_t test_shape_t;
    typedef typename trial_field_t::nset_t trial_shape_t;
 
    static size_t const nTest = test_shape_t::num_nodes;
    static size_t const nTrial = trial_shape_t::num_nodes;
 
    typedef Eigen::Matrix<std::complex<double>, nTest, nTrial> result_t;
 
    typedef typename trial_field_t::elem_t elem_t;
    typedef typename elem_t::domain_t domain_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename elem_t::x_t x_t;
 
    typedef regular_quad_store<domain_t, order> quadrature_t;
 
public:
    template <class wave_number_t>
    static result_t eval(elem_t const &elem, wave_number_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_DLP_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        typedef helmholtz_2d_DLP_kernel<wave_number_t> kernel_t;
        kernel_t kernel(k);
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto itv = quadrature_t::quadrature.begin(); itv != quadrature_t::quadrature.end(); ++itv) {
            xi_t v = itv->get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = itv->get_w() / 2.;
 
            // loop over u
            for (auto itu = quadrature_t::quadrature.begin(); itu != quadrature_t::quadrature.end(); ++itu) {
                xi_t u = itu->get_xi();
                double wu = itu->get_w();
 
                // loop over domains
                for (int d = 0; d < 2; ++d) {
                    xi_t xi, eta;
 
                    if (d == 0) {
                        xi(0) = u(0) + (+1. - u(0)) * v(0);
                        eta(0) = u(0) + (-1. - u(0)) * v(0);
                    }
                    else {
                        xi(0) = u(0) + (-1. - u(0)) * v(0);
                        eta(0) = u(0) + (+1. - u(0)) * v(0);
                    }
 
                    double Jxi = elem.get_normal(xi).norm();
                    std::complex<double> G = kernel(elem.get_x(xi), elem.get_x(eta), elem.get_normal(eta));
 
                    auto Nxi = test_shape_t::template eval_shape<0>(xi);
                    auto Neta = trial_shape_t::template eval_shape<0>(eta);
 
                    // evaluate integrand (Jeta is incorporated in the Kernel with the normal)
                    result_t F = G * (Jxi * 2 * (1. - v(0))) * (Nxi * Neta.transpose());
 
                    I += F * wv * wu;
                } // loop over domains
            } // loop over u
        } // loop over v
        return I;
    } // end fo function eval
}; // class helmholtz_2d_DLP_galerkin_face_general
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_DLPt_galerkin_face_general
{
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef typename test_field_t::nset_t test_shape_t;
    typedef typename trial_field_t::nset_t trial_shape_t;
 
    static size_t const nTest = test_shape_t::num_nodes;
    static size_t const nTrial = trial_shape_t::num_nodes;
 
    typedef Eigen::Matrix<std::complex<double>, nTest, nTrial> result_t;
 
    typedef typename trial_field_t::elem_t elem_t;
    typedef typename elem_t::domain_t domain_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename elem_t::x_t x_t;
 
    typedef regular_quad_store<domain_t, order> quadrature_t;
 
public:
    template <class wave_number_t>
    static result_t eval(elem_t const &elem, wave_number_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_DLPt_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        typedef helmholtz_2d_DLP_kernel<wave_number_t> kernel_t;
        kernel_t kernel(k);
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto itv = quadrature_t::quadrature.begin(); itv != quadrature_t::quadrature.end(); ++itv) {
            xi_t v = itv->get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = itv->get_w() / 2.;
 
            // loop over u
            for (auto itu = quadrature_t::quadrature.begin(); itu != quadrature_t::quadrature.end(); ++itu) {
                xi_t u = itu->get_xi();
                double wu = itu->get_w();
 
                // loop over domains
                for (int d = 0; d < 2; ++d) {
                    xi_t xi, eta;
 
                    if (d == 0) {
                        xi(0) = u(0) + (+1. - u(0)) * v(0);
                        eta(0) = u(0) + (-1. - u(0)) * v(0);
                    }
                    else {
                        xi(0) = u(0) + (-1. - u(0)) * v(0);
                        eta(0) = u(0) + (+1. - u(0)) * v(0);
                    }
 
                    double Jeta = elem.get_normal(eta).norm();
                    std::complex<double> G = kernel(elem.get_x(xi), elem.get_x(eta), elem.get_normal(xi));
 
                    auto Nxi = test_shape_t::template eval_shape<0>(xi);
                    auto Neta = trial_shape_t::template eval_shape<0>(eta);
 
                    // evaluate integrand (Jeta is incorporated in the Kernel with the normal)
                    result_t F = G * (Jeta * 2 * (1. - v(0))) * (Nxi * Neta.transpose());
 
                    I += F * wv * wu;
                } // loop over domains
            } // loop over u
        } // loop over v
        return I;
    } // end fo function eval
}; // class helmholtz_2d_DLPt_galerkin_face_general
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_HSP_galerkin_face_general
{
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef typename test_field_t::nset_t test_shape_t;
    typedef typename trial_field_t::nset_t trial_shape_t;
 
    static size_t const nTest = test_shape_t::num_nodes;
    static size_t const nTrial = trial_shape_t::num_nodes;
 
    typedef Eigen::Matrix<std::complex<double>, nTest, nTrial> result_t;
 
    typedef typename trial_field_t::elem_t elem_t;
    typedef typename elem_t::domain_t domain_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename elem_t::x_t x_t;
 
    typedef regular_quad_store<domain_t, order> quadrature_t;
 
    static result_t F2(xi_t const &u)
    {
        using namespace boost::math::double_constants;
        return test_shape_t::template eval_shape<0>(u)
            * trial_shape_t::template eval_shape<0>(u).transpose() / (2.0 * two_pi);
    }
 
public:
    template <class wave_number_t>
    static result_t eval(elem_t const &elem, wave_number_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_HSP_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        typedef helmholtz_2d_HSP_kernel<wave_number_t> kernel_t;
        kernel_t kernel(k);
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto it = quadrature_t::quadrature.begin(); it != quadrature_t::quadrature.end(); ++it) {
            xi_t v = it->get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = it->get_w() / 2.;
 
            // loop over u
            for (auto jt = quadrature_t::quadrature.begin(); jt != quadrature_t::quadrature.end(); ++jt) {
                xi_t u = jt->get_xi();
                double wu = jt->get_w();
 
                // loop over domains
                for (int d = 0; d < 2; ++d) {
                    xi_t xi, eta;
 
                    if (d == 0) {
                        xi(0) = u(0) + (+1. - u(0)) * v(0);
                        eta(0) = u(0) + (-1. - u(0)) * v(0);
                    }
                    else {
                        xi(0) = u(0) + (-1. - u(0)) * v(0);
                        eta(0) = u(0) + (+1. - u(0)) * v(0);
                    }
 
                    x_t x = elem.get_x(xi);
                    x_t y = elem.get_x(eta);
                    x_t Jeta = elem.get_normal(eta);
                    x_t Jxi = elem.get_normal(xi);
                    std::complex<double> G = kernel(x, y, Jxi, Jeta);
 
                    auto Nxi = test_shape_t::template eval_shape<0>(xi);
                    auto Neta = trial_shape_t::template eval_shape<0>(eta);
 
                    // evaluate integrand
                    result_t F = G * 2. * (1. - v(0)) * (Nxi * Neta.transpose());
 
                    I += F * wv * wu;
                } // loop over domains
 
                I -= 2.0 * F2(u) / (v(0) * v(0)) * wv * wu;
            } // loop over u
 
            I += 2.0 * (F2(xi_t(1.)) + F2(xi_t(-1.0))) / v(0) * wv;
        } // loop over v
 
        for (auto jt = quadrature_t::quadrature.begin(); jt != quadrature_t::quadrature.end(); ++jt) {
            xi_t u = jt->get_xi();
            double wu = jt->get_w();
            I -= 2.0 * F2(u) * wu;
        }
 
        I -= 2.0 * (F2(xi_t(1.0)) * std::log(2.0 * elem.get_dx(xi_t(1.0)).norm())
            + F2(xi_t(-1.0)) * std::log(2.0 * elem.get_dx(xi_t(-1.0)).norm()));
 
        return I;
    } // end fo function eval
}; // class helmholtz_2d_HSP_galerkin_face_general
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_HSP_galerkin_edge_general
{
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef typename test_field_t::nset_t test_shape_t;
    typedef typename test_shape_t::shape_t test_N_t;
    typedef typename trial_field_t::nset_t trial_shape_t;
    typedef typename trial_shape_t::shape_t trial_N_t;
 
    static size_t const nTest = test_shape_t::num_nodes;
    static size_t const nTrial = trial_shape_t::num_nodes;
 
    typedef Eigen::Matrix<std::complex<double> , nTest, nTrial> result_t;
 
    typedef typename trial_field_t::elem_t trial_elem_t;
    typedef typename test_field_t::elem_t test_elem_t;
 
    typedef typename trial_elem_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef typename trial_elem_t::x_t x_t;
 
    typedef regular_quad_store<domain_t, order> quadrature_t;
 
public:
    template<class wave_number_t>
    static result_t eval(
        test_elem_t const &test_elem,
        trial_elem_t const &trial_elem,
        wave_number_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_HSP_galerkin_edge_general" << std::endl;
            printed = true;
        }
#endif
        using namespace boost::math::double_constants;
 
        typedef helmholtz_2d_HSP_kernel<wave_number_t> kernel_t;
        kernel_t kernel(k);
 
        result_t I = result_t::Zero();
 
        int singular_domain;
        if (test_elem.get_nodes()(test_elem_t::num_nodes - 1) == trial_elem.get_nodes()(0))
            singular_domain = 1;
        else if (test_elem.get_nodes()(0) == trial_elem.get_nodes()(trial_elem_t::num_nodes - 1))
            singular_domain = 0;
        else throw std::logic_error("Invalid singular case detected");
 
        //  loop over two subdomains
        for (int d = 0; d < 2; ++d) {
            xi_t xi0, eta0;
            if (d == 0) {
                xi0 << -1;
                eta0 << +1;
            }
            else {
                xi0 << +1;
                eta0 << -1;
            }
 
            bool singular = singular_domain == d;
 
            x_t ax0, ay0, Jx0, Jy0;
            test_N_t Nx0;
            trial_N_t Ny0;
 
            // prepare quantities needed to compute F1(u)
            if (singular) {
                ax0 = test_elem.get_dx(xi0);
                ay0 = trial_elem.get_dx(eta0);
                Jx0 = test_elem.get_normal(xi0);
                Jy0 = trial_elem.get_normal(eta0);
                Nx0 = test_shape_t::template eval_shape<0>(xi0);
                Ny0 = trial_shape_t::template eval_shape<0>(eta0);
            }
 
            // loop over u
            for (auto itu = quadrature_t::quadrature.begin(); itu != quadrature_t::quadrature.end(); ++itu) {
                xi_t u = itu->get_xi();
                double wu = itu->get_w();
 
                // s and F1 are only computed if singular
                double s;
                result_t F1;
 
                if (singular) {
                    x_t svec = ay0 * (u - eta0)(0) - ax0 * (u - xi0)(0);
                    s = svec.norm();
                    x_t es = svec.normalized();
                    F1 = (Jx0.dot(Jy0) - 2 * es.dot(Jx0) * es.dot(Jy0)) / pi / s / s *
                        (Nx0 * Ny0.transpose());
                }
 
                // loop over v
                for (auto itv = quadrature_t::quadrature.begin(); itv != quadrature_t::quadrature.end(); ++itv) {
                    xi_t v = itv->get_xi();
                    // transform v to (0,1) from (-1,1)
                    v(0) = (v(0) + 1.) / 2.;
                    double wv = itv->get_w() / 2.;
 
                    // compute intrinsic variables
                    xi_t xi = xi0 + (u - xi0) * v;
                    xi_t eta = eta0 + (u - eta0) * v;
 
                    // compute kernel value
                    x_t x = test_elem.get_x(xi);
                    x_t y = trial_elem.get_x(eta);
                    x_t Jxi = test_elem.get_normal(xi);
                    x_t Jeta = trial_elem.get_normal(eta);
                    std::complex<double> G = kernel(x, y, Jxi, Jeta);
 
                    test_N_t Nxi = test_shape_t::template eval_shape<0>(xi);
                    trial_N_t Neta = trial_shape_t::template eval_shape<0>(eta);
 
                    // evaluate integrand
                    result_t F = G * 2. * v(0) * (Nxi * Neta.transpose());
 
                    I += F * wv * wu;
 
                    if (singular)
                        // subtract singular part
                        I -= F1 / v(0) * wv * wu;
                } // loop over v
 
                if (singular)
                    // add analytic (inner) integral of singular part
                    I += F1 * std::log(s) * wu;
 
            } // loop over u
        } // loop over domains
 
        return I;
    } // function eval
}; // class helmholtz_2d_HSP_galerkin_edge_general
 
 
template <class WaveNumber, class TestField, class TrialField>
class singular_integral_shortcut<
    helmholtz_2d_SLP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    !(is_constant_line<TrialField>::value && is_constant_line<TestField>::value)
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_SLP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_SLP_galerkin_face_general<TestField, TrialField, 10>::eval(
            trial_field.get_elem(), kernel.derived().get_wave_number());
        return result;
    }
};
 
 
template <class WaveNumber, class TestField, class TrialField>
class singular_integral_shortcut<
    helmholtz_2d_SLP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    is_constant_line<TrialField>::value &&
    is_constant_line<TestField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_SLP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        auto const &elem = trial_field.get_elem();
        double R = (elem.get_coords().col(1) - elem.get_coords().col(0)).norm()/2.;
        result(0, 0) = helmholtz_2d_SLP_galerkin_face_constant_line<5>::eval(
            R, kernel.derived().get_wave_number());
 
        return result;
    }
};
 
 
template <class TestField, class TrialField, class WaveNumber>
class singular_integral_shortcut<
    helmholtz_2d_DLP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    !std::is_same<typename TrialField::elem_t::lset_t, line_1_shape_set>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_DLP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_DLP_galerkin_face_general<TestField, TrialField, 11>::eval(
            test_field.get_elem(), kernel.derived().get_wave_number());
        return result;
    }
};
 
 
template <class TestField, class TrialField, class WaveNumber>
class singular_integral_shortcut<
    helmholtz_2d_DLPt_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    !std::is_same<typename TrialField::elem_t::lset_t, line_1_shape_set>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_DLPt_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_DLPt_galerkin_face_general<TestField, TrialField, 11>::eval(
            test_field.get_elem(), kernel.derived().get_wave_number());
        return result;
    }
};
 
 
template <class TestField, class TrialField, class WaveNumber>
class singular_integral_shortcut<
    helmholtz_2d_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_HSP_galerkin_face_general<TestField, TrialField, 11>::eval(
            test_field.get_elem(), kernel.derived().get_wave_number());
        return result;
    }
};
 
 
template <class TestField, class TrialField, class WaveNumber>
class singular_integral_shortcut<
    helmholtz_2d_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_0d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<helmholtz_2d_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_HSP_galerkin_edge_general<TestField, TrialField, 11>::eval(
            test_field.get_elem(),
            trial_field.get_elem(),
            kernel.derived().get_wave_number());
        return result;
    }
};
} // end of namespace NiHu
 
 
#endif // NIHU_HELMHOLTZ_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED