gsOverIntegrateRule.h

Go to the documentation of this file.

#pragma once

#include <gsAssembler/gsQuadRule.h>
#include <gsCore/gsBasis.h>

namespace gismo
{

template<class T>
class gsOverIntegrateRule GISMO_FINAL : public gsQuadRule<T>
{
public:

    typedef memory::unique_ptr<gsOverIntegrateRule> uPtr;

    gsOverIntegrateRule()
    :
    m_basis(nullptr),
    m_interior(),
    m_boundary()
    {}

    gsOverIntegrateRule(const  gsBasis<T> & basis,
                        const  std::vector<gsQuadRule<T> > & quadInterior,
                        const  std::vector<gsQuadRule<T> > & quadBoundary)
    :
    m_basis(&basis),
    m_interior(quadInterior),
    m_boundary(quadBoundary)
    {
        std::vector< gsVector<T> > nodes(m_basis->dim());
        m_start = m_basis->support().col(0);
        m_end = m_basis->support().col(1);
    };

    static uPtr make(   const  gsBasis<T> & basis,
                        const  std::vector<gsQuadRule<T> > & quadInterior,
                        const  std::vector<gsQuadRule<T> > & quadBoundary )
    { return uPtr( new gsOverIntegrateRule(basis,quadInterior,quadBoundary) ); }

    //const unsigned digits = std::numeric_limits<T>::digits10 );

    ~gsOverIntegrateRule() { }

public:
    index_t dim() const { return m_basis->dim(); }

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

        gsVector<T> bot = lower-m_start;
        gsVector<T> top = upper-m_end;
        std::vector<gsVector<T> > elNodes(m_basis->dim());
        std::vector<gsVector<T> > elWeights(m_basis->dim());
        gsMatrix<T> tmp;
        for (index_t d = 0; d!=dim(); d++)
        {
            if (bot[d]==0.0 || top[d]==0.0) // is the coordinate on a side?
            {
                m_boundary[d].mapTo( lower[d], upper[d], tmp, elWeights[d]);
                GISMO_ASSERT(tmp.rows()==1,"Dimension of the nodes is wrong!");
                elNodes[d] = tmp.transpose();
            }
            else
            {
                m_interior[d].mapTo( lower[d], upper[d], tmp, elWeights[d]);
                GISMO_ASSERT(tmp.rows()==1,"Dimension of the nodes is wrong!");
                elNodes[d] = tmp.transpose();
            }
        }
        this->computeTensorProductRule_into(elNodes,elWeights,nodes,weights);
    }

private:
    const gsBasis<T> * m_basis;
    std::vector<gsQuadRule<T> > m_interior, m_boundary;
    gsVector<T> m_start,m_end;

}; // class gsOverIntegrateRule

} // namespace gismo