gsBoundary.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>

#include <gsCore/gsExport.h>

namespace gismo
{
struct boundary
{
    enum side { west  = 1,  east = 2, south = 3, north= 4, front=5, back=6,
                stime = 7, etime = 8,
                left  = 1, right = 2, down  = 3, up   = 4 , none = 0 };


    // warning: naming southwest etc ambiguous for 3D (corresponds to an edge)
    enum corner { southwestfront = 1, southeastfront = 2, northwestfront = 3, northeastfront = 4,
                  southwestback  = 5, southeastback  = 6, northwestback  = 7, northeastback  = 8,
                  southwest      = 1, southeast      = 2, northwest      = 3, northeast      = 4
    };
};

struct interaction
{
    enum type { conforming = 0, contact = 2};
};

// forward declarations; defined later
struct boxCorner;
struct patchCorner;

class GISMO_EXPORT boxSide
{
public:

    short_t m_index;
public:
    boxSide (): m_index(0) {}
    boxSide (short_t dir, bool par) : m_index(index(dir,par))
    { GISMO_ASSERT(dir>=0,"invalid side");}
    boxSide (short_t a) : m_index(a) { GISMO_ASSERT(a>=0,"invalid side");}
    boxSide (boundary::side a) : m_index(a) { GISMO_ASSERT(a>=0,"invalid side");}

    // conversions
    operator boundary::side() const {return static_cast<boundary::side>(m_index);}

    short_t direction () const {return static_cast<short_t>((m_index-1) / 2);}

    bool    parameter () const {return (m_index-1)%2 != 0;}

    boxSide opposite() const {return boxSide(static_cast<short_t>(((m_index-1)^1)+1));}


    short_t    index () const {return m_index;}

    static inline short_t index (short_t dir, bool par) {return static_cast<short_t>(par?2*dir+2:2*dir+1);}


    void getContainedCorners (short_t dim, std::vector<boxCorner> &corners) const;

    static  boxSide    getFirst     (short_t)
    { return boxSide(1); }

    static  boxSide    getLast      (short_t dim) {return boxSide(static_cast<short_t>(2*dim));}
    static  boxSide    getEnd       (short_t dim) {return boxSide(static_cast<short_t>(2*dim+1));}

    boxSide& operator++ ()    { ++m_index; return *this;} //prefix
    //boxSide  operator++ (int) { boxSide temp(*this); ++m_index; return temp;} //postfix

    boxSide& operator-- ()    { --m_index; return *this;} //prefix
    //boxSide  operator-- (int) { boxSide temp(*this); --m_index; return temp;} //postfix

};

inline std::ostream &operator<<(std::ostream &os, const boxSide& o)
{
    switch (o)
    {
    case 0:
        os<<"none ";
        break;
    case 1:
        os<< "west ";
        break;
    case 2:
        os<< "east ";
        break;
    case 3:
        os<< "south";
        break;
    case 4:
        os<< "north";
        break;
    case 5:
        os<< "front";
        break;
    case 6:
        os<< "back ";
        break;
    default:
        os<< "side ";
        break;
    };
    os<<"("<<static_cast<int>(o)<<")";
    return os;
}

struct GISMO_EXPORT patchSide : public boxSide
{
public:
    index_t patch;              
protected:
    std::string m_label;        
public:

    patchSide() : boxSide(), patch(0), m_label("") { }

    patchSide(index_t p, boxSide s, const std::string & l="")
        : boxSide(s), patch(p), m_label(l) { }

    patchSide(index_t p, boundary::side s, const std::string & l="")
        : boxSide(s), patch(p), m_label(l) { }

    // Accessors
    boxSide& side()       {return *this;}
    const boxSide& side() const {return *this;}

    index_t patchIndex() const {return patch;}
    const std::string & label() const {return m_label;}
    std::string & label() {return m_label;}

    void setLabel(std::string label) { m_label = label; };

    void getContainedCorners (short_t dim, std::vector<patchCorner> &corners) const;

    using boxSide::getContainedCorners; // unhiding
};

inline std::ostream &operator<<(std::ostream &os, patchSide const & i)
{
    //os<<"Side: patch="<< i.patch<<", side="<< int(i.side)<<", ie. "<< i.side << ".\n";
    os<<i.patch<<":"<< i.side()<<" "<<i.label();
    return os;
}

// We do not need comparization operators for boxSide since it decays to a number anyway

inline bool operator== (const patchSide& a, const patchSide& b)
{ return a.patch == b.patch && a.m_index == b.m_index; }
inline bool operator<  (const patchSide& a, const patchSide& b)
{ return a.patch < b.patch || (a.patch == b.patch && a.m_index < b.m_index); }
GISMO_DELEGATING_COMPARISON_OPERATORS(patchSide)

// This might cause terrible bugs
GISMO_DELETE_COMPARISON_OPERATORS(boxSide,patchSide)



struct boxCorner
{
public:
    index_t m_index;
public:
    boxCorner (index_t a=0) : m_index(a) { GISMO_ASSERT(a>=0,"invalid side");}
    boxCorner (boundary::corner a) : m_index(a) { GISMO_ASSERT(a>=0,"invalid side");}
    boxCorner (gsVector<bool>  v) : m_index(1)
    {
        for (index_t p=0;p<v.rows();++p)
        {
            if(v(p))
                m_index+= 1<<p;
        }
    }

    // conversions
    operator boundary::corner() const {return static_cast<boundary::corner>(m_index);}
    boxCorner& operator= (boundary::corner a) {m_index= a; return *this;}

    void parameters_into (index_t dim, gsVector<bool> &param) const
    {
        param.resize(dim);
        for (short_t i=0; i<dim; ++i)
            param(i)=((m_index-1)>>i)&1;
    }

    gsVector<bool> parameters(index_t dim) const
    {
        gsVector<bool> r;
        parameters_into(dim,r);
        return r;
    }


    void getContainingSides (index_t dim, std::vector<boxSide> &sides) const
    {
        GISMO_ASSERT(dim>=0, "Dimension must be non negative");
        sides.resize(dim);
        gsVector<bool> param;
        parameters_into(dim,param);
        for (short_t i=0;i<dim;++i)
            sides[i]=boxSide(i,param(i));
    }

    static  boxCorner    getFirst     (short_t)
    { return boxCorner(1); }

    static  boxCorner    getLast      (short_t dim) {return boxCorner((1<<dim));}
    static  boxCorner    getEnd       (short_t dim) {return boxCorner((1<<dim)+1);}
    boxCorner& operator++ () { ++m_index; return *this;} //prefix
    //boxCorner operator++ (int ) { boxCorner temp(*this); ++m_index; return temp;} //postfix

    boxCorner& operator-- () { --m_index; return *this;} //prefix
    //boxCorner operator-- (int ) { boxCorner temp(*this); --m_index; return temp;} //postfix

};


struct patchCorner : public boxCorner
{
public:
    index_t patch;
public:
    patchCorner() : boxCorner(0) { }
    patchCorner(index_t p,boundary::corner c)
        : boxCorner(c), patch (p) { }

    patchCorner(index_t p, int c)
        : boxCorner(c), patch (p) { }

    patchCorner(index_t p, boxCorner c)
        : boxCorner(c), patch (p) { }

    // Accessors
    boxCorner& corner()       {return *this;}
    const boxCorner& corner() const {return *this;}

    void getContainingSides (short_t dim, std::vector<patchSide> &sides) const
    {
        GISMO_ASSERT(dim>=0, "Dimension must be non negative");
        sides.resize(dim);
        gsVector<bool> param;
        parameters_into(dim,param);
        for (short_t i=0;i<dim;++i)
            sides[i]=patchSide(patch, boxSide(i, param(i)));
    }
};

// We do not need comparization operators for boxCorner since it decays to a number anyway

inline bool operator== (const patchCorner& a, const patchCorner& b)
{ return a.patch == b.patch && a.m_index == b.m_index; }
inline bool operator<  (const patchCorner& a, const patchCorner& b)
{ return a.patch < b.patch || (a.patch == b.patch && a.m_index < b.m_index); }
GISMO_DELEGATING_COMPARISON_OPERATORS(patchCorner)

// This might cause terrible bugs
GISMO_DELETE_COMPARISON_OPERATORS(boxCorner,patchCorner)

struct GISMO_EXPORT boxComponent {
private:
    index_t m_index;

    short_t m_total_dim;    
public:

    boxComponent( index_t b, short_t total_dim ) : m_index(b), m_total_dim(total_dim) {}

    boxComponent( short_t total_dim ) : m_index(0), m_total_dim(total_dim) {}

    boxComponent( boxSide b, short_t total_dim );

    boxComponent( boxCorner b, short_t total_dim );

    short_t totalDim() const { return m_total_dim; }

    short_t dim() const;

    index_t index() const { return m_index; }

    std::vector<boxCorner> containedCorners() const;

    std::vector<boxSide> containingSides() const;

    boxSide asSide() const;

    boxCorner asCorner() const;

    boxComponent opposite() const;

    enum location {
        interior = 0,  
        begin = 1,     
        end = 2        
    };

    location locationForDirection(index_t direction) const;

    void setLocationForDirection(index_t direction, location par);

};

struct GISMO_EXPORT patchComponent : boxComponent {
private:
    index_t m_patch;    
public:

    patchComponent( index_t p, index_t b, short_t total_dim ) : boxComponent(b,total_dim), m_patch(p) {}

    patchComponent( index_t p, short_t total_dim ) : boxComponent(total_dim), m_patch(p) {}

    patchComponent( index_t p, boxComponent b ) : boxComponent(b), m_patch(p) {}

    patchComponent( patchSide p, short_t total_dim ) : boxComponent(p,total_dim), m_patch(p.patch) {}

    patchComponent( patchCorner p, short_t total_dim ) : boxComponent(p,total_dim), m_patch(p.patch) {}

    std::vector<patchCorner> containedCorners() const;

    std::vector<patchSide> containingSides() const;

    patchSide asSide() const
    {
        return patchSide( m_patch, boxComponent::asSide() );
    }

    patchCorner asCorner() const
    {
        return patchCorner( m_patch, boxComponent::asCorner() );
    }

    patchComponent opposite() const
    {
        return patchComponent( m_patch, boxComponent::opposite() );
    }

    index_t patch() const { return m_patch; }
};

inline bool operator== (const boxComponent& a, const boxComponent& b)
{ return a.index() == b.index(); }
inline bool operator<  (const boxComponent& a, const boxComponent& b)
{ return a.index() <  b.index(); }
GISMO_DELEGATING_COMPARISON_OPERATORS(boxComponent)

inline bool operator== (const patchComponent& a, const patchComponent& b)
{ return a.patch() == b.patch() && a.index() == b.index(); }
inline bool operator<  (const patchComponent& a, const patchComponent& b)
{ return a.patch() < b.patch() || (a.patch() == b.patch() && a.index() < b.index()); }
GISMO_DELEGATING_COMPARISON_OPERATORS(patchComponent)

// This might cause terrible bugs
GISMO_DELETE_COMPARISON_OPERATORS(boxComponent,patchComponent)

struct GISMO_EXPORT boundaryInterface
{
public:
    boundaryInterface() { }

    // special constructor for the 2d case
    boundaryInterface(patchSide const & _ps1,
                      patchSide const & _ps2,
                      bool o1,
                      const std::string & l = "")
        : m_label(l), ps1(_ps1), ps2(_ps2), m_type(interaction::conforming)
    {
        directionMap.resize(2);
        directionOrientation.resize(2);
        directionMap(ps1.direction())=ps2.direction();
        directionMap(1-ps1.direction())=1-ps2.direction();
        directionOrientation(ps1.direction())= (ps1.parameter()==ps2.parameter());
        directionOrientation(1-ps1.direction())=o1;
    }


    //
    boundaryInterface(patchSide const & _ps1,
                      patchSide const & _ps2,
                      short_t dim,
                      const std::string & l = "")
        : m_label(l), ps1(_ps1), ps2(_ps2), m_type(interaction::conforming)
    {
        directionMap.resize(dim);
        directionOrientation.resize(dim);

        for (int i = 1 ; i < dim; ++i)
        {
            const index_t o = (ps1.direction()+i)%dim;
            const index_t d = (ps2.direction()+i)%dim;

            directionMap(o)=d;
            directionOrientation(o)=true;
        }
        directionMap(ps1.direction())=ps2.direction();
        directionOrientation(ps1.direction())= (ps1.parameter()==ps2.parameter());

    }

    boundaryInterface(gsVector<short_t>     const & p,
                      gsVector<index_t> const & map_info,
                      gsVector<bool>    const & orient_flags,
                      const std::string & l = "")
    : m_label(l), ps1(p(0),p(1)), ps2(p(2),p(3)),
      directionMap(map_info),
      directionOrientation(orient_flags), m_type(interaction::conforming)
    {
        GISMO_ASSERT(p.size() == 4, "Expecting four integers");
        directionMap(ps1.direction())=ps2.direction();
        directionOrientation(ps1.direction())= (ps1.parameter()==ps2.parameter());
    }

    boundaryInterface(patchSide const & _ps1,
                      patchSide const & _ps2,
                      gsVector<index_t> const & map_info,
                      gsVector<bool>    const & orient_flags,
                      const std::string & l = "")
        : m_label(l), ps1(_ps1), ps2(_ps2), directionMap(map_info), directionOrientation(orient_flags)
        , m_type(interaction::conforming)
    {
        directionMap(ps1.direction())=ps2.direction();
        directionOrientation(ps1.direction())= (ps1.parameter()==ps2.parameter());
    }

    GISMO_DEPRECATED boundaryInterface(patchSide const & _ps1,
                      patchSide const & _ps2,
                      gsVector<bool>    const & orient_flags,
                      const std::string & l = "")
    {
        init(_ps1,_ps2,orient_flags,l);
    }

    GISMO_DEPRECATED boundaryInterface(gsVector<short_t>     const & p,
                      gsVector<bool>    const & orient_flags,
                      const std::string & l = "")
    {
        init(patchSide(p(0),boxSide(p(1))),patchSide(p(2),boxSide(p(3))) ,orient_flags,l);
    }

    //DEPRECATED
    void init (patchSide const & _ps1,
               patchSide const & _ps2,
               gsVector<bool>    const & orient_flags,
            const std::string & l = "")
    {
        ps1=_ps1;
        ps2=_ps2;
        m_label = l;

        const index_t dim = orient_flags.cols()+1;
        directionMap.resize(dim);
        directionOrientation.resize(dim);

        directionMap(ps1.direction())=ps2.direction();
        directionOrientation(ps1.direction())= (ps1.parameter()==ps2.parameter());

        directionMap(1-ps1.direction())=1-ps2.direction();
        directionOrientation(1-ps1.direction())= orient_flags(0);
        m_type = interaction::conforming;
    }

    const std::string & label() const {return m_label;}
    std::string & label() {return m_label;}


    bool operator== (const boundaryInterface & other) const
    {
        return ps1==other.ps1 && ps2==other.ps2
                && directionMap==other.directionMap
                && directionOrientation==other.directionOrientation;
    }

    inline bool operator< (const boundaryInterface& other) const
    {
        return ps1<other.ps1 || (ps1==other.ps1 && ps2<other.ps2);
    }

          patchSide& first ()        {return ps1;}
    const patchSide& first () const  {return ps1;}

          patchSide& second ()        {return ps2;}
    const patchSide& second () const  {return ps2;}

          patchSide&   other (const patchSide & ps)        {return ps==ps1 ? ps2 : ps1;}
    const patchSide&   other (const patchSide & ps) const  {return ps==ps1 ? ps2 : ps1;}

    // use boundaryInterface.first() and boundaryInterface.second()
    GISMO_DEPRECATED patchSide  operator [] (size_t i) const
    {
        if (i==0)
            return ps1;
        else if (i==1)
            return ps2;
        else
            GISMO_ERROR("Invalid index "<<i<<": Interface has 2 elements(sides).");
    }
    // use boundaryInterface.first() and boundaryInterface.second()
    GISMO_DEPRECATED patchSide & operator [] (size_t i)
    {
        if (i==0)
            return ps1;
        else if (i==1)
            return ps2;
        else
            GISMO_ERROR("Invalid index "<<i<<": Interface has 2 elements(sides).");
    }


    // this is a work in progress, the old code was 2D specific, look at mapInfo and mapOrientation
    GISMO_DEPRECATED gsVector<bool> orient () const
    {
        GISMO_ASSERT(directionOrientation.size()==2, "This is deprecated and does not work if dim>2");
        gsVector<bool> result(1);
        result(0)=directionOrientation(1-ps1.direction());
        return result;
    }


    boundaryInterface getInverse() const
    {
        boundaryInterface result;
        result.directionMap.resize(directionMap.rows());
        result.directionOrientation.resize(directionOrientation.rows());
        result.ps1=ps2;
        result.ps2=ps1;
        result.m_type = m_type;
        for (index_t i=0;i<directionMap.rows();++i)
        {
            result.directionMap(directionMap(i))=i;
            result.directionOrientation(directionMap(i))=directionOrientation(i);
        }
        return result;
    }

    patchCorner mapCorner ( const patchCorner c) const
    {
        gsVector<bool> par=c.parameters(directionMap.rows());
        const index_t dim=par.rows();
        gsVector<bool> new_par(dim);
        if (c.patch == ps1.patch && par(ps1.direction()) == ps1.parameter() )
        {
            index_t i=0;
            for (; i<ps1.direction();++i)
            {
                new_par(directionMap(i)) = directionOrientation(i) ?
                            par(i) : !par(i);
            }
            new_par(directionMap(i)) = ps2.parameter();
            for (++i; i<dim;++i)
            {
                new_par(directionMap(i)) = directionOrientation(i) ?
                            par(i) : !par(i);
            }
            return patchCorner(ps2.patch, boxCorner(new_par));
        }
        else if (c.patch == ps2.patch && par(ps2.direction()) == ps2.parameter() )
        {
            return getInverse().mapCorner(c);
        }
        else
        {
            gsWarn<<"cannot map corners that are not in the interface";
            return c;
        }
    }

    bool dirOrientation(const patchSide & ps,index_t dir) const
    {
        if(ps==ps1)
            return directionOrientation(dir);
        else
            return getInverse().dirOrientation(ps,dir);
    }

    index_t dirMap(const patchSide & ps,index_t dir) const
    {
        if(ps==ps1)
            return directionMap(dir);
        else
            return getInverse().dirMap(ps,dir);
    }

    gsVector<bool> dirOrientation(const patchSide & ps) const
    {
        if(ps==ps1)
            return directionOrientation;
        else
            return getInverse().dirOrientation(ps);
    }

    gsVector<index_t>  dirMap(const patchSide & ps) const
    {
        if(ps==ps1)
            return directionMap;
        else
            return getInverse().dirMap(ps);
    }

    const gsVector<index_t> & dirMap() const
    { return directionMap; }

    const gsVector<bool> & dirOrientation()  const
    { return directionOrientation; }

    void faceData(gsVector<bool> & flip, gsVector<index_t> & perm) const;

    void cornerMap(gsVector<index_t> & cmap) const;

    void reorderCorners(gsMatrix<index_t> & boundary) const;

    void setAsContact() { m_type = interaction::contact; }

    interaction::type type() const { return m_type; }

private:
    std::string m_label; 

    patchSide ps1; 
    patchSide ps2; 

    gsVector<index_t> directionMap;
    gsVector<bool>    directionOrientation;


    interaction::type m_type;

protected:
    friend std::ostream &operator<<(std::ostream &os, const boundaryInterface & i);
};


inline std::ostream &operator<<(std::ostream &os, const boundaryInterface & i)
{
    os << "interface between "<<i.ps1.patch<<":"<< i.ps1.side()<<" and "
       << i.ps2.patch<<":"<<i.ps2.side()<<" [ ";
    for (index_t j = 0; j<i.directionMap.size(); ++j)
    {
        if ( i.ps1.direction() == j )
            continue;
        os << j << "~" << (i.directionOrientation(j) ? "(+" : "(-") << i.directionMap(j)<<") ";
    }
    os <<"]";
    return os;
}


inline int sideOrientation(short_t s)
{
    // The orientation sign is determined by the formula:
    // (direction + parameter) + 1 mod 2
    // where 0 means "-" (CW) and 1 means "+" (CCW)
    return ( ( s + (s+1) / 2 ) % 2 ? 1 : -1 );
}

GISMO_DEPRECATED inline index_t direction (index_t s)
{
    GISMO_ASSERT( s>0, "Requested direction of none boundary.\n");
    return (s-1) / 2 ;
}


GISMO_DEPRECATED inline bool parameter (int s)
{
    GISMO_ASSERT( s>0, "Requested parameter of none boundary.\n");
    return ( (s+1) % 2 == 0 ? false : true ) ;
}

template <typename T>
gsMatrix<T> getFace (const boxSide side, const gsMatrix<T> &box)
{
    GISMO_ASSERT(side.direction()< box.rows(), "the specified side is not appropriate for the domain dimension");
    gsMatrix<T> temp=box;
    const index_t dir=side.direction();
    if (side.parameter()) // bottom face
        temp(dir,0)=box(dir,1);
    else    // top face
        temp(dir,1)=box(dir,0);
    return temp;
}

#ifdef GISMO_WITH_PYBIND11

  void pybind11_enum_gsBoundary(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo