duffy_quadrature.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 DUFFY_QUADRATURE_HPP_INCLUDED
#define DUFFY_QUADRATURE_HPP_INCLUDED
 
#include "quadrature.hpp"
#include "../library/lib_shape.hpp"
 
namespace NiHu
{
 
template <class LSet>
struct duffy_traits;
 
template <class QuadFamily, class LSet>
class duffy_quadrature
{
    // CRTP check
    static_assert(std::is_base_of<shape_set_base<LSet>, LSet>::value,
        "The LSet parameter must be derived from shape_set_base<LSet>");
public:
    typedef QuadFamily quadrature_family_tag;
    typedef LSet lset_t;
 
    typedef typename lset_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef typename domain_t::scalar_t scalar_t;
    typedef typename quadrature_type<quadrature_family_tag, domain_t>::type ret_quad_type;
 
private:
    typedef typename quadrature_type<quadrature_family_tag, quad_domain>::type source_quad_type;
    typedef Eigen::Matrix<scalar_t, quad_domain::num_corners, domain_t::dimension> coords_t;
 
    static ret_quad_type &duffy_impl(
        unsigned degree,
        unsigned const *array,
        xi_t const &sing_coord,
        ret_quad_type &result)
    {
        // the first array element is the number of triangles
        unsigned num_duffies = *array;
        array++;
 
        // create a regular quadrature
        source_quad_type source(degree);
 
        coords_t coords;
        coords.row(0) = coords.row(1) = sing_coord;
 
        for (size_t d = 0; d < num_duffies; ++d)
        {
            coords.row(2) = domain_t::get_corner(array[d]);
            coords.row(3) = domain_t::get_corner(array[d+1]);
 
            result += source.template transform<quad_1_shape_set>(coords);
        }
 
        return result;
    }
 
public:
    static ret_quad_type on_corner(unsigned degree, unsigned singular_corner)
    {
        ret_quad_type result;
        return duffy_impl(
            degree,
            duffy_traits<lset_t>::duffy_corner_indices[singular_corner],
            lset_t::corner_at(singular_corner),
            result);
    }
 
    static ret_quad_type on_face(unsigned degree, xi_t const &singular_point)
    {
        static unsigned const n = domain_t::num_corners;
        unsigned face_indices[n+2];
        face_indices[0] = n;
        for (unsigned c = 0; c <= n; ++c)
            face_indices[1+c] = c % n;
        ret_quad_type result;
        return duffy_impl(degree, face_indices, singular_point, result);
    }
}; // end of class duffy_quadrature
 
 
template <>
struct duffy_traits<tria_1_shape_set>
{
    static unsigned const duffy_corner_indices[3][2+1];
};
 
template <>
struct duffy_traits<quad_1_shape_set>
{
    static unsigned const duffy_corner_indices[4][3+1];
    static unsigned const duffy_face_indices[5+1];
};
 
template <>
struct duffy_traits<tria_2_shape_set>
{
    static unsigned const duffy_corner_indices[6][3+1];
};
 
template <>
struct duffy_traits<quad_2_shape_set>
{
    static unsigned const duffy_corner_indices[9][5+1];
};
 
}
 
#endif // DUFFY_QUADRATURE_HPP_INCLUDED