nearly_singular_collocational_telles.hpp

#ifndef NIHU_NEARLY_SINGULAR_COLLOCATIONAL_TELLES_HPP_INCLUDED
#define NIHU_NEARLY_SINGULAR_COLLOCATIONAL_TELLES_HPP_INCLUDED
 
#include "line_quad_store.hpp"
#include "location_normal.hpp"
#include "weighted_input.hpp"
#include "telles_1987.hpp"
#include "../core/element.hpp"
#include "../core/field.hpp"
#include "../core/inverse_mapping.hpp"
#include "../core/kernel.hpp"
#include "../core/nearly_singular_planar_constant_collocation_shortcut.hpp"
#include "../core/shapeset.hpp"
#include "../util/block_product.hpp"
 
namespace NiHu
{
 
namespace internal
{
 
template <class Domain>
struct telles_transformer
{
    typedef Domain domain_t;
    typedef gaussian_quadrature<domain_t> quad_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename domain_t::scalar_t scalar_t;
 
    telles_transformer(unsigned Order)
        : m_quad_orig(Order)
    {
 
    }
 
    quad_t transform(xi_t const& xi0, scalar_t const& zeta0)
    {
        return telles_transform(m_quad_orig, xi0, zeta0);
    }
 
private:
    quad_t const m_quad_orig;
};
 
template <>
struct telles_transformer<tria_domain>
{
    typedef tria_domain domain_t;
    typedef gaussian_quadrature<domain_t> tria_quad_t;
    typedef gaussian_quadrature<line_domain> line_quad_t;
 
    typedef typename domain_t::xi_t xi_t;
    typedef typename domain_t::scalar_t scalar_t;
 
    telles_transformer(unsigned Order)
        : m_quad_line(Order)
    {
    }
 
    tria_quad_t transform(xi_t const& xi0, scalar_t const& zeta0)
    {
        tria_quad_t quad_tria;
        size_t const n = m_quad_line.size();
        quad_tria.resize(n*n);
 
        typename line_domain::xi_t xi1_bar(2.0*xi0(0) - 1.0);
        // The jacobian of the re-transformation in xi_1 
        scalar_t jac1 = 2.0;
        
        // Transform the original line quadrature
        line_quad_t xi1_quad = telles_transform(m_quad_line, xi1_bar, zeta0 * jac1);
 
        for (size_t i = 0; i < xi1_quad.size(); ++i) {
            // Actual coordinate in xi1
            scalar_t xi1 = (1.0 + xi1_quad[i].get_xi()(0)) / 2.0;
            // Re-transform the singular point in xi_2 from [0, 1-xi1] to [-1, 1]
            typename line_domain::xi_t xi2_bar(2.0*xi0(1) / (1.0 - xi1) - 1.0);
            // The jacobian of the re-transformation in xi_2
            scalar_t jac2 = 2.0 / (1.0 - xi1);
 
            line_quad_t xi2_quad = telles_transform(m_quad_line, xi2_bar, zeta0 * jac2);
 
            // Fill the quadrature points on the triangle
            for (size_t j = 0; j < xi2_quad.size(); ++j) {
                scalar_t eta = (1.0 - xi1)*(xi2_quad[j].get_xi()(0) + 1.0)/2.0;
                quad_tria[i*xi2_quad.size() + j].set_xi(xi_t(xi1, eta));
                quad_tria[i*xi2_quad.size() + j].set_w(xi1_quad[i].get_w() / jac1 * xi2_quad[j].get_w() / jac2);
            }
        }
 
        return quad_tria;
    }
 
private:
    line_quad_t const m_quad_line;
};
 
 
} // of namespace internal
 
template <class TrialField, class Kernel, unsigned Order, class Enable = void>
class nearly_singular_collocational_telles;
 
template <class TrialField, class Kernel, unsigned Order>
class nearly_singular_collocational_telles<TrialField, Kernel, Order,
    typename std::enable_if<
    element_traits::is_surface_element<typename TrialField::elem_t>::value
    >::type>
{
public:
    enum {
        order = Order,          
    };
 
    typedef TrialField trial_field_t;
    typedef typename trial_field_t::elem_t elem_t;
 
    typedef typename elem_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef typename elem_t::scalar_t scalar_t;
 
    typedef typename elem_t::x_t x_t;
 
    typedef typename trial_field_t::nset_t trial_nset_t;
    typedef typename trial_nset_t::shape_t trial_n_shape_t;
 
    typedef Kernel kernel_t;
 
    typedef typename kernel_traits<kernel_t>::test_input_t test_input_t;
    typedef typename kernel_traits<kernel_t>::trial_input_t trial_input_t;
 
    typedef typename weighted_input<trial_input_t, elem_t>::type w_trial_input_t;
 
    typedef typename semi_block_product_result_type<
        typename kernel_t::result_t, trial_n_shape_t
    >::type total_result_t;
 
    // class for computing the transformed quadrature
    typedef internal::telles_transformer<domain_t> transformer_t;
 
 
    typedef gaussian_quadrature<domain_t> quad_t;
 
    nearly_singular_collocational_telles(
        field_base<trial_field_t> const &trial_field,
        kernel_base<kernel_t> const &kernel)
        : m_elem(trial_field.get_elem())
        , m_kernel(kernel)
        , m_transformer(Order)
    {
    }
 
    template <class result_t>
    void integrate(result_t &&I, test_input_t const &tsi)
    {
        // the reference point
        x_t x0 = tsi.get_x();
 
        // perform inverse mapping to obtain image of reference point on element
        double tol = 1e-6;
        unsigned max_iter = 100;
        // inverse mapping constrained inside element
        inverse_mapping<elem_t> im(m_elem, true);
        if (!im.eval(x0, tol, max_iter))
            throw std::runtime_error("Could not perform inverse mapping");
        auto res = im.get_result();
        m_xi0 = res.topRows(xi_t::RowsAtCompileTime);
        m_zeta = res(xi_t::RowsAtCompileTime);
 
        // perform telles transform
        quad_t quad_telles = m_transformer.transform(m_xi0, m_zeta);
 
        for (auto q : quad_telles)
        {
            xi_t xi = q.get_xi();
            double w = q.get_w();
            w_trial_input_t wtri(m_elem, xi);
            auto N = trial_nset_t::template eval_shape<0>(xi);
            I += semi_block_product(m_kernel(tsi, wtri) * wtri.get_jacobian() * w, N);
        }
 
    } // end of function integrate
 
private:
    elem_t const &m_elem;       
    kernel_base<kernel_t> const &m_kernel;  
    xi_t m_xi0;
    double m_zeta;
    transformer_t m_transformer;
};
 
} // end of namespace NiHu
 
#endif /* NIHU_NEARLY_SINGULAR_COLLOCATIONAL_HPP_INCLUDED */