integral_operator.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 INTEGRAL_OPERATOR_HPP_INCLUDED
#define INTEGRAL_OPERATOR_HPP_INCLUDED
 
#include "../util/crtp_base.hpp"
#include "function_space.hpp"
#include "single_integral.hpp"
#include "double_integral.hpp"
#include "integral_transform_forward.hpp"
 
namespace NiHu
{
 
 
template <class Derived>
struct integral_operator_traits;
 
template <class Derived>
class integral_operator_base
{
public:
    NIHU_CRTP_HELPERS
 
    typedef integral_operator_traits<Derived> traits_t;
 
    template <class TestField, class TrialField>
    struct wr_result_type : traits_t::template wr_result_type<TestField, TrialField> {};
 
    template <class TestField, class TrialField, class OnSameMesh>
    typename wr_result_type<TestField, TrialField>::type
        eval_on_fields(
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            OnSameMesh) const
    {
        return derived().derived_eval_on_fields(test_field, trial_field, OnSameMesh());
    }
 
    template <class FuncSpace>
    integral_transform<
        Derived,
        typename std::enable_if<is_function_space<FuncSpace>::value, FuncSpace>::type
    >
        operator[](FuncSpace &&funcspace)
    {
        return integral_transform<Derived, FuncSpace>(
            derived(),
            std::forward<FuncSpace>(funcspace));
    }
};
 
template <class IntOp>
struct is_integral_operator : std::is_base_of<
    integral_operator_base<typename std::decay<IntOp>::type>,
    typename std::decay<IntOp>::type
>{};
 
 
// forward declaration
template <class Scalar, class IntOp>
class scaled_integral_operator;
 
 
template <class Scalar, class IntOp>
struct integral_operator_traits<scaled_integral_operator<Scalar, IntOp> >
{
    template <class TestField, class TrialField>
    struct wr_result_type : plain_type<
        typename product_type<
        Scalar,
        typename integral_operator_traits<
            typename std::decay<IntOp>::type
        >::template wr_result_type<TestField, TrialField>::type
        >::type
    > {};
 
    static bool const is_local = integral_operator_traits<
        typename std::decay<IntOp>::type
    >::is_local;
};
 
 
template <class Scalar, class IntOp>
class scaled_integral_operator :
    public integral_operator_base<scaled_integral_operator<Scalar, IntOp> >
{
public:
    typedef integral_operator_base<scaled_integral_operator<Scalar, IntOp> > base_t;
 
 
    scaled_integral_operator(
        Scalar &&scalar,
        IntOp &&parent) :
        m_scalar(std::forward<Scalar>(scalar)),
        m_parent(std::forward<IntOp>(parent))
    {
    }
 
    template <class TestField, class TrialField, class OnSameMesh>
    typename base_t::template wr_result_type<TestField, TrialField>::type
        derived_eval_on_fields(
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            OnSameMesh) const
    {
        return m_scalar * m_parent.eval_on_fields(test_field, trial_field, OnSameMesh());
    }
 
private:
    Scalar m_scalar;
    IntOp m_parent;
};
 
 
template <class Scalar, class IntOp>
scaled_integral_operator<
    Scalar,
    typename std::enable_if<is_integral_operator<IntOp>::value, IntOp>::type
>
    operator*(Scalar &&scalar, IntOp &&intop)
{
    return scaled_integral_operator<Scalar, IntOp>(
        std::forward<Scalar>(scalar),
        std::forward<IntOp>(intop));
}
 
 
 
template <class LhsOp, class RhsOp>
class sum_integral_operator;
 
template <class LhsOp, class RhsOp>
struct integral_operator_traits<sum_integral_operator<LhsOp, RhsOp>>
{
    template <class TestField, class TrialField>
    struct wr_result_type : plain_type<
        typename sum_type<
            typename integral_operator_traits<
                typename std::decay<LhsOp>::type
            >::template wr_result_type<TestField, TrialField>::type,
            typename integral_operator_traits<
                typename std::decay<RhsOp>::type
            >::template wr_result_type<TestField, TrialField>::type
        >::type
    > {};
 
    static bool const is_local = integral_operator_traits<
        typename std::decay<LhsOp>::type
    >::is_local && integral_operator_traits<
        typename std::decay<RhsOp>::type
    >::is_local;
 
    static_assert(
        integral_operator_traits<typename std::decay<LhsOp>::type>::is_local ==
        integral_operator_traits<typename std::decay<RhsOp>::type>::is_local
    , "You cannot form the sum of operators with different locality");
};
 
 
template <class LhsOp, class RhsOp>
class sum_integral_operator :
    public integral_operator_base<sum_integral_operator<LhsOp, RhsOp> >
{
public:
    static_assert(
        integral_operator_traits<typename std::decay<LhsOp>::type>::is_local ==
        integral_operator_traits<typename std::decay<RhsOp>::type>::is_local
    , "You cannot form the sum of operators with different locality");
 
    typedef integral_operator_base<sum_integral_operator<LhsOp, RhsOp> > base_t;
 
 
    sum_integral_operator(
        LhsOp &&lhs,
        RhsOp &&rhs) :
        m_lhs(std::forward<LhsOp>(lhs)),
        m_rhs(std::forward<RhsOp>(rhs))
    {
    }
 
    template <class TestField, class TrialField, class OnSameMesh>
    typename base_t::template wr_result_type<TestField, TrialField>::type
        derived_eval_on_fields(
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            OnSameMesh) const
    {
        return m_lhs.eval_on_fields(test_field, trial_field, OnSameMesh()) + m_rhs.eval_on_fields(test_field, trial_field, OnSameMesh());
    }
 
private:
    LhsOp m_lhs;
    RhsOp m_rhs;
};
 
 
template <class LhsOp, class RhsOp>
sum_integral_operator<
    typename std::enable_if<is_integral_operator<LhsOp>::value, LhsOp>::type,
    typename std::enable_if<is_integral_operator<RhsOp>::value, RhsOp>::type
>
    operator+(LhsOp &&lhs, RhsOp &&rhs)
{
    return sum_integral_operator<LhsOp, RhsOp>(
        std::forward<LhsOp>(lhs),
        std::forward<RhsOp>(rhs));
}
 
 
 
// forward declaration
class identity_integral_operator;
 
template <>
struct integral_operator_traits<identity_integral_operator>
{
    template <class TestField, class TrialField>
    struct wr_result_type
    {
        typedef typename single_integral<TestField, TrialField>::result_t type;
    };
 
    static bool const is_local = true;
};
 
class identity_integral_operator :
    public integral_operator_base<identity_integral_operator>
{
public:
    typedef integral_operator_base<identity_integral_operator> base_t;
 
    template <class TestField, class TrialField, class OnSameMesh = std::false_type>
    typename base_t::template wr_result_type<TestField, TrialField>::type
        derived_eval_on_fields(
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            OnSameMesh) const
    {
        return single_integral<TestField, TrialField>::eval(test_field, trial_field);
    }
};
 
 
 
// forward declaration
template <class Kernel>
class integral_operator;
 
template <class Kernel>
struct integral_operator_traits<integral_operator<Kernel> >
{
    template <class TestField, class TrialField>
    struct wr_result_type
    {
        typedef typename double_integral<typename std::decay<Kernel>::type, TestField, TrialField>::result_t type;
    };
 
    static bool const is_local = false;
};
 
 
 
template <class Kernel>
class integral_operator :
    public integral_operator_base<integral_operator<Kernel> >
{
public:
    typedef integral_operator_base<integral_operator<Kernel> > base_t;
 
    typedef typename std::decay<Kernel>::type kernel_t;
 
    integral_operator(Kernel &&kernel) :
        m_kernel(std::forward<Kernel>(kernel))
    {
    }
 
    Kernel get_kernel(void) const
    {
        return m_kernel;
    }
 
    template <class TestField, class TrialField, class OnSameMesh>
    typename base_t::template wr_result_type<TestField, TrialField>::type
        derived_eval_on_fields(
            field_base<TestField> const &test_field,
            field_base<TrialField> const &trial_field,
            OnSameMesh) const
    {
        return double_integral<kernel_t, TestField, TrialField>::eval(
            m_kernel, test_field, trial_field, OnSameMesh());
    }
 
private:
    Kernel m_kernel;
};
 
 
template <class Kernel>
integral_operator<Kernel>
    create_integral_operator(Kernel &&kernel)
{
    return integral_operator<Kernel>(std::forward<Kernel>(kernel));
}
 
} // end of namespace NiHu
 
#endif