double_integral.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 NIHU_DOUBLE_INTEGRAL_HPP_INCLUDED
#define NIHU_DOUBLE_INTEGRAL_HPP_INCLUDED
 
#include "../util/matrix_traits.hpp"
#include "../util/product_type.hpp"
#include "../util/plain_type.hpp"
#include "../util/block_product.hpp"
#include "../library/location_normal.hpp"
#include "../library/weighted_input.hpp"
#include "../tmp/control.hpp"
#include "../util/store_pattern.hpp"
#include "complexity_estimator.hpp"
#include "element_match.hpp"
#include "field_type_accelerator.hpp"
#include "formalism.hpp"
#include "global_definitions.hpp"
#include "kernel.hpp"
#include "match_types.hpp"
#include "nearly_singular_integral.hpp"
#include "singular_accelerator.hpp"
#include "singular_integral_shortcut.hpp"
 
 
namespace NiHu
{
 
 
/*
template <class Elem>
struct weighted_brick;
 
template <class LSet, class Scalar>
struct weighted_brick<volume_element<LSet, Scalar> >
{
    typedef volume_jacobian<typename volume_element<LSet, Scalar>::space_t> type;
};
 
template <class LSet, class Scalar>
struct weighted_brick<surface_element<LSet, Scalar> >
{
    typedef normal_jacobian<typename surface_element<LSet, Scalar>::space_t> type;
};
*/
 
 
template <class SingularityDimension>
struct singular_shortcut_switch
{
    struct type
    {
        template <class result_t, class Kernel, class TestField, class TrialField>
        bool operator()(
            result_t &result,
            kernel_base<Kernel> const &kernel,
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            element_match const &mtch)
        {
            // if the parameter singularity is valid, evaluate shortcut
            if (mtch.get_match_dimension() == SingularityDimension::value)
            {
#if NIHU_DEBUGGING
                static bool printed = false;
                if (!printed)
                {
                    std::cout <<
                        "Singular shortcut switch called for mtch dim: " <<
                        mtch.get_match_dimension() <<
                        ", sing val: " << SingularityDimension::value << std::endl;
                    printed = true;
                }
#endif
                singular_integral_shortcut<Kernel, TestField, TrialField, SingularityDimension>::eval(
                    result, kernel, test_field, trial_field, mtch);
                return true;
            }
 
            return false;
        }
    };
};
 
 
 
template <class Kernel, class TestField, class TrialField>
struct double_integral_traits
{
    typedef typename block_product_result_type<
        typename TestField::nset_t::shape_t,
        typename Kernel::result_t,
        typename TrialField::nset_t::shape_t
    >::type result_t;
};
 
 
 
template <
    class Kernel, class TestField, class TrialField,
    class Formalism = typename get_formalism<TestField, TrialField>::type
>
class double_integral;
 
 
template <class Kernel, class TestField, class TrialField>
class double_integral<Kernel, TestField, TrialField, formalism::general>
{
    typedef double_integral_traits<Kernel, TestField, TrialField> traits_t;
 
public:
    typedef std::true_type WITH_SINGULARITY_CHECK;
    typedef std::false_type WITHOUT_SINGULARITY_CHECK;
 
    typedef typename TestField::elem_t test_elem_t;
    typedef typename TrialField::elem_t trial_elem_t;
    typedef typename kernel_traits<Kernel>::test_input_t test_input_t;
    typedef typename kernel_traits<Kernel>::trial_input_t trial_input_t;
    typedef typename weighted_input<test_input_t, test_elem_t>::type w_test_input_t;
    typedef typename weighted_input<trial_input_t, trial_elem_t>::type w_trial_input_t;
 
    typedef typename kernel_traits<Kernel>::quadrature_family_t quadrature_family_t;
 
    typedef typename traits_t::result_t result_t;
 
protected:
    template <class dual_iterator_t>
    static result_t &eval_on_accelerator(
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        dual_iterator_t it,
        dual_iterator_t end
    )
    {
        for (; it != end; ++it)
        {
            w_test_input_t test_input(test_field.get_elem(), it.get_first()->get_xi());
            w_trial_input_t trial_input(trial_field.get_elem(), it.get_second()->get_xi());
 
            result += block_product(
                test_field.eval(it.get_first()->get_xi()) *
                (test_input.get_jacobian() * it.get_first()->get_w()),
                kernel(test_input, trial_input),
                trial_field.eval(it.get_second()->get_xi()) *
                (trial_input.get_jacobian() * it.get_second()->get_w())
            );
        }
 
        return result;
    }
 
public:
    template <class singular_accelerator_t, class dummy>
    struct eval_singular_on_accelerator
    {
        template <class singular_iterator_t>
        static result_t &eval(
            result_t &result,
            kernel_base<Kernel> const &kernel,
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            singular_iterator_t begin,
            singular_iterator_t end)
        {
            for (; begin != end; ++begin)
            {
                w_test_input_t test_input(test_field.get_elem(), begin.get_test_quadrature_elem().get_xi());
                w_trial_input_t trial_input(trial_field.get_elem(), begin.get_trial_quadrature_elem().get_xi());
 
                result += block_product(
                    test_field.eval(begin.get_test_quadrature_elem().get_xi()) *
                    (test_input.get_jacobian() * begin.get_test_quadrature_elem().get_w()),
                    kernel(test_input, trial_input),
                    trial_field.eval(begin.get_trial_quadrature_elem().get_xi()) *
                    (trial_input.get_jacobian() * begin.get_trial_quadrature_elem().get_w())
                );
            }
 
            return result;
        }
    };
 
    template <class dummy>
    struct eval_singular_on_accelerator<invalid_singular_accelerator, dummy>
    {
        template <class singular_iterator_t>
        static result_t &eval(
            result_t &result,
            kernel_base<Kernel> const &,
            field_base<TestField> const &,
            field_base<TrialField> const &,
            singular_iterator_t,
            singular_iterator_t)
        {
            throw std::runtime_error("Invalid general singular quadrature");
            return result;
        }
    };
 
protected:
    static result_t &eval(
        WITHOUT_SINGULARITY_CHECK,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field)
    {
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("double_integral::eval without singularity check called on elements %d <- %d\n",
                test_field.get_dofs()(0), trial_field.get_dofs()(0));
            printed = true;
        }
#endif
        // store type of the regular test quadratures
        typedef store<field_type_accelerator_pool<
            TestField, quadrature_family_t, GLOBAL_ACCELERATION, GLOBAL_MAX_ORDER
        > > test_store_t;
 
        // store type of the regular trial quadratures
        typedef store<field_type_accelerator_pool<
            TrialField, quadrature_family_t, GLOBAL_ACCELERATION, GLOBAL_MAX_ORDER
        > > trial_store_t;
 
        // determine integration degree based on the complexity estimator
        unsigned degree = complexity_estimator<
            TestField, TrialField,
            typename Kernel::estimator_t
        >::eval(test_field, trial_field);
 
        auto acc(create_dual_field_type_accelerator(
            test_store_t::get_data()[degree], trial_store_t::get_data()[degree], iteration::diadic()));
 
        return eval_on_accelerator(
            result, kernel, test_field, trial_field, acc.begin(), acc.end());
    }
 
 
    static result_t &eval(
        WITH_SINGULARITY_CHECK,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field)
    {
        // evaluate element match
        auto mtch(element_match_eval(test_field, trial_field));
 
        // if no match simple regular integral is computed
        if (mtch.get_match_dimension() == -1)
            return eval(WITHOUT_SINGULARITY_CHECK(), result, kernel, test_field, trial_field);
 
        typedef typename match_type_vector<TestField, TrialField>::type possible_match_types;
 
        // traverse possible singular integral shortcuts with tmp::call_until
        if (!tmp::call_until<
            possible_match_types,
            singular_shortcut_switch<tmp::_1>,
            result_t &,
            kernel_base<Kernel> const &,
            field_base<TestField> const &,
            field_base<TrialField> const &,
            element_match const &
        >(result, kernel, test_field, trial_field, mtch))
            std::cerr << "UNHANDLED GALERKIN SINGULARITY TYPE: " << mtch.get_match_dimension() << std::endl;
 
        return result;
    }
 
public:
    template <class OnSameMesh>
    static result_t eval(
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        OnSameMesh)
    {
        static bool const sing_check_needed =
            kernel_traits<Kernel>::is_singular && std::is_same<OnSameMesh, std::true_type>::value;
 
        result_t result;
        result.setZero();   // clear result
 
        return eval(std::integral_constant<bool, sing_check_needed>(),
            result, kernel, test_field, trial_field);
    }
};
 
 
 
 
template <class Kernel, class TestField, class TrialField>
class double_integral<Kernel, TestField, TrialField, formalism::collocational>
{
 
    typedef double_integral_traits<Kernel, TestField, TrialField> traits_t;
 
public:
    typedef std::true_type WITH_SINGULARITY_CHECK;
    typedef std::false_type WITHOUT_SINGULARITY_CHECK;
 
    typedef typename TrialField::elem_t trial_elem_t;
    typedef typename kernel_traits<Kernel>::test_input_t test_input_t;
    typedef typename kernel_traits<Kernel>::trial_input_t trial_input_t;
    typedef typename kernel_traits<Kernel>::result_t kernel_result_t;
 
    typedef typename weighted_input<trial_input_t, trial_elem_t>::type w_trial_input_t;
 
    typedef typename kernel_traits<Kernel>::quadrature_family_t quadrature_family_t;
 
    typedef typename TestField::nset_t test_nset_t;
    typedef typename TrialField::nset_t trial_nset_t;
 
    typedef typename traits_t::result_t result_t;
 
    enum {
        kernel_rows = num_rows<kernel_result_t>::value,
        kernel_cols = num_cols<kernel_result_t>::value
    };
 
protected:
    template <class dual_iterator_t>
    static result_t &eval_on_accelerator(
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        dual_iterator_t first,
        dual_iterator_t last)
    {
        for (; first != last; ++first)
        {
            test_input_t test_input(test_field.get_elem(), first.get_first()->get_xi());
            w_trial_input_t trial_input(trial_field.get_elem(), first.get_second()->get_xi());
 
            /*
            result += block_product(
                first.get_first()->get_N() * (first.get_first()->get_w()),
                kernel(test_input, trial_input),
                first.get_second()->get_N() * (trial_input.get_jacobian() * first.get_second()->get_w())
            );
            */
            result += block_product(
                //test_field.eval(first.get_first()->get_xi()) * (first.get_first()->get_w()),
                first.get_first()->get_N() * (first.get_first()->get_w()),
                kernel(test_input, trial_input),
                trial_field.eval(first.get_second()->get_xi()) * (trial_input.get_jacobian() * first.get_second()->get_w())
            );
 
        }
 
        return result;
    }
 
public:
    template <class singular_accelerator_t, class dummy>
    struct eval_singular_on_accelerator
    {
        static result_t & eval(
            result_t &result,
            kernel_base<Kernel> const &kernel,
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            singular_accelerator_t const &sa)
        {
            for (unsigned idx = 0; idx < test_nset_t::num_nodes; ++idx)
            {
                test_input_t collocation_point(test_field.get_elem(), test_nset_t::corner_at(idx));
 
                for (auto const &q : sa.get_trial_quadrature(idx))
                {
                    w_trial_input_t trial_input(trial_field.get_elem(), q.get_xi());
 
                    result.template block<kernel_rows, kernel_cols*trial_nset_t::num_nodes>(idx*kernel_rows, 0)
                        += semi_block_product(
                            kernel(collocation_point, trial_input),
                            trial_field.eval(q.get_xi()) * trial_input.get_jacobian() * q.get_w()
                        );
                }
            }
 
            return result;
        }
    };
 
    template <class dummy>
    struct eval_singular_on_accelerator<invalid_singular_accelerator, dummy>
    {
        static result_t & eval(
            result_t &result,
            kernel_base<Kernel> const &,
            field_base<TestField> const &,
            field_base<TrialField> const &,
            invalid_singular_accelerator const &)
        {
            throw std::runtime_error("Invalid quadrature returned by eval_singular_on_accelerator");
            return result;
        }
    };
 
protected:
    static result_t &eval(
        WITHOUT_SINGULARITY_CHECK,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field)
    {
        // evaluate nearly with nearly singular shortcut if needed
        if (nearly_singular_integral<Kernel, TestField, TrialField>::needed(
            kernel, test_field, trial_field))
            return nearly_singular_integral<Kernel, TestField, TrialField>::eval(
                result, kernel, test_field, trial_field);
 
        // evaluate regular integral
        unsigned degree = complexity_estimator<
            TestField, TrialField,
            typename Kernel::estimator_t
        >::eval(test_field, trial_field);
 
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("double_integral::eval without singularity check called on elements %d <- %d\n",
                test_field.get_dofs()(0), trial_field.get_dofs()(0));
            mexPrintf("Regular integral with degree %d\n", degree);
            printed = true;
        }
#endif
 
        if (degree > GLOBAL_MAX_ORDER)
            throw std::out_of_range("Too high quadrature degree selected for collocational integration");
 
        typedef store<field_type_accelerator_pool<
            TestField, quadrature_family_t, GLOBAL_ACCELERATION, GLOBAL_MAX_ORDER
        > > test_store_t;
 
        typedef store<field_type_accelerator_pool<
            TrialField, quadrature_family_t, GLOBAL_ACCELERATION, GLOBAL_MAX_ORDER
        > > trial_store_t;
 
        auto acc(create_dual_field_type_accelerator(
            test_store_t::get_data()[degree], trial_store_t::get_data()[degree], iteration::diadic()));
 
        return eval_on_accelerator(
            result, kernel, test_field, trial_field, acc.begin(), acc.end());
    }
 
 
    static result_t &eval(
        WITH_SINGULARITY_CHECK,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field)
    {
 
        // perform singularity check and call regular integral if not singular
        auto mtch(element_match_eval(test_field, trial_field));
        if (mtch.get_match_dimension() == -1)
            return eval(WITHOUT_SINGULARITY_CHECK(), result, kernel, test_field, trial_field);
 
 
#if NIHU_DEBUGGING
        static bool printed = false;
        if (!printed)
            std::cout << "Now checking singular shortcuts" << std::endl;
#endif
 
        // call singular_shortcut_switch for each possible match type
        typedef typename match_type_vector<TestField, TrialField>::type possible_match_types;
        if (!tmp::call_until<
            possible_match_types,
            singular_shortcut_switch<tmp::_1>,
            result_t &,
            kernel_base<Kernel> const &,
            field_base<TestField> const &,
            field_base<TrialField> const &,
            element_match const &
        >(result, kernel, test_field, trial_field, mtch))
        {
            std::cerr << "UNHANDLED COLLOCATIONAL SINGULARITY TYPE: " << mtch.get_match_dimension() << std::endl;
        }
 
#if NIHU_DEBUGGING
        if (!printed)
            printed = true;
#endif
 
        return result;
    }
 
public:
    template <class OnSameMesh>
    static result_t eval(
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        OnSameMesh)
    {
        typedef typename TrialField::nset_t trial_nset_t;
 
        static unsigned const num_test_corners = shape_set_traits::num_nodes<typename TestField::nset_t>::value;
 
        static bool const sing_check_needed =
            kernel_traits<Kernel>::is_singular && std::is_same<OnSameMesh, std::true_type>::value;
 
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("double_integral<collocation>::eval called\nSingular kernel: %d\tSing needed: %d\n", kernel_traits<Kernel>::is_singular, sing_check_needed);
            printed = true;
        }
#endif
 
        result_t result = result_t::Zero();
 
        eval(
            std::integral_constant<bool, sing_check_needed>(),
            result, kernel, test_field, trial_field);
 
        // Weighting by corner angles
        for (unsigned corner_idx = 0; corner_idx < num_test_corners; ++corner_idx) 
        {
            double weight = test_field.get_corner_angle(corner_idx);
            // Block for each corner in result
            result.template block<kernel_rows, kernel_cols*trial_nset_t::num_nodes>(corner_idx*kernel_rows, 0) *= weight;
        }
 
        return result;
    }
};
 
 
template <class Kernel, class TestField, class TrialField, class SingularityDimension, class Enable>
class singular_integral_shortcut
{
public:
    typedef void unspecialised;
 
private:
 
    typedef typename get_formalism<TestField, TrialField>::type formalism_t;
    typedef typename select_singular_accelerator<
        Kernel, TestField, TrialField, SingularityDimension
    >::type singular_accelerator_t;
    typedef store<singular_accelerator_t> store_t;
    typedef double_integral<Kernel, TestField, TrialField> double_integral_t;
 
    template <class result_t>
    static result_t &eval_impl(
        formalism::general,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &match)
    {
        return double_integral_t::template eval_singular_on_accelerator<singular_accelerator_t, void>::eval(
            result, kernel, test_field, trial_field,
            store_t::get_data().begin(match),
            store_t::get_data().end(match));
    }
 
    template <class result_t>
    static result_t &eval_impl(
        formalism::collocational,
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &)
    {
#if NIHU_DEBUGGING
        static bool printed = false;
        if (!printed) {
            std::cout << "General version of singular_integral_shortcut called" << std::endl;
            printed = true;
        }
#endif
#if NIHU_MEX_DEBUGGING
    {
 
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("General singular integral shortcut for collocation called\n");
            printed = true;
        }
    }
#endif
        return double_integral_t::template eval_singular_on_accelerator<singular_accelerator_t, void>::eval(
            result, kernel, test_field, trial_field, store_t::get_data());
    }
 
public:
    template <class result_t>
    static result_t &eval(
        result_t &result,
        kernel_base<Kernel> const &kernel,
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field,
        element_match const &match)
    {
#if NIHU_MEX_DEBUGGING
        static bool printed = false;
        if (!printed)
        {
            mexPrintf("General singular integral shortcut called\n");
            printed = true;
        }
#endif
        return eval_impl(formalism_t(), result, kernel, test_field, trial_field, match);
    }
};
 
}
 
#endif /* NIHU_DOUBLE_INTEGRAL_HPP_INCLUDED */