laplace_3d_singular_integrals.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-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_3D_SINGULAR_INTEGRALS_HPP_INCLUDED
#define LAPLACE_3D_SINGULAR_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 "guiggiani_1992.hpp"
#include "lib_element.hpp"
#include "laplace_kernel.hpp"
#include "normal_derivative_singular_integrals.hpp"
#include "plane_element_helper.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
{
 
class laplace_3d_SLP_collocation_constant_plane
{
public:
    template <class elem_t>
    static double eval(elem_t const &elem, typename elem_t::x_t const &x0)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_3d_SLP_collocation_constant_plane" << std::endl;
            printed = true;
        }
#endif
 
        enum { N = elem_t::domain_t::num_corners };
        double r[N], theta[N], alpha[N], result = 0.;
        plane_element_helper(elem, x0, r, theta, alpha);
 
        for (unsigned i = 0; i < N; ++i)
            result += r[i] * std::sin(alpha[i]) *
            std::log(std::tan((alpha[i] + theta[i]) / 2.) / std::tan(alpha[i] / 2.));
 
        return result / (4. * pi);
    }
};
 
 
class laplace_3d_HSP_collocation_constant_plane
{
public:
    template <class elem_t>
    static double eval(elem_t const &elem, typename elem_t::x_t const &x0)
    {
        using namespace boost::math::double_constants;
 
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using laplace_3d_HSP_collocation_constant_plane" << std::endl;
            printed = true;
        }
#endif      
 
        enum { N = elem_t::domain_t::num_corners };
        double r[N], theta[N], alpha[N], result = 0.;
        plane_element_helper(elem, x0, r, theta, alpha);
 
        for (unsigned i = 0; i < N; ++i)
            result += (std::cos(alpha[i] + theta[i]) - std::cos(alpha[i])) / (r[i] * std::sin(alpha[i]));
 
        return result / (4. * pi);
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_3d_SLP_kernel, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::value &&
    is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_3d_SLP_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = laplace_3d_SLP_collocation_constant_plane::eval(
            trial_field.get_elem(),
            trial_field.get_elem().get_center());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_3d_HSP_kernel, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::value &&
    is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_3d_HSP_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
        result(0, 0) = laplace_3d_HSP_collocation_constant_plane::eval(
            trial_field.get_elem(),
            trial_field.get_elem().get_center());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_3d_Gxx_kernel, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::value &&
    is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_3d_Gxx_kernel> const &,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
 
        result(0, 0) = -1. * laplace_3d_HSP_collocation_constant_plane::eval(
            trial_field.get_elem(),
            trial_field.get_elem().get_center());
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class singular_integral_shortcut<
    laplace_3d_HSP_kernel, TestField, TrialField, match::match_2d_type,
    typename std::enable_if<
    is_collocational<TestField, TrialField>::value &&
    !is_constant_tria<TrialField>::value
    >::type
>
{
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<laplace_3d_HSP_kernel> const &kernel,
        field_base<TestField> const &,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
#ifdef NIHU_PRINT_DEBUG
        static bool printed = false;
        if (!printed) {
            std::cout << "Using guiggiani for laplace_3d_HSP_kernel" << std::endl;
            printed = true;
        }
#endif
        typedef guiggiani<TrialField, laplace_3d_HSP_kernel, 5, 9> guiggiani_t;
        auto const &elem = trial_field.get_elem();
        guiggiani_t gui(elem, kernel.derived());
        for (unsigned r = 0; r < TestField::num_dofs; ++r)
        {
            auto const &xi0 = TestField::nset_t::corner_at(r);
            gui.integrate(result.row(r), xi0, elem.get_normal(xi0));
        }
        return result;
    }
};
 
} // namespace NiHu
 
#endif /* LAPLACE_3D_SINGULAR_INTEGRALS_HPP_INCLUDED */