gsMultiPatch.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsBoxTopology.h>
#include <gsCore/gsGeometry.h>
 
#include <gsCore/gsMultiBasis.h>
 
namespace gismo
{
 
namespace internal
{
 
struct ElementBlock
{
    index_t numElements;    // NE
    std::list<gsMatrix<index_t>> actives;       // Nodes ( size = NE )
    std::list<gsMatrix<real_t>>  coefVectors;   // Bezier Coefficient Vectors ID ( size = NE )
    index_t PR; // Polynomial degree in R direction  
    index_t PS; // Polynomial degree in S direction  
    index_t PT; // Polynomial degree in T direction  
};
} // end namespace internal
 
 
template<class T>
class gsMultiPatch : public gsBoxTopology, public gsFunctionSet<T>
{
 
public:
    typedef gsBoxTopology    BaseA;
    typedef gsFunctionSet<T> BaseB;
 
    typedef memory::shared_ptr< gsMultiPatch > Ptr;
    typedef memory::unique_ptr< gsMultiPatch > uPtr;
 
    typedef std::vector<gsGeometry<T> *> PatchContainer;
    typedef std::map<boundaryInterface,typename gsGeometry<T>::uPtr>  InterfaceRep;
    typedef std::map<patchSide,typename gsGeometry<T>::uPtr>         BoundaryRep;
 
    typedef typename PatchContainer::iterator iterator;
    typedef typename PatchContainer::const_iterator const_iterator;
 
public:
 
    gsMultiPatch() { }
 
    gsMultiPatch( const gsMultiPatch& other );
 
#if ! EIGEN_HAS_RVALUE_REFERENCES
 
    gsMultiPatch& operator= ( gsMultiPatch other )
    {
        this->swap(other);
        return *this;
    }
 
#else
    gsMultiPatch( gsMultiPatch&& other )
        : BaseA( give(other) ), m_patches( give(other.m_patches) )
    {}
 
    gsMultiPatch& operator= ( const gsMultiPatch& other );
 
    gsMultiPatch& operator= ( gsMultiPatch&& other );
 
#endif
 
    explicit gsMultiPatch( PatchContainer & patches );
 
    gsMultiPatch( const gsGeometry<T> & geo );
 
    gsMultiPatch( PatchContainer & patches,
                  const std::vector<patchSide>& boundary,
                  const std::vector<boundaryInterface>& interfaces );
 
    ~gsMultiPatch();
 
    GISMO_CLONE_FUNCTION(gsMultiPatch)
 
public:
 
    const_iterator begin() const
    { return m_patches.begin(); }
 
    const_iterator end() const
    { return m_patches.end(); }
 
    iterator begin()
    { return m_patches.begin(); }
 
    iterator end()
    { return m_patches.end(); }
 
public:
 
    const gsGeometry<T> & piece(const index_t i) const { return patch(i); }
 
    gsMultiPatch<T> coord(const index_t c) const;
 
    index_t nPieces() const { return static_cast<index_t>(m_patches.size()); }
 
    index_t size() const { return 1; }
 
    index_t coefsSize() const
    {
        index_t sz = 0;
        for (typename PatchContainer::const_iterator it =
                 m_patches.begin(); it != m_patches.end(); ++it )
            sz += (*it)->coefsSize();
        return sz;
    }
 
    gsMatrix<T> coefs() const
    {
        gsMatrix<T> result(this->coefsSize(),this->geoDim());
        result.setZero();
        index_t offset = 0;
        for (typename PatchContainer::const_iterator it =
                 m_patches.begin(); it != m_patches.end(); ++it )
        {
            result.block(offset,0,(*it)->coefsSize(),(*it)->geoDim()) = (*it)->coefs();
            offset += (*it)->coefsSize();
        }
        return result;
    }
 
public:
    gsAffineFunction<T> getMapForInterface(const boundaryInterface &bi, T scaling=0) const;
 
    void swap(gsMultiPatch& other)
    {
        BaseA::swap( other );
        m_patches.swap( other.m_patches );
    }
 
    std::ostream& print( std::ostream& os ) const;
 
    std::string detail() const;
 
    short_t parDim() const
    {
        //GISMO_ASSERT( m_patches.size() > 0 , "Empty multipatch object.");
        return m_dim;
    }
    short_t domainDim () const {return parDim();}
 
    short_t geoDim() const;
    short_t targetDim () const {return geoDim();}
 
    short_t coDim() const;
 
    bool isClosed() { return this->nBoundary() == 0; }
 
    bool empty() const { return m_patches.empty(); }
 
    gsMatrix<T> parameterRange(int i = 0) const;
 
    size_t nPatches() const { return m_patches.size(); }
 
    PatchContainer const& patches() const { return m_patches; }
 
    const BaseA & topology() const { return *this; }
 
    const gsGeometry<T> & operator []( size_t i ) const { return *m_patches[i]; }
 
    std::vector<gsBasis<T> *> basesCopy(bool numeratorOnly = false) const;
 
    gsGeometry<T>& patch( size_t i ) const
    {
        GISMO_ASSERT( i < m_patches.size(), "Invalid patch index "<<i<<" requested from gsMultiPatch" );
        return *m_patches[i];
    }
 
    void permute(const std::vector<short_t> & perm);
 
    gsBasis<T> & basis( const size_t i ) const;
 
    index_t addPatch(typename gsGeometry<T>::uPtr g);
 
    index_t addPatch(const gsGeometry<T> & g);
 
    size_t findPatchIndex( gsGeometry<T>* g ) const;
 
    GISMO_DEPRECATED void addInterface( gsGeometry<T>* g1, boxSide s1,
            gsGeometry<T>* g2, boxSide s2 );
 
    using BaseA::addInterface; // unhide base function
 
    void addPatchBoundary( gsGeometry<T>* g, boxSide s ) {
        size_t p = findPatchIndex( g );
        BaseA::addBoundary( patchSide( p, s ) );
    }
 
    gsMatrix<T> pointOn( const patchCorner& pc );
 
    gsMatrix<T> pointOn( const patchSide& ps );
 
    void uniformRefine(int numKnots = 1, int mul = 1, short_t const dir = -1);
 
    void degreeElevate(short_t const elevationSteps = 1, short_t const dir = -1);
    void degreeIncrease(short_t const elevationSteps = 1, short_t const dir = -1);
 
    void degreeReduce(int elevationSteps = 1);
 
    void degreeDecrease(int elevationSteps = 1);
 
    void uniformCoarsen(int numKnots = 1);
 
    void embed(const index_t N)
    {
        for ( typename PatchContainer::const_iterator it = m_patches.begin();
              it != m_patches.end(); ++it )
        {
            ( *it )->embed(N);
        }
    }
 
    bool computeTopology( T tol = 1e-4, bool cornersOnly = false, bool tjunctions = false);
 
    void fixOrientation();
 
    void closeGaps( T tol = 1e-4 );
 
    gsDofMapper getMapper(T tol) const;
 
    gsSurfMesh toMesh() const;
    
    void clear()
    {
        BaseA::clearAll();
        freeAll(m_patches);
        m_patches.clear();
    }
 
    void boundingBox(gsMatrix<T> & result) const;
 
    gsMultiPatch<T> uniformSplit(index_t dir =-1) const;
 
 
    void repairInterfaces()
    {
        std::vector< boundaryInterface > bivec = interfaces();
        size_t kmax = 2*bivec.size();
        size_t k = 0;
        bool sthChanged = false;
        bool change = false;
        do
        {
            sthChanged = false;
            for( size_t i = 0; i < bivec.size(); i++ )
            {
                if ( bivec[i].type() != interaction::contact)
                {
                    change = repairInterface( bivec[i] );
                    sthChanged = sthChanged || change;
                }
            }
            k++; // just to be sure this cannot go on infinitely
        }
        while( sthChanged && k <= kmax );
    }
 
    bool repairInterface( const boundaryInterface & bi );
 
    gsMultiPatch<T> approximateLinearly(index_t nsamples) const;
 
    void locatePoints(const gsMatrix<T> & points, gsVector<index_t> & pids, gsMatrix<T> & preim, const T accuracy = 1e-6) const;
 
    void locatePoints(const gsMatrix<T> & points, index_t pid1, gsVector<index_t> & pid2, gsMatrix<T> & preim) const;
 
    T closestDistance(const gsVector<T> & pt,std::pair<index_t,gsVector<T> > & result,
                                                   const T accuracy = 1e-6) const;
 
    std::pair<index_t,gsVector<T> > closestPointTo(const gsVector<T> & pt,
                                                   const T accuracy = 1e-6) const;
 
    std::vector<T> HausdorffDistance(   const gsMultiPatch<T> & other,
                                        const index_t nsamples = 1000,
                                        const T accuracy = 1e-6,
                                        const bool directed=false);
 
    T averageHausdorffDistance(         const gsMultiPatch<T> & other,
                                        const index_t nsamples = 1000,
                                        const T accuracy = 1e-6,
                                        const bool directed=false);
 
    void constructInterfaceRep();
    void constructBoundaryRep();
    void constructSides();
 
    void constructInterfaceRep(const std::string l);
    void constructBoundaryRep(const std::string l);
 
    const InterfaceRep & interfaceRep() const { return m_ifaces; }
    const BoundaryRep & boundaryRep() const { return m_bdr; }
    const BoundaryRep & sides() const { return m_sides; }
 
    gsMultiPatch<T> extractBezier() const;
    
protected:
 
    void setIds();
 
    // Data members
private:
 
    PatchContainer m_patches;
 
    InterfaceRep m_ifaces;
    BoundaryRep m_bdr, m_sides;
 
private:
    // implementation functions
 
    // match the vertices in ci1 starting from start to the end with the vertices
    // in ci2 that are still non matched
    // cc1 and cc2 are the physical coordinates of the vertices
    // ci1 and ci2 are the indices corner indices
    // start is the index in ci1 of the vertex to match
    // reference is the image of the vertex ci1(0) and it is used to compute orientation
    // tol is the allowed distance between two vertexes in physical domain
    // matched keeps track of the already matched vertices
    // dirMap and dirO are the output orientation of the match
    // return true if all the vertices starting from start are matched
    static bool matchVerticesOnSide (
           const gsMatrix<T> &cc1, const std::vector<boxCorner> &ci1, index_t start,
           const gsMatrix<T> &cc2, const std::vector<boxCorner> &ci2,
           const gsVector<bool> &matched,
           gsVector<index_t> &dirMap, gsVector<bool>    &dirO,
           T tol, index_t reference=0);
 
public:
    std::map< std::array<size_t, 4>, internal::ElementBlock> BezierOperator() const;
 
}; // class gsMultiPatch
 
 
template<class T>
std::ostream& operator<<( std::ostream& os, const gsMultiPatch<T>& b )
{
    return b.print( os );
}
 
#ifdef GISMO_WITH_PYBIND11
 
  void pybind11_init_gsMultiPatch(pybind11::module &m);
 
#endif // GISMO_WITH_PYBIND11
 
 
} // namespace gismo
 
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsMultiPatch.hpp)
#endif