gsDPatchBase.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsBoxTopology.h>
#include <gsCore/gsMultiPatch.h>

#include <gsIO/gsOptionList.h>

#include <gsMSplines/gsMappedBasis.h>
#include <gsMSplines/gsMappedSpline.h>

// #include <gsUnstructuredSplines/src/gsDPatchBase.hpp>

namespace gismo
{


template<short_t d,class T>
class gsDPatchBase
{
protected:
    typedef typename std::vector<std::tuple<index_t,index_t,T>> sparseEntry_t;

public:

    typedef memory::shared_ptr< gsDPatchBase > Ptr;

    typedef memory::unique_ptr< gsDPatchBase > uPtr;

    gsDPatchBase()
    :
    gsDPatchBase(gsMultiBasis<T>(),gsMultiPatch<T>())
    { }

    gsDPatchBase(const gsMultiBasis<T> & mb, const gsMultiPatch<T> & mp)
    :
    m_patches(mp),
    m_Bbases(mb),
    m_topology(m_Bbases.topology()),
    m_computed(false)             
    {
    }

    gsDPatchBase(const gsMultiBasis<T> & mb)
    :
    gsDPatchBase(mb, gsMultiPatch<T>())
    {
    }

    gsDPatchBase(const gsMultiPatch<T> & mp)
    :
    gsDPatchBase(give(gsMultiBasis<T>(mp)), mp)
    {
    }

public:
//----------------------------------------------------------------------------------------------------------------------------
    virtual void compute();

    virtual void defaultOptions();
    // {
    //     GISMO_NO_IMPLEMENTATION;
    // }

    virtual gsOptionList & options() { return m_options; }


//----------------------------------------------------------------------------------------------------------------------------
// Getters for basis, geometry and map
    virtual const gsMultiBasis<T> & localBasis() const
    {
        return m_bases;
    }

    virtual void localBasis_into(gsMultiBasis<T> & localBasis) const
    {
        localBasis = m_bases;
    }

    virtual void localGeometry_into(gsMultiPatch<T> & localGeometry)
    {
        localGeometry = exportToPatches(m_patches);
    }

    virtual const gsMultiPatch<T> & getGeometry() const
    {
        return m_patches;
    }

    virtual void globalBasis_into(gsMappedBasis<d,T> & mbasis) const
    {
        GISMO_ASSERT(m_computed,"The method has not been computed! Call compute().");
        gsSparseMatrix<T> matrix = m_matrix.transpose();
        mbasis.init(m_bases,matrix);
        gsDebugVar("hi");
    }

    virtual void globalGeometry_into(const gsMultiPatch<T> & patches, gsMappedSpline<d,T> & mspline)
    {
        GISMO_ASSERT(!patches.empty() && patches.nPatches()==m_bases.nBases(),"The reference multipatch is empty!");
        gsMappedBasis<d,T> mbasis;
        this->globalBasis_into(mbasis);
        gsMatrix<T> localCoefs = this->_preCoefficients(patches);
        mspline.init(mbasis,localCoefs);
    }

    virtual void globalGeometry_into(gsMappedSpline<d,T> & mspline)
    {
        this->globalGeometry_into(m_patches,mspline);
    }

    virtual void update( gsMappedBasis<d,T> & mbasis ) const
    {
        mbasis.init(m_bases,m_matrix.transpose());
    }

    virtual gsGeometry<T>* exportPatch(index_t patch, bool computeCoefs=true);

    virtual gsMultiPatch<T> exportToPatches(const gsMultiPatch<T> & patches)
    {
        GISMO_ASSERT(m_computed,"The method has not been computed! Call compute().");
        GISMO_ASSERT(!patches.empty(),"The reference multipatch is empty!");
        m_coefs = this->_preCoefficients(patches);
        m_coefs = m_matrix.transpose() * m_coefs;

        std::vector<gsGeometry<T> *> PatchContainer(patches.nPatches());
        for (size_t p=0; p!=patches.nPatches(); p++)
            PatchContainer[p]= this->exportPatch(p,false);

        return gsMultiPatch<T>(PatchContainer,m_topology.boundaries(),m_topology.interfaces());
    }

    virtual gsMultiPatch<T> exportToPatches()
    {
        GISMO_ASSERT(!m_patches.empty(),"The reference multipatch is empty!");
        return this->exportToPatches(m_patches);
    }


    virtual const void matrix_into(gsSparseMatrix<T> & matrix) const
    {
        matrix = this->matrix();
    }

    virtual const gsSparseMatrix<T> & Tmatrix() const
    {
        return m_tMatrix;
    }

    virtual const gsSparseMatrix<T> & matrix() const
    {
        GISMO_ASSERT(m_computed,"The method has not been computed! Call compute().");
        return m_matrix;
    }

//----------------------------------------------------------------------------------------------------------------------------
// Info functions
    virtual void mapperInfo() const;

    virtual void vertexInfo(patchCorner corner) const;

    virtual void sideInfo(patchSide side) const;

    virtual void sideInfo() const;

    virtual void cornerInfo() const;


    gsMatrix<T> preCoefficients() { return _preCoefficients();};

protected:
//----------------------------------------------------------------------------------------------------------------------------
// Pure virtual functions (to be overloaded)
    virtual gsMatrix<T> _preCoefficients(const gsMultiPatch<T> & patches) = 0;
    virtual gsMatrix<T> _preCoefficients()
    {
        return _preCoefficients(m_patches);
    }

    virtual void _initialize();


    virtual void _computeMapper();

    virtual void _computeInterfaceMapper(boundaryInterface iface);

    virtual void _computeBoundaryMapper(patchSide boundary);

    virtual void _computeVertexMapper(patchCorner pcorner);

    virtual void _computeSmoothMatrix();

    void _push(sparseEntry_t entries)
    {
        index_t rowIdx,colIdx;
        T weight;
        for (typename sparseEntry_t::const_iterator it=entries.begin(); it!=entries.end(); it++)
        {
            std::tie(rowIdx,colIdx,weight) = *it;
            m_matrix(rowIdx,colIdx) = weight;
        }
    }

    void _pushAndCheck(sparseEntry_t entries)
    {
        index_t rowIdx,colIdx;
        T weight;
        for (typename sparseEntry_t::const_iterator it=entries.begin(); it!=entries.end(); it++)
        {
            std::tie(rowIdx,colIdx,weight) = *it;
            m_matrix(rowIdx,colIdx) = weight;
            m_basisCheck[rowIdx] = true;
        }
    }

    virtual void _handleVertex(patchCorner pcorner);
    // interior vertices
    virtual void _handleInteriorVertex(patchCorner pcorner, index_t valence);

    virtual void _handleInterface(boundaryInterface iface);
    virtual void _handleBoundary(patchSide side);
    virtual void _handleInterior();

protected:
    virtual void _handleRegularCorner(patchCorner pcorner);

    virtual void _handleRegularBoundaryVertexSmooth(patchCorner pcorner, index_t valence);

    virtual void _handleRegularBoundaryVertexNonSmooth(patchCorner pcorner, index_t valence);

    virtual void _handleIrregularBoundaryVertexSmooth(patchCorner pcorner, index_t valence);

    virtual void _handleIrregularBoundaryVertexNonSmooth(patchCorner pcorner, index_t valence);

protected:

    virtual void _makeTHB() = 0;

    virtual void _computeEVs() = 0;

protected:
    // Boundary vertex of valence 1
    // template<bool _boundary, index_t _v> // valence=2
    // virtual typename std::enable_if<  _boundary && _v==1, void>::type
    // SAME
    virtual void _computeMapperRegularCorner_v1(patchCorner pcorner, index_t valence);

    // Boundary vertex of valence 2 with C1 smoothness
    // template<bool _boundary, index_t _v, bool _smooth> // valence=2
    // virtual  typename std::enable_if<  _boundary && _v==2 && _smooth, void>::type
    // SAME
    virtual void _computeMapperRegularBoundaryVertexSmooth_v2(patchCorner pcorner, index_t valence);

    // Boundary vertex of valence 2 with C0 smoothness
    // template<bool _boundary, index_t _v, bool _smooth> // valence=2
    // virtual  typename std::enable_if<  _boundary && _v==2 && (!_smooth), void>::type
    // DIFFERENT
    virtual void _computeMapperRegularBoundaryVertexNonSmooth_v2(patchCorner pcorner, index_t valence);

    // Boundary vertex of valence !(1,2,3) with C1 smoothness
    // template<bool _boundary, index_t _v, bool _smooth>
    // virtual  typename std::enable_if<  _boundary && _v==-1 && _smooth, void>::type
    // DIFFERENT
    virtual void _computeMapperIrregularBoundaryVertexSmooth_v(patchCorner pcorner, index_t valence);

    // Boundary vertex of valence !(1,2,3) with C0 smoothness
    // template<bool _boundary, index_t _v, bool _smooth>
    // virtual  typename std::enable_if<  _boundary && _v==-1 && (!_smooth), void>::type
    // DIFFERENT
    virtual void _computeMapperIrregularBoundaryVertexNonSmooth_v(patchCorner pcorner, index_t valence);

    // Interior vertex
    // template<bool _boundary, index_t _v>
    // virtual  typename std::enable_if<  (!_boundary) && _v==-1, void>::type
    // DIFFERENT
    virtual void _computeMapperInteriorVertex_v(patchCorner pcorner, index_t valence);


protected:
//----------------------------------------------------------------------------------------------------------------------------
// Helper functions
    virtual std::vector<bool> getSharpCorners(T tol = 1e-2) const;

    virtual const index_t _indexFromSides(index_t index1, const patchSide side1, index_t index2, const patchSide side2);


    virtual const index_t _indexFromVert(const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
    virtual const index_t _indexFromVert(const gsMultiBasis<T> & bases, const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
    virtual const index_t _indexFromVert(const gsBasis<T> * basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
private:
    template<class U>
    typename util::enable_if<util::is_same<U, const gsHTensorBasis<d,T> *>::value,const index_t>::type
    _indexFromVert_impl(U basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;

    template<class U>
    typename util::enable_if<util::is_same<U, const gsTensorBSplineBasis<d,T> *>::value,const index_t>::type
    _indexFromVert_impl(U basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
protected:

    virtual const std::pair<index_t,bool> _vertexData(const patchCorner corner) const;

    virtual const index_t _getValence( patchCorner corner) const
    { return this->_vertexData(corner).first; }

    virtual const bool _isInteriorVertex( patchCorner corner) const
    { return this->_vertexData(corner).second; }

    virtual const index_t _sideIndex( index_t patch,  boxSide bside)     const
    { return 4*patch + bside - 1; }
    virtual const index_t _sideIndex( patchSide pside)     const
    { return _sideIndex( pside.patch , pside.side() ); }

    virtual const index_t _vertIndex( index_t patch,  boxCorner corner)  const
    { return 4*patch + corner -1; }

    virtual const index_t _vertIndex( patchCorner pcorner)     const
    { return _vertIndex( pcorner.patch , pcorner.corner() ); }


    virtual void _getLowestCorners(std::vector<patchCorner> & pcorners, index_t n = 3) const;

    virtual void _removeLowestCorners(std::vector<patchCorner> & pcorners, index_t n = 3) const;

    virtual void _getLowestIndices(std::vector<std::pair<index_t,index_t>> & indices, index_t n = 3) const;

    virtual void _removeLowestIndices(std::vector<std::pair<index_t,index_t>> & indices, index_t n = 3) const;


    virtual std::vector<std::pair<index_t,index_t>> _getInterfaceIndices(patchCorner pcorner, index_t depth, const gsMultiBasis<T> & mbasis) const;
    virtual std::vector<std::pair<index_t,index_t>> _getAllInterfaceIndices(patchCorner pcorner, index_t depth, const gsMultiBasis<T> & mbasis) const;


    virtual bool _checkMatrix(const gsSparseMatrix<T> & matrix) const;

// protected:
//     /**
//      * @brief      Prepares the THB basis if needed.
//      *
//      * This function constructs THB refinements on the places where they are needed, i.e. around EVs. It also constructs the transfer matrix (m_tMatrix) forms the transformation between the original B-spline basis and the THB-Spline basis.
//      */
//     void _makeTHB();

//     /**
//      * @brief      Computes D-Patch smoothing
//      *
//      * Given a basis with THB refinement around the EVs, this function computes the D-Patch smoothing
//      */
//     void _computeEVs();

//     /**
//      * @brief      Makes the Pi matrix
//      *
//      * This matrix is used to transform the coefficients of the D-Patch smoothing matrix
//      *
//      * @param[in]  valence  The valence
//      *
//      * @return     Matrix for smoothing around an EV}
//      */
//     gsMatrix<T> _makePi(index_t valence);

    virtual void _initChecks();

    virtual void _initTHB();

    virtual void _initBasis();

    virtual void _initMappers();

    virtual void _countDoFs()
    {
        GISMO_NO_IMPLEMENTATION;
    }

    virtual void _initMatrix();

    virtual void _performChecks(bool basis);

    virtual void _resetChecks(bool basis);

    virtual void _whichHandled();

// protected:
//     bool _checkMatrix(const gsSparseMatrix<T> & matrix) const // ! makes a deep copy (otherwise the contents of m_matrix get destroyed somehow...)
//     {
//         GISMO_ASSERT(matrix.cols()==matrix.outerSize(),"is the matrix ColMajor?");
//         gsVector<T> colSums(matrix.cols());
//         colSums.setZero();
//         for (index_t i = 0; i<matrix.outerSize(); ++i)
//             for (typename gsSparseMatrix<T>::iterator it = matrix.begin(i); it; ++it)
//                 colSums.at(i) += it.value();

//         return (colSums.array() < 1+1e-8).any() && (colSums.array() > 1-1e-8).any();
//     }

protected:
    const gsMultiPatch<T> & m_patches;
    const gsMultiBasis<T> m_Bbases; // reference?
    
    gsMultiBasis<T> m_bases;

    gsBoxTopology m_topology;

    bool m_computed;

    mutable gsSparseMatrix<T> m_tMatrix;
    mutable std::vector<bool> m_sideCheck;
    mutable std::vector<bool> m_vertCheck;
    mutable std::vector<bool> m_basisCheck;
    mutable std::vector<bool> m_C0s;

    mutable gsDofMapper m_mapModified,m_mapOriginal;

    mutable gsSparseMatrix<T> m_matrix;

    mutable size_t m_size;

    mutable gsMatrix<T> m_coefs;

    gsOptionList m_options;

    size_t m_nSides;
    size_t m_nVerts;

};



}

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