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