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
{

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);

    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 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++ )
            {
                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; }
    
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);

}; // 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