single_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_SINGLE_INTEGRAL_HPP_INCLUDED
#define NIHU_SINGLE_INTEGRAL_HPP_INCLUDED
 
#include "global_definitions.hpp"
#include "../util/store_pattern.hpp"
#include "../util/plain_type.hpp"
#include "../util/product_type.hpp"
#include "../util/block_product.hpp"
#include "../util/materialize_sequence.hpp"
#include "gaussian_quadrature.hpp"
#include "field_type_accelerator.hpp"
#include "formalism.hpp"
#include "jacobian_computer.hpp"
 
namespace NiHu
{
 
template <class TestField, class TrialField>
struct single_integral_traits
{
    typedef typename block_product_result_type<
        typename TestField::nset_t::shape_t,
        Eigen::Matrix<double, TestField::quantity_dimension, TrialField::quantity_dimension>,
        typename TrialField::nset_t::shape_t
    >::type result_t;
};
 
 
template <class TestField, class TrialField, class Formalism = typename get_formalism<TestField, TrialField>::type>
class single_integral_impl
{
public:
    typedef TestField test_field_t;
    typedef TrialField trial_field_t;
 
    typedef single_integral_traits<test_field_t, trial_field_t> traits_t;
    typedef typename traits_t::result_t result_t;
 
    typedef typename test_field_t::elem_t elem_t;
    typedef typename elem_t::lset_t lset_t;
 
    typedef gauss_family_tag quadrature_family_t;
 
    typedef typename test_field_t::nset_t test_nset_t;
    typedef typename trial_field_t::nset_t trial_nset_t;
 
    static unsigned const degree
        = test_nset_t::polynomial_order
        + trial_nset_t::polynomial_order
        + lset_t::jacobian_order;
 
public:
    static result_t eval(
        field_base<test_field_t> const &test_field,
        field_base<trial_field_t> const &trial_field)
    {
        result_t result;
        result.setZero();
 
        typedef store<field_type_accelerator_pool<
                test_field_t, quadrature_family_t, GLOBAL_ACCELERATION, GLOBAL_MAX_ORDER
        > > test_store_t;
 
        typedef store<field_type_accelerator_pool<
                trial_field_t, 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::diagonal());
            
        // identity matrix instance for block product integration
        enum {Rows = TestField::quantity_dimension, Cols = TrialField::quantity_dimension};
        auto mat = Eigen::Matrix<double, Rows, Cols>::Identity();
 
        static bool printed = false;
 
        if (!printed) {
            std::cout << "single_integral_impl::eval() called" << std::endl;
            printed = true;
        }
 
        for (auto it = acc.begin(); it != acc.end(); ++it)
        {
            auto jac = NiHu::jacobian_computer<elem_t>::eval(test_field.get_elem(), it.get_first()->get_xi());
            
            /*
            result += block_product(it.get_first()->get_N(),
                (mat * it.get_first()->get_w()*jac).eval(),
                it.get_second()->get_N());
            */
 
            result += block_product(
                test_field.eval(it.get_first()->get_xi()),
                (mat * it.get_first()->get_w()*jac).eval(),
                trial_field.eval(it.get_second()->get_xi())
            );
        }
 
        return result;
    }
};
 
 
template <class TestField, class TrialField>
class single_integral_impl<TestField, TrialField, formalism::collocational>
{
public:
    typedef typename TestField::nset_t test_nset_t;
    typedef typename TrialField::nset_t trial_nset_t;
 
    typedef typename TrialField::elem_t trial_elem_t;
    typedef typename element_traits::lset<trial_elem_t>::type trial_lset_t;
 
    typedef single_integral_traits<TestField, TrialField> traits_t;
    typedef typename traits_t::result_t result_t;
 
    static result_t eval(
        field_base<TestField> const &test_field,
        field_base<TrialField> const &trial_field)
    {
        result_t result;
        result.setZero();
        
        enum {Rows = TestField::quantity_dimension, Cols = TrialField::quantity_dimension};
        auto mat = Eigen::Matrix<double, Rows, Cols>::Identity();
        
        for (unsigned row = 0; row < test_nset_t::num_nodes; ++row)
        {
            auto N = trial_field.eval(test_nset_t::corner_at(row));
            double weight = test_field.get_corner_angle(row);
 
            for (unsigned col = 0; col < trial_nset_t::num_nodes; ++col)
                result.template block<Rows, Cols>(row*Rows, col*Cols) = weight * mat * N(col);
        }
        
        return result;
    }
};
 
 
template <class TestField, class TrialField, class = void>
class single_integral
{
public:
    typedef typename single_integral_traits<TestField, TrialField>::result_t result_t;
 
    static constexpr result_t eval(
        field_base<TestField> const &,
        field_base<TrialField> const &)
    {
        return result_t::Zero();
    }
};
 
 
template <class TestField, class TrialField>
class single_integral<TestField, TrialField,
    typename std::enable_if<
        std::is_same<
            typename TestField::elem_t,
            typename TrialField::elem_t
        >::value
    >::type
> : public single_integral_impl<TestField, TrialField> {};
 
}
 
#endif /* NIHU_SINGLE_INTEGRAL_HPP_INCLUDED */