gsQuadRule.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>

namespace gismo
{

template<class T>
class gsQuadRule
{
public:

    typedef memory::unique_ptr<gsQuadRule> uPtr;

    gsQuadRule()
    { }

    virtual ~gsQuadRule() { }

    virtual void setNodes( gsVector<index_t> const &,
                           unsigned = 0)
    { GISMO_NO_IMPLEMENTATION }

    void setNodes( index_t numNodes,
                   unsigned digits = 0 )
    {
        gsVector<index_t> nn(1);
        nn[0] = numNodes;
        this->setNodes(nn, digits);
    }

    const gsMatrix<T> & referenceNodes() const { return m_nodes; }

    const gsVector<T> & referenceWeights() const { return m_weights; }

    // Reference element is [-1,1]^d
    //const gsMatrix<T> & referenceElement() { }

    index_t numNodes() const { return m_weights.size(); }

    index_t dim() const { return m_nodes.rows(); }


    virtual inline void mapTo( const gsVector<T>& lower, const gsVector<T>& upper,
                       gsMatrix<T> & nodes, gsVector<T> & weights ) const;

    void mapTo(const gsMatrix<T>& ab, gsMatrix<T> & nodes) const;
    
    virtual void mapTo( T startVal, T endVal, //mapToInterval
                gsMatrix<T> & nodes, gsVector<T> & weights ) const;

    void mapToAll( const std::vector<T> & breaks,
                   gsMatrix<T> & nodes, gsVector<T> & weights ) const;

protected:

    void computeTensorProductRule(const std::vector<gsVector<T> > & nodes,
                                  const std::vector<gsVector<T> > & weights);

    void computeTensorProductRule_into( const std::vector<gsVector<T> > & nodes,
                                        const std::vector<gsVector<T> > & weights,
                                        gsMatrix<T> & targetNodes,
                                        gsVector<T> & targetWeights
                                        ) const;

protected:

    gsMatrix<T> m_nodes;

    gsVector<T> m_weights;

}; // class gsQuadRule


// Note: left here for inlining
template<class T> void
gsQuadRule<T>::mapTo( const gsVector<T>& lower, const gsVector<T>& upper,
                      gsMatrix<T> & nodes, gsVector<T> & weights ) const
{
    const index_t d = lower.size();
    GISMO_ASSERT( d == m_nodes.rows(), "Inconsistent quadrature mapping");

    nodes.resize( m_nodes.rows(), m_nodes.cols() );
    weights.resize( m_weights.size() );
    nodes.setZero();
    weights.setZero();

    const gsVector<T> h = (upper-lower) / (T)(2) ;
    // Linear map from [-1,1]^d to [lower,upper]
    nodes.noalias() = ( h.asDiagonal() * (m_nodes.array()+1).matrix() ).colwise() + lower;

    T hprod(1.0); //volume of the cube.
    for ( index_t i = 0; i!=d; ++i)
    {
        // the factor 0.5 is due to the reference interval is [-1,1].
        hprod *= ( 0 == h[i] ? (T)(0.5) : h[i] );
    }

    // Adjust the weights (multiply by the Jacobian of the linear map)
    weights.noalias() = hprod * m_weights;
}

} // namespace gismo

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsQuadRule.hpp)
#endif