laplace_2d_singular_double_integrals.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 LAPLACE_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED
#define LAPLACE_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED
 
#include "../core/match_types.hpp"
#include "../core/singular_integral_shortcut.hpp"
#include "../util/math_functions.hpp"
#include "field_type_helpers.hpp"
#include "lib_element.hpp"
#include "laplace_kernel.hpp"
#include "normal_derivative_singular_integrals.hpp"
#include "quadrature_store_helper.hpp"
 
#include <boost/math/constants/constants.hpp>
 
#define NIHU_PRINT_DEBUG
 
#ifdef NIHU_PRINT_DEBUG
#include <iostream>
#endif
 
namespace NiHu
{
 
template <class TestField, class TrialField, size_t order>
class laplace_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<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 size_t const inner_order =
        elem_t::lset_t::jacobian_order * 2 +        // Jy(v) (possible overkill)
        elem_t::lset_t::jacobian_order * 2 +        // Jx(v) (possible overkill)
        trial_field_t::nset_t::polynomial_order +   // Ny(v) exact
        test_field_t::nset_t::polynomial_order +    // Nx(v) exact
        1;                                          // 1-v (exact)
    typedef log_quad_store<inner_order> log_quadrature_t;
 
public:
    static result_t eval(elem_t const &elem)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_SLP_galerkin_face_general, inner order: " << inner_order << std::endl;
            printed = true;
        }
#endif
 
        result_t I = result_t::Zero();
 
        // loop over u
        for (auto const &qu : quadrature_t::quadrature) {
            double u = qu.get_xi()(0);
            double wu = qu.get_w();
 
            // loop over Gauss-Legendre base points
            for (auto const &qv : quadrature_t::quadrature) {
                double v = (qv.get_xi()(0) + 1.) / 2.;
                double wv = qv.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 + (+1. - u) * v;
                        eta(0) = u + (-1. - u) * v;
                    }
                    else {
                        xi(0) = u + (-1. - u) * v;
                        eta(0) = u + (+1. - u) * v;
                    }
 
                    x_t x = elem.get_x(xi);
                    x_t y = elem.get_x(eta);
                    double r = (y - x).norm();
                    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);
 
                    // The regular part of the kernel
                    double Greg = -std::log(r / v) / two_pi;
                    // evaluate regular integrand
                    I += Greg * Jxi * Jeta * (2. * (1. - v)) * wu * wv
                        * (Nxi * Neta.transpose());
                } // loop over domains
            } // loop over v
 
            // loop over logarithmic Gaussian base points
            for (auto const &quad : log_quadrature_t::quadrature) {
                double v = quad.get_xi()(0);
                double wv = quad.get_w();
 
                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 + (+1. - u) * v;
                        eta(0) = u + (-1. - u) * v;
                    }
                    else {
                        xi(0) = u + (-1. - u) * v;
                        eta(0) = u + (+1. - u) * v;
                    }
 
                    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);
 
                    double Gsing = 1.0 / two_pi; //  -log(v) / two_pi;
 
                    // evaluate regular integrand
                    result_t F = Gsing * Jxi * Jeta * (Nxi * Neta.transpose());
                    I += F * (2. * (1. - v)) * wu * wv;
                } // loop over domains
            } // loop over v
        } // loop over u
        return I;
    } // function eval
}; // class laplace_2d_SLP_galerkin_face_general
 
 
class laplace_2d_SLP_galerkin_face_constant_line
{
public:
    static double eval(line_1_elem const &elem)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_SLP_galerkin_face_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        auto const &C = elem.get_coords();
        double d = (C.col(1) - C.col(0)).norm();
        return d * d * (1.5 - std::log(d)) / two_pi;
    }
};
 
class laplace_2d_SLP_galerkin_face_linear_line
{
public:
    static void eval(line_1_elem const &elem, double &i1, double &i2)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_SLP_galerkin_face_linear_line" << std::endl;
            printed = true;
        }
#endif
 
        auto const &C = elem.get_coords();
        double d = (C.col(1) - C.col(0)).norm();    // element length
        double c = d * d / (8. * pi);
        double lnd = std::log(d);
        i1 = c * (1.75 - lnd);
        i2 = c * (1.25 - lnd);
    }
};
 
 
template <class TestField, class TrialField, size_t order>
class laplace_2d_SLP_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<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> reg_quadrature_t;
    static size_t const inner_order =
        trial_elem_t::lset_t::jacobian_order * 2 +  // Jy(v) (possible overkill)
        test_elem_t::lset_t::jacobian_order * 2 +   // Jx(v) (possible overkill)
        trial_field_t::nset_t::polynomial_order +   // Ny(v) exact
        test_field_t::nset_t::polynomial_order +    // Nx(v) exact
        1;                                          // v (exact)
    typedef log_quad_store<inner_order> log_quadrature_t;
 
public:
    static result_t eval(
        test_elem_t const &test_elem,
        trial_elem_t const &trial_elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_SLP_galerkin_edge_general, inner order: " << inner_order << std::endl;
            printed = true;
        }
#endif
        using namespace boost::math::double_constants;
 
        laplace_2d_SLP_kernel kernel;
 
        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;
 
            // loop over u
            for (auto const &qu : reg_quadrature_t::quadrature) {
                xi_t u = qu.get_xi();
                double wu = qu.get_w();
 
                // loop over v
                for (auto const &qv : reg_quadrature_t::quadrature) {
                    double v = (qv.get_xi()(0) + 1.) / 2.;
                    double wv = qv.get_w() / 2.;
 
                    // compute intrinsic variables
                    xi_t xi = xi0 + (u - xi0) * v;
                    xi_t eta = eta0 + (u - eta0) * v;
 
                    // compute kernel value
                    double r = (trial_elem.get_x(eta) - test_elem.get_x(xi)).norm();
                    double G = -std::log(r) / two_pi;
                    if (singular)
                        G -= -std::log(v) / two_pi;
 
                    double Jxi = test_elem.get_normal(xi).norm();
                    double Jeta = trial_elem.get_normal(eta).norm();
                    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
                    I += G * Jxi * Jeta * (Nxi * Neta.transpose()) * 2. * v * wv * wu;
                } // loop over regular v
 
                if (!singular)
                    continue;
 
                // loop over logarithmic v
                for (auto const &qv : log_quadrature_t::quadrature) {
                    double v = qv.get_xi()(0);
                    double wv = qv.get_w();
 
                    // compute intrinsic variables
                    xi_t xi = xi0 + (u - xi0) * v;
                    xi_t eta = eta0 + (u - eta0) * v;
 
                    // compute kernel value
                    double G = 1.0 / two_pi;
                    double Jxi = test_elem.get_normal(xi).norm();
                    double Jeta = trial_elem.get_normal(eta).norm();
                    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
                    I += G * Jxi * Jeta * (Nxi * Neta.transpose()) * 2. * v * wv * wu;
                } // loop over singular v
            } // loop over u
        } // loop over domains
 
        return I;
    } // function eval
}; // class laplace_2d_SLP_galerkin_edge_general
 
 
class laplace_2d_SLP_galerkin_edge_constant_line
{
    static double qfunc(double a, double phi)
    {
        using namespace boost::math::double_constants;
 
        if (std::abs(phi) < 1e-3)
            return a / (a + 1.);
        double cotphi = std::tan(half_pi - phi);
        return std::atan(a / std::sin(phi) + cotphi) - std::atan(cotphi);
    }
 
public:
    static double eval(line_1_elem const &elem1, line_1_elem const &elem2)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_SLP_galerkin_edge_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        // get the element corner coordinates
        auto const &C1 = elem1.get_coords();
        auto const &C2 = elem2.get_coords();
        // and side vectors
        auto r1vec = (C1.col(1) - C1.col(0)).eval(), r2vec = (C2.col(1) - C2.col(0)).eval();
        // and side lengths
        double r1 = r1vec.norm(), r2 = r2vec.norm();
        // and signed angle between them
        double phi = std::asin(r1vec(0) * r2vec(1) - r2vec(0) * r1vec(1)) / (r1 * r2);
        // third side length
        double r3 = std::sqrt(r1 * r1 + 2 * r1 * r2 * std::cos(phi) + r2 * r2);
        return (
            r1 * r2 * (3. - 2. * std::log(r3))
            + std::cos(phi) * (r1 * r1 * std::log(r1 / r3) + r2 * r2 * std::log(r2 / r3))
            - std::sin(phi) * (r1 * r1 * qfunc(r2 / r1, phi) + r2 * r2 * qfunc(r1 / r2, phi))
            ) / (4. * pi);
    }
};
 
 
template <class TestField, class TrialField, size_t order>
class laplace_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<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;
 
    typedef laplace_2d_DLP_kernel kernel_t;
 
public:
    static result_t eval(elem_t const &elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_DLP_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        kernel_t kernel;
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto const &qv : quadrature_t::quadrature) {
            xi_t v = qv.get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = qv.get_w() / 2.;
 
            // loop over u
            for (auto const &qu : quadrature_t::quadrature) {
                xi_t u = qu.get_xi();
                double wu = qu.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();
                    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 laplace_2d_DLP_galerkin_face_general
 
 
class laplace_2d_DLP_galerkin_edge_constant_line
{
    static double qfunc(double a, double phi)
    {
        using namespace boost::math::double_constants;
 
        if (std::abs(phi) < 1e-3)
            return a / (a + 1.);
        double cotphi = std::tan(half_pi - phi);
        return std::atan(a / std::sin(phi) + cotphi) - std::atan(cotphi);
    }
 
public:
    static double eval(line_1_elem const &elem1, line_1_elem const &elem2)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_DLP_galerkin_edge_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        // get element corners
        auto const &C1 = elem1.get_coords();
        auto const &C2 = elem2.get_coords();
        // get side vectors
        auto r1vec = (C1.col(1) - C1.col(0)).eval(), r2vec = (C2.col(1) - C2.col(0)).eval();
        // get side lengths
        double r1 = r1vec.norm(), r2 = r2vec.norm();
        // get angle between elements
        double phi = std::asin(r1vec(0) * r2vec(1) - r2vec(0) * r1vec(1)) / (r1 * r2);
        // general expression
        double r3 = std::sqrt(r1 * r1 + 2 * r1 * r2 * std::cos(phi) + r2 * r2);
        double res = (
            r2 * std::cos(phi) * qfunc(r1 / r2, phi)
            - r1 * qfunc(r2 / r1, phi)
            + r2 * std::sin(phi) * std::log(r2 / r3)
            ) / two_pi;
 
        if (elem1.get_nodes()(0) == elem2.get_nodes()(1))
            res *= -1;
 
        return res;
    }
};
 
 
template <class TestField, class TrialField, size_t order>
class laplace_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<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;
 
    typedef laplace_2d_DLPt_kernel kernel_t;
 
public:
    static result_t eval(elem_t const &elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_DLPt_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        kernel_t kernel;
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto const &qv : quadrature_t::quadrature) {
            xi_t v = qv.get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = qv.get_w() / 2.;
 
            // loop over u
            for (auto const &qu : quadrature_t::quadrature) {
                xi_t u = qu.get_xi();
                double wu = qu.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(xi).norm();
                    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 laplace_2d_DLPt_galerkin_face_general
 
 
template <class TestField, class TrialField, size_t order>
class laplace_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<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;
 
    typedef laplace_2d_HSP_kernel kernel_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:
    static result_t eval(elem_t const &elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_HSP_galerkin_face_general" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
 
        kernel_t kernel;
 
        result_t I = result_t::Zero();
 
        // loop over v
        for (auto const &qv : quadrature_t::quadrature) {
            xi_t v = qv.get_xi();
            v(0) = (v(0) + 1.) / 2.;
            double wv = qv.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);
                    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 laplace_2d_HSP_galerkin_face_general
 
 
class laplace_2d_HSP_galerkin_face_constant_line
{
    typedef line_1_elem elem_t;
public:
    static double eval(elem_t const &elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_HSP_galerkin_face_constant_line" << std::endl;
            printed = true;
        }
#endif
 
        using namespace boost::math::double_constants;
        double J = elem.get_normal().norm();
        return -(1.0 + std::log(2.0 * J)) / pi;
    } // end of function eval
}; // class laplace_2d_HSP_galerkin_face_constant_line
 
 
template <class TestField, class TrialField, size_t order>
class laplace_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<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;
 
    typedef laplace_2d_HSP_kernel kernel_t;
 
public:
    static result_t eval(
        test_elem_t const &test_elem,
        trial_elem_t const &trial_elem)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_2d_HSP_galerkin_edge_general" << std::endl;
            printed = true;
        }
#endif
        using namespace boost::math::double_constants;
 
        kernel_t kernel;
 
        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 const &qu : quadrature_t::quadrature) {
                xi_t u = qu.get_xi();
                double wu = qu.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 const &qv : quadrature_t::quadrature) {
                    xi_t v = qv.get_xi();
                    // transform v to (0,1) from (-1,1)
                    v(0) = (v(0) + 1.) / 2.;
                    double wv = qv.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);
                    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 laplace_2d_HSP_galerkin_edge_general
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_SLP_kernel, 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<TestField>::value &&
    is_constant_line<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_2d_SLP_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = laplace_2d_SLP_galerkin_face_constant_line::eval(trial_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_SLP_kernel, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    is_linear_line<TrialField>::value &&
    is_linear_line<TestField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_2d_SLP_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        laplace_2d_SLP_galerkin_face_linear_line::eval(trial_field.get_elem(), result(0, 0), result(0, 1));
        result(1, 0) = result(0, 1);
        result(1, 1) = result(0, 0);
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_SLP_kernel, TestField, TrialField, match::match_1d_type,
    typename std::enable_if<
    std::is_same<typename get_formalism<TestField, TrialField>::type, formalism::general>::value &&
    !(is_linear_line<TrialField>::value &&is_linear_line<TestField>::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<laplace_2d_SLP_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = laplace_2d_SLP_galerkin_face_general<TestField, TrialField, 11>::eval(trial_field.get_elem());
        return result;
    }
};
 
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_SLP_kernel, TestField, TrialField, match::match_0d_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<laplace_2d_SLP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = laplace_2d_SLP_galerkin_edge_constant_line::eval(
            test_field.get_elem(), trial_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_SLP_kernel, TestField, TrialField, match::match_0d_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<laplace_2d_SLP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = laplace_2d_SLP_galerkin_edge_general<
            TestField, TrialField, 8
        >::eval(test_field.get_elem(), trial_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_DLP_kernel, 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<laplace_2d_DLP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = laplace_2d_DLP_galerkin_face_general<TestField, TrialField, 11>::eval(
            test_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_DLP_kernel, TestField, TrialField, match::match_0d_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<laplace_2d_DLP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = laplace_2d_DLP_galerkin_edge_constant_line::eval(
            test_field.get_elem(), trial_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_DLPt_kernel, 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<laplace_2d_DLPt_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &,
        element_match const &)
    {
        result = laplace_2d_DLPt_galerkin_face_general<TestField, TrialField, 11>::eval(
            test_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_HSP_kernel, 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<TestField>::value &&is_constant_line<TrialField>::value)
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_2d_HSP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &,
        element_match const &)
    {
        result = laplace_2d_HSP_galerkin_face_general<TestField, TrialField, 10>::eval(
            test_field.get_elem()
        );
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_HSP_kernel, 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<TestField>::value &&is_constant_line<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_2d_HSP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &,
        element_match const &)
    {
        result(0, 0) = laplace_2d_HSP_galerkin_face_constant_line::eval(test_field.get_elem());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_2d_HSP_kernel, 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<laplace_2d_HSP_kernel> const &,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result = laplace_2d_HSP_galerkin_edge_general<
            TestField, TrialField, 11>::eval(
                test_field.get_elem(), trial_field.get_elem());
        return result;
    }
};
 
} // namespace NiHu
 
#endif /* LAPLACE_2D_SINGULAR_DOUBLE_INTEGRALS_HPP_INCLUDED */