x2p_integral.hpp

Go to the documentation of this file.

 
#ifndef NIHU_X2P_INTEGRAL_HPP_INCLUDED
#define NIHU_X2P_INTEGRAL_HPP_INCLUDED
 
#include "core/gaussian_quadrature.hpp"
#include "core/field.hpp"
#include "fmm_operator.hpp"
#include "integral_operator_expression.hpp"
#include "kron_identity.hpp"
#include "core/jacobian_computer.hpp"
#include "util/matrix_traits.hpp"
#include "util/type2tag.hpp"
 
#include <type_traits> // std::enable_if
 
 
namespace NiHu
{
namespace fmm
{
 
template <class Operator, class TestField>
class x2p_integral;
 
template <class Operator, class TestField>
struct integral_operator_expression_traits<x2p_integral<Operator, TestField> >
{
    typedef typename std::decay<Operator>::type operator_t;
    typedef TestField test_input_t;
    typedef typename operator_t::trial_input_t trial_input_t;
    typedef Eigen::Matrix<
        typename scalar<typename operator_t::result_t>::type,
        num_rows<typename operator_t::result_t>::value *TestField::nset_t::num_nodes,
        num_cols<typename operator_t::result_t>::value
    > result_t;
};
 
 
template <class Operator, class TestField>
class x2p_integral
    : public integral_operator_expression<x2p_integral<Operator, TestField> >
    , public fmm_operator<typename std::decay<Operator>::type::fmm_tag>
{
public:
    typedef integral_operator_expression<x2p_integral<Operator, TestField> > base_t;
 
    typedef typename std::decay<Operator>::type operator_t;
    typedef TestField test_field_t;
 
    typedef typename base_t::trial_input_t trial_input_t;
    typedef typename base_t::test_input_t test_input_t;
    typedef typename base_t::result_t result_t;
 
    typedef typename test_field_t::elem_t elem_t;
    typedef typename test_field_t::nset_t nset_t;
 
    typedef typename elem_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef NiHu::gaussian_quadrature<domain_t> quadrature_t;
 
    static size_t const op_num_rows = num_rows<typename operator_t::result_t>::value;
 
    static size_t const result_rows =
        op_num_rows * nset_t::num_nodes;
 
    x2p_integral(Operator &&op, size_t order)
        : m_op(std::forward<Operator>(op))
        , m_quadrature(unsigned(order))
    {
    }
 
    size_t rows(test_input_t const &) const
    {
        return result_rows;
    }
 
    size_t cols(trial_input_t const &from) const
    {
        return m_op.cols(from);
    }
 
    result_t operator()(test_input_t const &to, trial_input_t const &from) const
    {
#if 0
        static bool printed = false;
#endif
 
        result_t mat = result_t::Zero(result_rows, cols(from));
        elem_t const &elem = to.get_elem();
        // traverse quadrature elements
        for (auto const &q : m_quadrature)
        {
            xi_t const &xi = q.get_xi();
            double w = q.get_w();
            typename operator_t::test_input_t tsi(elem, xi);
            double jac = jacobian_computer<elem_t>::eval(elem, xi);
            auto N = (nset_t::template eval_shape<0>(xi) * (w * jac)).eval();
 
            mat += create_kron_identity<op_num_rows>(N) * m_op(tsi, from);
#if 0
            if (!printed) {
                printed = true;
                mexPrintf("Size of kmat in m2p integral: %d x %d\n", mat.rows(), mat.cols());
            }
#endif
#if 0
            auto op = m_op(tsi, from);
            for (Eigen::Index i = 0; i < N.rows(); ++i)
                mat.block(i * op_num_rows, 0, op_num_rows, cols(from)) += N(i, 0) * op;
 
            if (!printed) {
                printed = true;
                mexPrintf("Size of op in m2p integral: %d x %d\n", op.rows(), op.cols());
            }
#endif
        }
        return mat;
    }
 
private:
    Operator m_op;
    quadrature_t m_quadrature;
};
 
 
 
template <class Operator, class TestField>
class x2p_integral<Operator, NiHu::dirac_field<TestField> >
    : public integral_operator_expression<x2p_integral<Operator, NiHu::dirac_field<TestField> > >
    , public fmm_operator<typename std::decay<Operator>::type::fmm_tag>
{
public:
    typedef typename std::decay<Operator>::type operator_t;
    typedef TestField test_field_t;
 
    typedef integral_operator_expression<x2p_integral<Operator, NiHu::dirac_field<TestField> > > base_t;
 
    typedef typename base_t::trial_input_t trial_input_t;
    typedef typename base_t::test_input_t test_input_t;
    typedef typename base_t::result_t result_t;
 
    typedef typename test_field_t::elem_t elem_t;
    typedef typename test_field_t::nset_t nset_t;
 
    typedef typename elem_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef NiHu::gaussian_quadrature<domain_t> quadrature_t;
 
    static size_t const result_rows =
        num_rows<typename operator_t::result_t>::value * nset_t::num_nodes;
 
    x2p_integral(Operator &&op, size_t order = 0)
        : m_op(std::forward<Operator>(op))
    {
    }
 
    size_t rows(test_input_t const &) const
    {
        return result_rows;
    }
 
    size_t cols(trial_input_t const &from) const
    {
        return m_op.cols(from);
    }
 
    result_t operator()(test_input_t const &to, trial_input_t const &from) const
    {
        result_t mat = result_t::Zero(result_rows, cols(from));
        auto const M = num_rows<typename operator_t::result_t>::value;
        for (size_t N = 0; N < nset_t::num_nodes; ++N)
        {
            // get the node in intrinsic coordinates
            xi_t const &xi = nset_t::corner_at(N);
            // form the test input
            typename operator_t::test_input_t tsi(to.get_elem(), xi);
            // compute the operator and place in result's block
            mat.block(N * M, 0, M, cols(from)) = m_op(tsi, from);
        }
 
        return mat;
    }
 
private:
    Operator m_op;
};
 
template <class Operator, class FieldTag>
x2p_integral<Operator, typename tag2type<FieldTag>::type>
create_x2p_integral(Operator &&op, size_t order, FieldTag)
{
    typedef typename tag2type<FieldTag>::type test_field_t;
    return x2p_integral<Operator, test_field_t>(std::forward<Operator>(op), order);
}
 
} // end of namespace fmm
} // namespace NiHu
 
#endif /* NIHU_X2P_INTEGRAL_HPP_INCLUDED */