gsPatchRule.h

Go to the documentation of this file.

#pragma once

#include <gsAssembler/gsQuadRule.h>
#include <gsNurbs/gsKnotVector.h>
#include <gsNurbs/gsBSplineBasis.h>

namespace gismo
{

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

    typedef memory::unique_ptr<gsPatchRule> uPtr;

    gsPatchRule()
    :
    m_basis(nullptr),
    m_deg(0),
    m_reg(0),
    m_over(false),
    m_fixDir(-1)
    {};

    gsPatchRule(const gsBasis<T> & basis,
                const index_t degree,
                const index_t regularity,
                const bool overintegrate,
                const short_t fixDir = -1);


    static uPtr make(   const gsBasis<T> & basis,
                        const index_t degree,
                        const index_t regularity,
                        const bool overintegrate,
                        const short_t fixDir = -1)
    { return uPtr( new gsPatchRule(basis,degree,regularity,overintegrate,fixDir) ); }


    ~gsPatchRule() { };

public:

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

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

    void mapToAll( const std::vector<T> &,
                   gsMatrix<T> &, gsVector<T> &) const
    { GISMO_NO_IMPLEMENTATION }

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

    gsVector<T> nodes(const index_t d=0) const { return m_nodes.at(d); }
    gsVector<T> weights(const index_t d=0) const { return m_weights.at(d); }

protected:

    index_t _numQuads(const gsKnotVector<T> knots) const
    {
        index_t ndof = (m_deg + 1)*knots.numElements() - (m_reg+1)*(knots.numElements()-1);
        ndof = static_cast<index_t>(std::ceil(ndof/2.0));
        return ndof;
    }

    bool _isSymmetric(const gsKnotVector<T> knots) const
    {
        gsKnotVector<T> copy = knots;
        copy.reverse();
        std::vector<T> result;
        knots.symDifference(copy,result);
        return result.size()==0 ? true : false;
    }


    gsKnotVector<T> _init(const gsBSplineBasis<T> * Bbasis) const;

    std::pair<gsMatrix<T>,gsVector<T> > _integrate(const gsKnotVector<T> & knots ) const;

    std::pair<gsVector<T>,gsVector<T> > _compute(const gsKnotVector<T> & knots,
                                                const gsMatrix<T> & greville,
                                                const gsVector<T> & integrals,
                                                const T tol = 1e-10) const;

private:
    const gsBasis<T> * m_basis;
    const index_t m_deg,m_reg;
    const bool m_over;
    const short_t m_fixDir;

    std::vector<gsVector<T> > m_nodes;
    std::vector<gsVector<T> > m_weights;

    gsVector<T> m_end;

    mutable typename gsSparseSolver<T>::QR m_solver;

    size_t m_dim;

    std::vector<std::map<T,T> > m_maps;
}; // class gsPatchRule

} // namespace gismo


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