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 QUADRATURE_HPP_INCLUDED
#define QUADRATURE_HPP_INCLUDED
 
#include "shapeset.hpp"
#include "../library/lib_shape.hpp"
#include "../util/eigen_utils.hpp"
#include "../util/crtp_base.hpp"
 
namespace NiHu
{
 
template <class XiType, class ScalarType>
class quadrature_elem
{
public:
    typedef XiType xi_t;
    typedef ScalarType scalar_t;
 
    quadrature_elem(xi_t const &xi = xi_t(), scalar_t const &w = scalar_t())
        : m_xi(xi), m_w(w)
    {
    }
 
    xi_t &get_xi(void)
    {
        return m_xi;
    }
    
    xi_t const &get_xi(void) const
    {
        return m_xi;
    }
 
    scalar_t const &get_w(void) const
    {
        return m_w;
    }
 
    void set_w(scalar_t const &w)
    {
        m_w = w;
    }
 
    void set_xi(xi_t const& xi)
    {
        m_xi = xi;
    }
 
    quadrature_elem &operator *=(scalar_t const &c)
    {
        m_w *= c;
        return *this;
    }
 
    template <class LSet, bool Signed = true>
    quadrature_elem &transform_inplace(
        Eigen::Matrix<scalar_t, LSet::num_nodes, LSet::domain_t::dimension> const &coords)
    {
        // CRTP check
        static_assert(std::is_base_of<shape_set_base<LSet>, LSet>::value,
            "LSet must be derived from shape_set_base<LSet>");
        // compute Jacobian and multiply with original weight
        auto j = (LSet::template eval_shape<1>(m_xi).transpose() * coords).determinant();
        if (Signed)
            m_w *= j;
        else
            m_w *= std::abs(j);
        // compute new quadrature location
        m_xi = LSet::template eval_shape<0>(m_xi).transpose() * coords;
        return *this;
    }
 
    // struct is used in a std::vector, therefore this declaration is necessary
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 
protected:
    xi_t m_xi;
    scalar_t m_w;
};
 
 
// forward declaration of quadrature traits
template <class Derived>
struct quadrature_traits;
 
template <class Derived>
struct quadr_elem
{
    typedef quadrature_elem<
        typename quadrature_traits<Derived>::domain_t::xi_t,
        typename quadrature_traits<Derived>::domain_t::scalar_t
    > type;
};
 
template <class Derived>
class quadrature_base :
    public eigen_std_vector<typename quadr_elem<Derived>::type>::type
{
public:
    typedef typename eigen_std_vector<typename quadr_elem<Derived>::type>::type base_t;
 
    NIHU_CRTP_HELPERS
 
public:
    typedef quadrature_traits<Derived> traits_t;
    typedef typename traits_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
    typedef typename domain_t::scalar_t scalar_t;
    typedef quadrature_elem<xi_t, scalar_t> quadrature_elem_t;
 
    quadrature_base(size_t N = 0)
    {
        base_t::reserve(N);
    }
 
    scalar_t sum_of_weights(void) const
    {
        scalar_t res = 0.0;
        for (auto it = base_t::begin(); it != base_t::end(); ++it)
            res += it->get_w();
        return res;
    }
 
    std::ostream & print(std::ostream & os) const
    {
        os << base_t::size() << std::endl;
        for (auto it = base_t::begin(); it != base_t::end(); ++it)
            os << it->get_xi().transpose() << "\t\t" << it->get_w() << std::endl;
        return os;
    }
 
    Derived &operator *=(scalar_t const &c)
    {
        for (size_t i = 0; i < base_t::size(); ++i)
            (*this)[i] *= c;
        return derived();
    }
 
    template <class LSet, bool Signed = true>
    Derived transform(
        Eigen::Matrix<scalar_t, LSet::num_nodes, LSet::domain_t::dimension> const &coords) const
    {
        // CRTP check
        static_assert(std::is_base_of<shape_set_base<LSet>, LSet>::value,
            "LSet must be derived from shape_set_base<LSet>");
        // dimension check
        static_assert(std::is_same<xi_t, typename LSet::xi_t>::value,
            "Quadrature and shape set dimensions must match");
        Derived result(derived());
        return result.template transform_inplace<LSet, Signed>(coords);
    }
 
    template <class LSet, bool Signed = true>
    Derived &transform_inplace(
        const Eigen::Matrix<scalar_t, LSet::num_nodes, LSet::domain_t::dimension> &coords)
    {
        // CRTP check
        static_assert(std::is_base_of<shape_set_base<LSet>, LSet>::value,
            "LSet must be derived from shape_set_base<LSet>");
        // dimension check
        static_assert(std::is_same<xi_t, typename LSet::xi_t>::value,
            "Quadrature and shape set dimensions must match");
        for (auto it = base_t::begin(); it != base_t::end(); ++it)
            it->template transform_inplace<LSet, Signed>(coords);
        return derived();
    }
 
    template <class otherDerived>
    Derived operator +(const quadrature_base<otherDerived> &other) const
    {
        static_assert(std::is_same<xi_t, typename otherDerived::xi_t>::value,
            "Quadrature domain dimensions must match");
        Derived result(base_t::size()+other.size());
        result.insert(result.end(), base_t::begin(), base_t::end());
        result.insert(result.end(), other.begin(), other.end());
        return result;
    }
 
    template <class otherDerived>
    Derived &operator +=(const quadrature_base<otherDerived> &other)
    {
        static_assert(std::is_same<xi_t, typename otherDerived::xi_t>::value,
            "Quadrature domain dimensions must match");
        base_t::reserve(base_t::size()+other.size());
        base_t::insert(base_t::end(), other.begin(), other.end());
        return derived();
    }
}; // class quadrature_base
 
template<class Derived>
std::ostream & operator << (std::ostream & os, const quadrature_base<Derived>& Q)
{
    return Q.print(os);
}
 
template <class Family, class Domain>
struct quadrature_type;
 
}
 
 
#endif