gsPlanarDomain.h

Go to the documentation of this file.

#pragma once

#include <iostream>
#include <deque>

#include <gsModeling/gsCurveLoop.h> 
#include <gsModeling/gsTemplate.h>

// #include <gsUtils/gsSortedVector.h>


namespace gismo
{
template <class T>
class gsBemSolution;

// TODO Add gsWriteParaview using PolygonalData....
template<class T>
class gsPlanarDomain
{
public:

    typedef memory::shared_ptr<gsPlanarDomain> Ptr;
    typedef memory::unique_ptr<gsPlanarDomain> uPtr;
    
    gsPlanarDomain() { }

    gsPlanarDomain( gsCurveLoop<T> * boundary)
    {
        if ( boundary->is_ccw() )
        {
            m_loops.push_back(boundary);
        }
        else
        {
            boundary->reverse();
            m_loops.push_back(boundary);
        }

        updateBoundingBox();
    }

    gsPlanarDomain( std::vector< gsCurveLoop<T> *> const & loops);

    gsPlanarDomain( gsCurve<T> * boundary)
    {
        m_loops.push_back( new gsCurveLoop<T>(boundary) );
        updateBoundingBox();
    }

    ~gsPlanarDomain()
    {
        freeAll( m_loops );
    }

    gsPlanarDomain( const gsPlanarDomain & other)
        : m_loops( other.m_loops.size() )
    {
        m_bbox = other.m_bbox;
        cloneAll( other.m_loops.begin(), other.m_loops.end(),
                  this->m_loops.begin() );
    }


    gsPlanarDomain& operator=( const gsPlanarDomain& other)
    {
        freeAll( m_loops );
        m_loops.resize( other.m_loops.size() );
        m_bbox = other.m_bbox;
        cloneAll( other.m_loops.begin(), other.m_loops.end(),
                  this->m_loops.begin() );
        return *this;
    }

    //GISMO_CLONE_FUNCTION(gsPlanarDomain)
    uPtr clone() const
    {
        return uPtr(new gsPlanarDomain(*this));
    }

public:

    void check()
    {
        if ( !m_loops[0]->is_ccw() )
            gsWarn<< "Wrong orientation in outer loop of planar domain.";
        for(size_t i=1; i< m_loops.size(); i++)
        {
            if( m_loops[i]->is_ccw())
                gsWarn<< "Wrong orientation in loop["<< i <<"] of planar domain.";
        }
    }

    void insertHole( gsCurveLoop<T> * hole )
    {
        if ( hole->is_ccw() )
        {
            hole->reverse();
        }
        m_loops.push_back( hole );
    }

    int numLoops() const { return m_loops.size();    }
    int numHoles() const { return m_loops.size() -1; }

    gsCurveLoop<T> & outer()             { return loop(0); }
    const gsCurveLoop<T> & outer() const { return loop(0); }

    gsCurveLoop<T> & loop(unsigned loopNumber)
    {
        GISMO_ASSERT( loopNumber<m_loops.size(), "Loop does not exist" );
        return *m_loops[loopNumber];
    }
    const gsCurveLoop<T> & loop(unsigned loopNumber) const
    {
        GISMO_ASSERT( loopNumber<m_loops.size(), "Loop does not exist" );
        return *m_loops[loopNumber];
    }

    gsCurve<T> & curve(unsigned loopNumber, unsigned curveNumber)
    {
        GISMO_ASSERT( loopNumber<m_loops.size(), "Loop does not exist" );
        return m_loops[loopNumber]->curve(curveNumber);
    }
    const gsCurve<T> & curve(unsigned loopNumber, unsigned curveNumber) const
    {
        GISMO_ASSERT( loopNumber<m_loops.size(), "Loop does not exist" );
        return m_loops[loopNumber]->curve(curveNumber);
    }

    bool contains( gsVector<T> const &, T)
    {
        GISMO_NO_IMPLEMENTATION
    }

    gsMatrix<T> boundingBox() const
    {
        gsMatrix<T> res(2,2);

        res = m_loops[0]->getBoundingBox();
        return res;
    }

    void translate(gsVector<T> const & v)
    {
        for ( typename std::vector< gsCurveLoop<T> *>::iterator it =
              m_loops.begin();  it != m_loops.end(); ++it)
            (*it)->translate(v);
    }

    //gsMatrix<T> averageValue( std::vector<gsFunction<T>*> const &f, std::vector<T> const & breaks);

    bool inDomain( gsMatrix<T> const & u, int direction = 0);

    bool onBoundary(gsMatrix<T> const & u);

    std::ostream &print(std::ostream &os) const;

    friend std::ostream& operator<<(std::ostream& os, const gsPlanarDomain& pd) { return pd.print(os); }

    void sampleLoop_into( int loopID, int npoints, int numEndPoints, gsMatrix<T> & u );

    gsMatrix<T> sampleLoop( int loopID, int npoints=50, int numEndPoints=2)
    {
        gsMatrix<T> u;
        sampleLoop_into( loopID, npoints, numEndPoints, u );
        return u;
    }

    void sampleCurve_into( int loopID, int curveID, int npoints, gsMatrix<T> & u );

    gsMatrix<T> sampleCurve(int loopID, int curveID, int npoints = 50)
    {
        gsMatrix<T> u;
        sampleCurve_into( loopID, curveID, npoints, u );
        return u;
    }

    memory::unique_ptr<gsMesh<T> > toMesh(int npoints = 50) const;     // FOR NOW ONLY ONE LOOP

    uPtr split(int startIndex, int endIndex,
                          gsCurve<T> * newCurveThisFace, gsCurve<T> * newCurveNewFace)
    {
        typename gsCurveLoop<T>::uPtr newCurveLoop = this->m_loops[0]->split(startIndex, endIndex, newCurveThisFace, newCurveNewFace);
        updateBoundingBox();

        return uPtr(new gsPlanarDomain<T>(newCurveLoop.release()));
    }

    // TODO: make this function private or protected. In order to do this we
    // would need to replace the functions "outer", "loop" and "curve" with
    // new public functions for modifying the loops that call this function
    // when they are done.
    void updateBoundingBox()
    {
        assert(!m_loops.empty()); // outer loop does not exist
        m_bbox = m_loops[0]->getBoundingBox();
    }

    gsMatrix<T> splitCurve(size_t loopId, size_t curveId, T lengthRatio=.5)
    {
        return m_loops[loopId]->splitCurve(curveId,lengthRatio);
    }



    // Data members
private:

    // m_loops[0] is CCW, all others CW holes
    std::vector< gsCurveLoop<T> *> m_loops;

    // The lower left and upper right corner of a square bounding the
    // domain
    gsMatrix<T,2,2> m_bbox;

public:
#   define Eigen gsEigen
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
#   undef Eigen

}; // class gsPlanarDomain

} // namespace gismo


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