helmholtz_2d_singular_collocation_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_2D_SINGULAR_COLLOCATION_INTEGRALS_HPP_INCLUDED
#define NIHU_HELMHOLTZ_2D_SINGULAR_COLLOCATION_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_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_collocation_curved
{
    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> reg_quadrature_t;
    typedef log_quad_store<order> log_quadrature_t;
 
    
public:
    template <class WaveNumber>
    static result_t eval(elem_t const &elem, WaveNumber const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_SLP_collocation_curved" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        helmholtz_2d_SLP_kernel<WaveNumber> kernel(k);
 
        result_t result = result_t::Zero();
 
        // traverse collocation points (rows)
        for (size_t i = 0; i < nTest; ++i) {
            xi_t const &xi0 = test_shape_t::corner_at(i);
 
            // the singular location
            x_t x = elem.get_x(xi0);
 
            // traverse integration domains
            for (int dom = 0; dom < 2; ++dom) {
                xi_t xi_lim = domain_t::get_corner(dom);
                double jac = std::abs(xi_lim(0) - xi0(0));
 
                // traverse regular quadrature points
                for (auto const &q : reg_quadrature_t::quadrature) {
                    // transform quadrature to [0 1];
                    double v = (q.get_xi()(0) + 1.) / 2.;
                    double wv = q.get_w() / 2;
 
                    xi_t eta = xi0 + v * (xi_lim - xi0);
                    double r = (elem.get_x(eta) - x).norm();
                    double Jeta = elem.get_dx(eta).norm();
 
                    // evaluate Green's function
                    std::complex<double> Gtot;
                    kernel.template eval<0>(r, &Gtot);
                    double Glog = -std::log(v) / two_pi;
                    std::complex<double> Greg = Gtot - Glog;
 
                    // integrate difference numerically
                    result.row(i) += Greg * Jeta * jac * wv
                        * trial_shape_t::template eval_shape<0>(eta);
                } // loop over regular quadrature
 
                // traverse log quadrature points
                for (auto const &q : log_quadrature_t::quadrature) {
                    double v = q.get_xi()(0);
                    double wv = q.get_w();
 
                    xi_t eta = xi0 + v * (xi_lim - xi0);
                    double Jeta = elem.get_dx(eta).norm();
 
                    // evaluate Green's function
                    double Glog = 1.0 / two_pi;
 
                    // integrate difference numerically
                    result.row(i) += Glog * Jeta * jac * wv
                        * trial_shape_t::template eval_shape<0>(eta);
                } // loop over log quadrature
            } // loop over subdomains
        } // loop over collocation points
 
        return result;
    } // function eval
}; //class helmholtz_2d_SLP_collocation_curved
 
 
template <unsigned expansion_length>
class helmholtz_2d_SLP_collocation_constant_line
{
public:
    template <class wavenumber_t>
    static std::complex<double> eval(
        line_1_elem const &elem,
        line_1_elem::x_t const &x0,
        wavenumber_t const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_SLP_collocation_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
        
        std::complex<double> const c(euler, half_pi);
 
        // compute elem radius
        auto R = (x0 - elem.get_coords().col(0)).norm();
 
        auto Q = k * R / 2.;
        auto clnq = c + std::log(Q);
 
        auto res(clnq - 1.);    // initial (k=0) value of the result
        decltype(Q) B(1.);      // the power term in the series
        double Cn = 0.0;        // the harmonic series up to n = 0
        for (unsigned n = 1; n <= expansion_length; ++n) {
            B *= -Q*Q / n / n;      // the actual power term
            Cn += 1. / n;       // the actual harmonic term
            res += B / (2 * n + 1) * (clnq - Cn - 1. / (2 * n + 1));
        }
 
        return -R / pi * res;
    } // function eval
}; // helmholtz_2d_SLP_collocation_constant_line
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_DLP_collocation_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> quadr_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_collocation_general" << std::endl;
            printed = true;
        }
#endif
 
        result_t result = result_t::Zero();
        
        helmholtz_2d_DLP_kernel<wave_number_t> kernel(k);
        
        // traverse collocation points
        for (size_t i = 0; i < nTest; ++i) {
            xi_t xi0 = test_shape_t::corner_at(i);
            
            // singular point and normal at singular point
            x_t x = elem.get_x(xi0);
 
            for (size_t dom = 0; dom < 2; ++dom) {
                xi_t const &a = domain_t::get_corner(dom);
            
                // traverse quadrature points
                for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it) {
                    // transform quadrature to [a xi0]
                    xi_t xi = it->get_xi() * ((xi0(0)-a(0))/2.) + (xi0+a)/2.;
                    double w = it->get_w() * std::abs(xi0(0)-a(0))/2.;
                    double rho = xi(0) - xi0(0);
                    
                    // get trial location, jacobian and normal
                    x_t y = elem.get_x(xi);
                    x_t Jyvec = elem.get_normal(xi);
                    
                    // evaluate Green's function
                    std::complex<double> G = kernel(x, y, Jyvec);
                    // multiply by Shape function and Jacobian
                    auto N = trial_shape_t::template eval_shape<0>(xi);
                    auto F = G * N;
                    
                    // integrate difference numerically
                    result.row(i) += F * w;
                } // quadrature points
            } // domains
        } // collocation points
        return result;
    } // function eval
}; // class helmholtz_2d_DLP_collocation_general
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_DLPt_collocation_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> quadr_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_collocation_general" << std::endl;
            printed = true;
        }
#endif
 
        result_t result = result_t::Zero();
        
        helmholtz_2d_DLPt_kernel<wave_number_t> kernel(k);
        
        // traverse collocation points
        for (size_t i = 0; i < nTest; ++i)
        {
            xi_t xi0 = test_shape_t::corner_at(i);
            
            // singular point and normal at singular point
            x_t x = elem.get_x(xi0);
            x_t Jxvec = elem.get_normal(xi0);
            x_t nx = Jxvec.normalized();
 
            for (size_t dom = 0; dom < 2; ++dom) {
                xi_t const &a = domain_t::get_corner(dom);
            
                // traverse quadrature points
                for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it) {
                    // transform quadrature to [a xi0]
                    xi_t xi = it->get_xi() * ((xi0(0)-a(0))/2.) + (xi0+a)/2.;
                    double w = it->get_w() * std::abs(xi0(0)-a(0)) / 2.;
                    double rho = xi(0) - xi0(0);
                    
                    // get trial location, jacobian and normal
                    x_t y = elem.get_x(xi);
                    x_t Jyvec = elem.get_normal(xi);
                    double jac = Jyvec.norm();
                    
                    // evaluate Green's function
                    std::complex<double> G = kernel(x, y, nx);
                    // multiply by Shape function and Jacobian
                    auto N = trial_shape_t::template eval_shape<0>(xi);
                    auto F = G * jac * N;
                    
                    // integrate numerically
                    result.row(i) += F * w;
                } // quadrature points
            } // domains
        } // collocation points
        return result;
    } // function eval
}; // class helmholtz_2d_DLPt_collocation_general
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_HSP_collocation_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> quadr_t;
    
public:
    
    template <class WaveNumber>
    static result_t eval(elem_t const &elem, WaveNumber  const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_HSP_collocation_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        // instantiate and clear result matrix
        result_t result = result_t::Zero();
        
        // instantiate the kernel
        helmholtz_2d_HSP_kernel<WaveNumber> kernel(k);
        
        // integration limits in the reference domain
        xi_t const &a = domain_t::get_corner(0);
        xi_t const &b = domain_t::get_corner(1);
        
        // traverse collocation points
        for (size_t i = 0; i < nTest; ++i)
        {
            // the collocation point
            xi_t const &xi0 = test_shape_t::corner_at(i);
            
            // trial shape function and its derivative at the singular point
            auto N0 = trial_shape_t::template eval_shape<0>(xi0);
            auto N1 = trial_shape_t::template eval_shape<1>(xi0);
            
            // singular point and normal at singular point in physical domain
            x_t x = elem.get_x(xi0);
            x_t Jxvec = elem.get_normal(xi0);
            x_t nx = Jxvec.normalized();
            double jac0 = Jxvec.norm();
            
            // traverse quadrature points
            for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
            {
                // transform quadrature to [a xi0]
                xi_t xi = it->get_xi() * ((xi0(0)-a(0))/2.) + (xi0+a)/2.;
                double w = it->get_w() * ((xi0(0)-a(0))/2.);
                double rho = xi(0) - xi0(0);
                
                // get trial location, jacobian and normal
                x_t y = elem.get_x(xi);
                x_t Jyvec = elem.get_normal(xi);
                x_t ny = Jyvec.normalized();
                double jac = Jyvec.norm();
                
                // evaluate Green's function
                std::complex<double> G = kernel(x, y, nx, ny);
                // multiply by Shape function and Jacobian
                auto N = trial_shape_t::template eval_shape<0>(xi);
                auto F = (G * jac) * N;
                // evaluate singular part
                auto F0 = (
                    (N0/rho + N1)/rho/jac0
                    - k*k/2. * N0 * jac0 * std::log(std::abs(rho) * jac0)
                ) / two_pi;
                
                // integrate difference numerically
                result.row(i) += (F - F0) * w;
            }
            
            // traverse quadrature points
            for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
            {
                // transform quadrature to [xi0 b];
                xi_t xi = it->get_xi() * ((b(0)-xi0(0))/2.) + (xi0+b)/2.;
                double w = it->get_w() * ((b(0)-xi0(0))/2.);
                double rho = xi(0) - xi0(0);
                
                // get trial location, jacobian and normal
                x_t y = elem.get_x(xi);
                x_t Jyvec = elem.get_normal(xi);
                x_t ny = Jyvec.normalized();
                double jac = Jyvec.norm();
                
                // evaluate Green's function
                std::complex<double> G = kernel(x, y, nx, ny);
                // multiply by Shape function and Jacobian
                auto N = trial_shape_t::template eval_shape<0>(xi);
                auto F = (G * jac) * N;
                // evaluate singular part
                auto F0 = (
                    (N0/rho + N1)/rho/jac0
                    - k*k/2. * N0 * jac0 * std::log(std::abs(rho) * jac0)
                ) / two_pi;
                
                // integrate difference numerically
                result.row(i) += (F - F0) * w;
            }
            
            // add analytic integral of singular part
            double d1 = std::abs(a(0) - xi0(0)) * jac0;
            double d2 = std::abs(b(0) - xi0(0)) * jac0;
            result.row(i) += (
                ((1./(a(0)-xi0(0)) - 1./(b(0)-xi0(0))) * N0 +
                std::log(std::abs((b(0)-xi0(0)) / (a(0)-xi0(0)))) * N1) / jac0
                - k*k/2. *
                N0 * (d1 * (std::log(d1) - 1.) + d2 * (std::log(d2) - 1.))
                ) / two_pi;
        }
        return result;
    }
};
 
 
template <class TestField, class TrialField, size_t order>
class helmholtz_2d_HSP_collocation_straight_line
{
    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 NiHu::line_1_elem 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> quadr_t;
    
public:
    
    template <class WaveNumber>
    static result_t eval(elem_t const &elem, WaveNumber  const &k)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using helmholtz_2d_HSP_collocation_straight_line" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        // instantiate and clear result matrix
        result_t result = result_t::Zero();
        
        // instantiate the kernel
        helmholtz_2d_HSP_kernel<WaveNumber> kernel(k);
        
        // integration limits in the reference domain
        xi_t const &a = domain_t::get_corner(0);
        xi_t const &b = domain_t::get_corner(1);
        
        // traverse collocation points
        for (size_t i = 0; i < nTest; ++i)
        {
            // the collocation point
            xi_t const &xi0 = test_shape_t::corner_at(i);
            
            // trial shape function and its derivative at the singular point
            auto N0 = trial_shape_t::template eval_shape<0>(xi0);
            auto N1 = trial_shape_t::template eval_shape<1>(xi0);
            
            // singular point and normal in physical domain
            x_t x = elem.get_x(xi0);
            x_t Jvec = elem.get_normal(xi0);
            x_t n = Jvec.normalized();
            double jac = Jvec.norm();
            
            // traverse quadrature points
            for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
            {
                // transform quadrature to [a xi0]
                xi_t xi = it->get_xi() * ((xi0(0)-a(0))/2.) + (xi0+a)/2.;
                double w = it->get_w() * ((xi0(0)-a(0))/2.);
                double rho = xi(0) - xi0(0);
                
                // get trial location, jacobian and normal
                x_t y = elem.get_x(xi);
                
                // evaluate Green's function
                std::complex<double> G = kernel(x, y, n, n);
                // multiply by Shape function and Jacobian
                auto N = trial_shape_t::template eval_shape<0>(xi);
                auto F = (G * jac) * N;
                // evaluate singular part
                auto F0 = (
                    (N0/rho + N1)/rho/jac
                    - k*k/2. * N0 * jac * std::log(std::abs(rho) * jac)
                ) / two_pi;
                
                // integrate difference numerically
                result.row(i) += (F - F0) * w;
            }
            
            // traverse quadrature points
            for (auto it = quadr_t::quadrature.begin(); it != quadr_t::quadrature.end(); ++it)
            {
                // transform quadrature to [xi0 b];
                xi_t xi = it->get_xi() * ((b(0)-xi0(0))/2.) + (xi0+b)/2.;
                double w = it->get_w() * ((b(0)-xi0(0))/2.);
                double rho = xi(0) - xi0(0);
                
                // get trial location, jacobian and normal
                x_t y = elem.get_x(xi);
                
                // evaluate Green's function
                std::complex<double> G = kernel(x, y, n, n);
                // multiply by Shape function and Jacobian
                auto N = trial_shape_t::template eval_shape<0>(xi);
                auto F = (G * jac) * N;
                // evaluate singular part
                auto F0 = (
                    (N0/rho + N1)/rho/jac
                    - k*k/2. * N0 * jac * std::log(std::abs(rho) * jac)
                ) / two_pi;
                
                // integrate difference numerically
                result.row(i) += (F - F0) * w;
            }
            
            // add analytic integral of singular part
            double d1 = std::abs(a(0) - xi0(0)) * jac;
            double d2 = std::abs(b(0) - xi0(0)) * jac;
            result.row(i) += (
                ((1./(a(0)-xi0(0)) - 1./(b(0)-xi0(0))) * N0 +
                std::log(std::abs((b(0)-xi0(0)) / (a(0)-xi0(0)))) * N1) / jac
                - k*k/2. *
                N0 * (d1 * (std::log(d1) - 1.) + d2 * (std::log(d2) - 1.))
                ) / two_pi;
        }
        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<
    is_collocational<TestField, TrialField>::value &&
    is_constant_line<TrialField>::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(0, 0) = helmholtz_2d_SLP_collocation_constant_line<9>::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_2d_SLP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::value &&
    !is_constant_line<TrialField>::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_collocation_curved<
            TestField, TrialField, 30
        >::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_DLP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::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 &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_DLP_collocation_general<TestField, TrialField, 30>::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_DLPt_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::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 &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_DLPt_collocation_general<TestField, TrialField, 30>::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_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::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_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_HSP_collocation_straight_line<TestField, TrialField, 30>::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_HSP_kernel<WaveNumber>, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::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_HSP_kernel<WaveNumber> > const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = helmholtz_2d_HSP_collocation_general<TestField, TrialField, 30>::eval(
            trial_field.get_elem(),
            kernel.derived().get_wave_number());
            
        return result;
    }
};
 
} // end of namespace NiHu
 
#endif // NIHU_HELMHOLTZ_2D_SINGULAR_COLLOCATION_INTEGRALS_HPP_INCLUDED