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 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 index_t _indexFromSides(index_t index1, const patchSide side1, index_t index2, const patchSide side2);
 
 
    virtual index_t _indexFromVert(const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
    virtual index_t _indexFromVert(const gsMultiBasis<T> & bases, const index_t index, const patchCorner corner, const patchSide side, const index_t offsets = 0) const;
    virtual 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 index_t _getValence( patchCorner corner) const
    { return this->_vertexData(corner).first; }
 
    virtual bool _isInteriorVertex( patchCorner corner) const
    { return this->_vertexData(corner).second; }
 
    virtual index_t _sideIndex( index_t patch,  boxSide bside)     const
    { return 4*patch + bside - 1; }
    virtual index_t _sideIndex( patchSide pside)     const
    { return _sideIndex( pside.patch , pside.side() ); }
 
    virtual index_t _vertIndex( index_t patch,  boxCorner corner)  const
    { return 4*patch + corner -1; }
 
    virtual 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