gsSurfMesh.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsIO/gsXml.h>

namespace gismo // typedefs
{
typedef gsEigen::Vector<real_t,3> Point;
typedef gsEigen::Vector<float,2> Vec2f;
typedef Point Normal;
typedef Point Color;
typedef Point Texture_coordinate;
}

#define Eigen gsEigen
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(gismo::Point)
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(gismo::Vec2f)
#undef Eigen

#include <gsMesh2/gsProperty.h>


namespace gismo {

//== CLASS DEFINITION =========================================================


class GISMO_EXPORT gsSurfMesh
{
public:
typedef real_t Scalar;
typedef gismo::Point Point;

private:
typedef gsEigen::Vector<Scalar,3> Vec3;

public: //------------------------------------------------------ topology types

    class Base_handle
    {
    public:

        explicit Base_handle(int _idx=-1) : idx_(_idx) {}

        int idx() const { return idx_; }

        void reset() { idx_=-1; }

        bool is_valid() const { return idx_ != -1; }

        bool operator==(const Base_handle& _rhs) const {
            return idx_ == _rhs.idx_;
        }

        bool operator!=(const Base_handle& _rhs) const {
            return idx_ != _rhs.idx_;
        }

        bool operator<(const Base_handle& _rhs) const {
            return idx_ < _rhs.idx_;
        }

    private:
        friend class Vertex_iterator;
        friend class Halfedge_iterator;
        friend class Edge_iterator;
        friend class Face_iterator;
        friend class gsSurfMesh;
        int idx_;
    };


    struct GISMO_EXPORT Vertex : public Base_handle
    {
        explicit Vertex(int _idx=-1) : Base_handle(_idx) {}
        std::ostream& operator<<(std::ostream& os) const { return os << 'v' << idx(); }
    };


    struct GISMO_EXPORT Halfedge : public Base_handle
    {
        explicit Halfedge(int _idx=-1) : Base_handle(_idx) {}
    };


    struct Edge : public Base_handle
    {
        explicit Edge(int _idx=-1) : Base_handle(_idx) {}
    };


    struct Face : public Base_handle
    {
        explicit Face(int _idx=-1) : Base_handle(_idx) {}
    };




public: //-------------------------------------------------- connectivity types

    struct Vertex_connectivity
    {
        Halfedge  halfedge_;
    };


    struct Halfedge_connectivity
    {
        Face      face_;
        Vertex    vertex_;
        Halfedge  next_halfedge_;
        Halfedge  prev_halfedge_;
    };


    struct Face_connectivity
    {
        Halfedge  halfedge_;
    };




public: //------------------------------------------------------ property types

    template <class T> class Vertex_property : public gsProperty<T>
    {
    public:

        explicit Vertex_property() {}
        explicit Vertex_property(gsProperty<T> p) : gsProperty<T>(p) {}

        typename gsProperty<T>::reference operator[](Vertex v)
        {
            return gsProperty<T>::operator[](v.idx());
        }

        typename gsProperty<T>::const_reference operator[](Vertex v) const
        {
            return gsProperty<T>::operator[](v.idx());
        }
    };


    template <class T> class Halfedge_property : public gsProperty<T>
    {
    public:

        explicit Halfedge_property() {}
        explicit Halfedge_property(gsProperty<T> p) : gsProperty<T>(p) {}

        typename gsProperty<T>::reference operator[](Halfedge h)
        {
            return gsProperty<T>::operator[](h.idx());
        }

        typename gsProperty<T>::const_reference operator[](Halfedge h) const
        {
            return gsProperty<T>::operator[](h.idx());
        }
    };


    template <class T> class Edge_property : public gsProperty<T>
    {
    public:

        explicit Edge_property() {}
        explicit Edge_property(gsProperty<T> p) : gsProperty<T>(p) {}

        typename gsProperty<T>::reference operator[](Edge e)
        {
            return gsProperty<T>::operator[](e.idx());
        }

        typename gsProperty<T>::const_reference operator[](Edge e) const
        {
            return gsProperty<T>::operator[](e.idx());
        }
    };


    template <class T> class Face_property : public gsProperty<T>
    {
    public:

        explicit Face_property() {}
        explicit Face_property(gsProperty<T> p) : gsProperty<T>(p) {}

        typename gsProperty<T>::reference operator[](Face f)
        {
            return gsProperty<T>::operator[](f.idx());
        }

        typename gsProperty<T>::const_reference operator[](Face f) const
        {
            return gsProperty<T>::operator[](f.idx());
        }
    };


    template <class T> class Mesh_property : public gsProperty<T>
    {
    public:

        explicit Mesh_property() {}
        explicit Mesh_property(gsProperty<T> p) : gsProperty<T>(p) {}

        typename gsProperty<T>::reference operator[](size_t idx)
        {
            return gsProperty<T>::operator[](idx);
        }

        typename gsProperty<T>::const_reference operator[](size_t idx) const
        {
            return gsProperty<T>::operator[](idx);
        }
    };



public: //------------------------------------------------------ iterator types

    class Vertex_iterator
    {
    public:

        Vertex_iterator(Vertex v=Vertex(), const gsSurfMesh* m=NULL) : hnd_(v), mesh_(m)
        {
            if (mesh_ && mesh_->garbage()) while (mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
        }

        Vertex operator*()  const { return  hnd_; }

        const Vertex * operator->()  const { return  &hnd_; }

        // (!)
        size_t operator-(const Vertex_iterator& rhs) const
        {
            //assumes no deleted vertices..
            return (hnd_.idx_ - rhs.hnd_.idx_);
        }

        // (!)
        Vertex_iterator& operator+=(const size_t i)
        {
            //assumes no deleted vertices..
            hnd_.idx_+= i;
            return *this;
        }

        bool operator==(const Vertex_iterator& rhs) const
        {
            return (hnd_==rhs.hnd_);
        }

        bool operator<(const Vertex_iterator& rhs) const
        {
            return (hnd_<rhs.hnd_);
        }

        bool operator!=(const Vertex_iterator& rhs) const
        {
            return !operator==(rhs);
        }

        Vertex_iterator& operator++()
        {
            ++hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
            return *this;
        }

        Vertex_iterator& operator--()
        {
            --hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) --hnd_.idx_;
            return *this;
        }

    private:
        Vertex  hnd_;
        const gsSurfMesh* mesh_;
    };


    class Halfedge_iterator
    {
    public:

        Halfedge_iterator(Halfedge h=Halfedge(), const gsSurfMesh* m=NULL) : hnd_(h), mesh_(m)
        {
            if (mesh_ && mesh_->garbage()) while (mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
        }

        Halfedge operator*()  const { return  hnd_; }

        const Halfedge * operator->()  const { return  &hnd_; }

        bool operator==(const Halfedge_iterator& rhs) const
        {
            return (hnd_==rhs.hnd_);
        }

        bool operator!=(const Halfedge_iterator& rhs) const
        {
            return !operator==(rhs);
        }

        Halfedge_iterator& operator++()
        {
            ++hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
            return *this;
        }

        Halfedge_iterator& operator--()
        {
            --hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) --hnd_.idx_;
            return *this;
        }

    private:
        Halfedge  hnd_;
        const gsSurfMesh* mesh_;
    };


    class Edge_iterator
    {
    public:

        Edge_iterator(Edge e=Edge(), const gsSurfMesh* m=NULL) : hnd_(e), mesh_(m)
        {
            if (mesh_ && mesh_->garbage()) while (mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
        }

        Edge operator*()  const { return  hnd_; }

        const Edge * operator->()  const { return  &hnd_; }

        bool operator==(const Edge_iterator& rhs) const
        {
            return (hnd_==rhs.hnd_);
        }

        bool operator!=(const Edge_iterator& rhs) const
        {
            return !operator==(rhs);
        }

        Edge_iterator& operator++()
        {
            ++hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
            return *this;
        }

        Edge_iterator& operator--()
        {
            --hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) --hnd_.idx_;
            return *this;
        }

    private:
        Edge  hnd_;
        const gsSurfMesh* mesh_;
    };


    class Face_iterator
    {
    public:

        Face_iterator(Face f=Face(), const gsSurfMesh* m=NULL) : hnd_(f), mesh_(m)
        {
            if (mesh_ && mesh_->garbage()) while (mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
        }

        Face operator*()  const { return  hnd_; }

        const Face * operator->()  const { return  &hnd_; }

        // (!)
        size_t operator-(const Face_iterator& rhs) const
        {
            //assumes no deleted faces..
            return (hnd_.idx_ - rhs.hnd_.idx_);
        }

        // (!)
        Face_iterator& operator+=(const size_t i)
        {
            //assumes no deleted faces..
            hnd_.idx_+= i;
            return *this;
        }

        bool operator==(const Face_iterator& rhs) const
        {
            return (hnd_==rhs.hnd_);
        }

        bool operator<(const Face_iterator& rhs) const
        {
            return (hnd_<rhs.hnd_);
        }

        bool operator!=(const Face_iterator& rhs) const
        {
            return !operator==(rhs);
        }

        Face_iterator& operator++()
        {
            ++hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) ++hnd_.idx_;
            return *this;
        }

        Face_iterator& operator--()
        {
            --hnd_.idx_;
            assert(mesh_);
            while (mesh_->garbage() && mesh_->is_valid(hnd_) && mesh_->is_deleted(hnd_)) --hnd_.idx_;
            return *this;
        }

    private:
        Face  hnd_;
        const gsSurfMesh* mesh_;
    };



public: //-------------------------- containers for C++11 range-based for loops

    class Vertex_container
    {
    public:
        Vertex_container(Vertex_iterator _begin, Vertex_iterator _end) : begin_(_begin), end_(_end) {}
        Vertex_iterator begin() const { return begin_; }
        Vertex_iterator end()   const { return end_;   }
    private:
        Vertex_iterator begin_, end_;
    };



    class Halfedge_container
    {
    public:
        Halfedge_container(Halfedge_iterator _begin, Halfedge_iterator _end) : begin_(_begin), end_(_end) {}
        Halfedge_iterator begin() const { return begin_; }
        Halfedge_iterator end()   const { return end_;   }
    private:
        Halfedge_iterator begin_, end_;
    };



    class Edge_container
    {
    public:
        Edge_container(Edge_iterator _begin, Edge_iterator _end) : begin_(_begin), end_(_end) {}
        Edge_iterator begin() const { return begin_; }
        Edge_iterator end()   const { return end_;   }
    private:
        Edge_iterator begin_, end_;
    };



    class Face_container
    {
    public:
        Face_container(Face_iterator _begin, Face_iterator _end) : begin_(_begin), end_(_end) {}
        Face_iterator begin() const { return begin_; }
        Face_iterator end()   const { return end_;   }
    private:
        Face_iterator begin_, end_;
    };





public: //---------------------------------------------------- circulator types

    class Vertex_around_vertex_circulator
    {
    public:

        Vertex_around_vertex_circulator(const gsSurfMesh* m=NULL, Vertex v=Vertex())
        : mesh_(m), active_(true)
        {
            if (mesh_) halfedge_ = mesh_->halfedge(v);
        }

        bool operator==(const Vertex_around_vertex_circulator& rhs) const
        {
            assert(mesh_);
            return (active_ && (mesh_==rhs.mesh_) && (halfedge_==rhs.halfedge_));
        }

        bool operator!=(const Vertex_around_vertex_circulator& rhs) const
        {
            return !operator==(rhs);
        }

        Vertex_around_vertex_circulator& operator++()
        {
            assert(mesh_);
            halfedge_ = mesh_->ccw_rotated_halfedge(halfedge_);
            active_ = true;
            return *this;
        }

        Vertex_around_vertex_circulator& operator--()
        {
            assert(mesh_);
            halfedge_ = mesh_->cw_rotated_halfedge(halfedge_);
            return *this;
        }

        Vertex operator*()  const
        {
            assert(mesh_);
            return mesh_->to_vertex(halfedge_);
        }

        operator bool() const { return halfedge_.is_valid(); }

        Halfedge halfedge() const { return halfedge_; }

        // helper for C++11 range-based for-loops
        Vertex_around_vertex_circulator& begin() { active_=!halfedge_.is_valid(); return *this; }
        // helper for C++11 range-based for-loops
        Vertex_around_vertex_circulator& end()   { active_=true;  return *this; }

    private:
        const gsSurfMesh*  mesh_;
        Halfedge         halfedge_;
        // helper for C++11 range-based for-loops
        bool active_;
    };


    class Halfedge_around_vertex_circulator
    {
    public:

        Halfedge_around_vertex_circulator(const gsSurfMesh* m=NULL, Vertex v=Vertex())
        : mesh_(m), active_(true)
        {
            if (mesh_) halfedge_ = mesh_->halfedge(v);
        }

        bool operator==(const Halfedge_around_vertex_circulator& rhs) const
        {
            assert(mesh_);
            return (active_ && (mesh_==rhs.mesh_) && (halfedge_==rhs.halfedge_));
        }

        bool operator!=(const Halfedge_around_vertex_circulator& rhs) const
        {
            return !operator==(rhs);
        }

        Halfedge_around_vertex_circulator& operator++()
        {
            assert(mesh_);
            halfedge_ = mesh_->ccw_rotated_halfedge(halfedge_);
            active_ = true;
            return *this;
        }

        Halfedge_around_vertex_circulator& operator--()
        {
            assert(mesh_);
            halfedge_ = mesh_->cw_rotated_halfedge(halfedge_);
            return *this;
        }

        Halfedge operator*() const { return halfedge_; }

        const Halfedge * operator->()  const { return  &halfedge_; }

        operator bool() const { return halfedge_.is_valid(); }

        // helper for C++11 range-based for-loops
        Halfedge_around_vertex_circulator& begin() { active_=!halfedge_.is_valid(); return *this; }
        // helper for C++11 range-based for-loops
        Halfedge_around_vertex_circulator& end()   { active_=true;  return *this; }

    private:
        const gsSurfMesh*  mesh_;
        Halfedge         halfedge_;
        // helper for C++11 range-based for-loops
        bool active_;
    };


    class Face_around_vertex_circulator
    {
    public:

        Face_around_vertex_circulator(const gsSurfMesh* m=NULL, Vertex v=Vertex())
        : mesh_(m), active_(true)
        {
            if (mesh_)
            {
                halfedge_ = mesh_->halfedge(v);
                if (halfedge_.is_valid() && mesh_->is_boundary(halfedge_))
                    operator++();
            }
        }

        bool operator==(const Face_around_vertex_circulator& rhs) const
        {
            assert(mesh_);
            return (active_ && (mesh_==rhs.mesh_) && (halfedge_==rhs.halfedge_));
        }

        bool operator!=(const Face_around_vertex_circulator& rhs) const
        {
            return !operator==(rhs);
        }

        Face_around_vertex_circulator& operator++()
        {
            assert(mesh_ && halfedge_.is_valid());
            do {
                halfedge_ = mesh_->ccw_rotated_halfedge(halfedge_);
            } while (mesh_->is_boundary(halfedge_));
            active_ = true;
            return *this;
        }

        Face_around_vertex_circulator& operator--()
        {
            assert(mesh_ && halfedge_.is_valid());
            do
                halfedge_ = mesh_->cw_rotated_halfedge(halfedge_);
            while (mesh_->is_boundary(halfedge_));
            return *this;
        }

        Face operator*() const
        {
            assert(mesh_ && halfedge_.is_valid());
            return mesh_->face(halfedge_);
        }

        operator bool() const { return halfedge_.is_valid(); }

        // helper for C++11 range-based for-loops
        Face_around_vertex_circulator& begin() { active_=!halfedge_.is_valid(); return *this; }
        // helper for C++11 range-based for-loops
        Face_around_vertex_circulator& end()   { active_=true;  return *this; }

    private:
        const gsSurfMesh*  mesh_;
        Halfedge         halfedge_;
        // helper for C++11 range-based for-loops
        bool active_;
    };


    class Vertex_around_face_circulator
    {
    public:

        Vertex_around_face_circulator(const gsSurfMesh* m=NULL, Face f=Face())
        : mesh_(m), active_(true)
        {
            if (mesh_) halfedge_ = mesh_->halfedge(f);
        }

        bool operator==(const Vertex_around_face_circulator& rhs) const
        {
            assert(mesh_);
            return (active_ && (mesh_==rhs.mesh_) && (halfedge_==rhs.halfedge_));
        }

        bool operator!=(const Vertex_around_face_circulator& rhs) const
        {
            return !operator==(rhs);
        }

        Vertex_around_face_circulator& operator++()
        {
            assert(mesh_ && halfedge_.is_valid());
            halfedge_ = mesh_->next_halfedge(halfedge_);
            active_ = true;
            return *this;
        }

        Vertex_around_face_circulator& operator--()
        {
            assert(mesh_ && halfedge_.is_valid());
            halfedge_ = mesh_->prev_halfedge(halfedge_);
            return *this;
        }

        Vertex operator*() const
        {
            assert(mesh_ && halfedge_.is_valid());
            return mesh_->to_vertex(halfedge_);
        }

        Halfedge he() const
        {
            assert(mesh_ && halfedge_.is_valid());
            return halfedge_;
        }
        
        // helper for C++11 range-based for-loops
        Vertex_around_face_circulator& begin() { active_=false; return *this; }
        // helper for C++11 range-based for-loops
        Vertex_around_face_circulator& end()   { active_=true;  return *this; }

    private:
        const gsSurfMesh*  mesh_;
        Halfedge         halfedge_;
        // helper for C++11 range-based for-loops
        bool active_;
    };


    class Halfedge_around_face_circulator
    {
    public:

        Halfedge_around_face_circulator(const gsSurfMesh* m=NULL, Face f=Face())
        : mesh_(m), active_(true)
        {
            if (mesh_) halfedge_ = mesh_->halfedge(f);
        }

        bool operator==(const Halfedge_around_face_circulator& rhs) const
        {
            assert(mesh_);
            return (active_ && (mesh_==rhs.mesh_) && (halfedge_==rhs.halfedge_));
        }

        bool operator!=(const Halfedge_around_face_circulator& rhs) const
        {
            return !operator==(rhs);
        }

        Halfedge_around_face_circulator& operator++()
        {
            assert(mesh_ && halfedge_.is_valid());
            halfedge_ = mesh_->next_halfedge(halfedge_);
            active_ = true;
            return *this;
        }

        Halfedge_around_face_circulator& operator--()
        {
            assert(mesh_ && halfedge_.is_valid());
            halfedge_ = mesh_->prev_halfedge(halfedge_);
            return *this;
        }

        Halfedge operator*() const { return halfedge_; }

        const Halfedge * operator->()  const { return  &halfedge_; }

        // helper for C++11 range-based for-loops
        Halfedge_around_face_circulator& begin() { active_=false; return *this; }
        // helper for C++11 range-based for-loops
        Halfedge_around_face_circulator& end()   { active_=true;  return *this; }

    private:
        const gsSurfMesh*  mesh_;
        Halfedge         halfedge_;
        // helper for C++11 range-based for-loops
        bool active_;
    };



public: //-------------------------------------------- constructor / destructor



    gsSurfMesh();

    gsSurfMesh(const gsMatrix<Scalar> & pts);

    // destructor (is virtual, since we inherit from Geometry_representation)
    virtual ~gsSurfMesh();

    gsSurfMesh(const gsSurfMesh& rhs) { operator=(rhs); }

    gsSurfMesh& operator=(const gsSurfMesh& rhs);

    gsSurfMesh& assign(const gsSurfMesh& rhs);





public: //------------------------------------------------------------- file IO



    bool read(const std::string& filename);

    bool write(const std::string& filename) const;





public: //----------------------------------------------- add new vertex / face



    Vertex add_vertex(const Point& p);

    Face add_face(const std::vector<Vertex>& vertices);

    Face add_triangle(Vertex v1, Vertex v2, Vertex v3);

    Face add_quad(Vertex v1, Vertex v2, Vertex v3, Vertex v4);





public: //--------------------------------------------------- memory management



    unsigned int vertices_size() const { return (unsigned int) vprops_.size(); }
    unsigned int halfedges_size() const { return (unsigned int) hprops_.size(); }
    unsigned int edges_size() const { return (unsigned int) eprops_.size(); }
    unsigned int faces_size() const { return (unsigned int) fprops_.size(); }


    unsigned int n_vertices() const { return vertices_size() - deleted_vertices_; }
    unsigned int n_halfedges() const { return halfedges_size() - 2*deleted_edges_; }
    unsigned int n_edges() const { return edges_size() - deleted_edges_; }
    unsigned int n_faces() const { return faces_size() - deleted_faces_; }


    unsigned int empty() const { return n_vertices() == 0; }


    void clear();

    void free_memory();

    void reserve(unsigned int nvertices,
                 unsigned int nedges,
                 unsigned int nfaces );


    void garbage_collection();


    bool is_deleted(Vertex v) const
    {
        return vdeleted_[v];
    }
    bool is_deleted(Halfedge h) const
    {
        return edeleted_[edge(h)];
    }
    bool is_deleted(Edge e) const
    {
        return edeleted_[e];
    }
    bool is_deleted(Face f) const
    {
        return fdeleted_[f];
    }


    bool is_valid(Vertex v) const
    {
        return (0 <= v.idx()) && (v.idx() < (int)vertices_size());
    }
    bool is_valid(Halfedge h) const
    {
        return (0 <= h.idx()) && (h.idx() < (int)halfedges_size());
    }
    bool is_valid(Edge e) const
    {
        return (0 <= e.idx()) && (e.idx() < (int)edges_size());
    }
    bool is_valid(Face f) const
    {
        return (0 <= f.idx()) && (f.idx() < (int)faces_size());
    }





public: //---------------------------------------------- low-level connectivity



    Halfedge halfedge(Vertex v) const
    {
        return vconn_[v].halfedge_;
    }

    void set_halfedge(Vertex v, Halfedge h)
    {
        vconn_[v].halfedge_ = h;
    }

    bool is_boundary(Vertex v) const
    {
        Halfedge h(halfedge(v));
        return (!(h.is_valid() && face(h).is_valid()));
    }

    bool is_isolated(Vertex v) const
    {
        return !halfedge(v).is_valid();
    }

    bool is_manifold(Vertex v) const
    {
        // The vertex is non-manifold if more than one gap exists, i.e.
        // more than one outgoing boundary halfedge.
        int n(0);
        Halfedge_around_vertex_circulator hit = halfedges(v), hend=hit;
        if (hit) do
        {
            if (is_boundary(*hit))
                ++n;
        }
        while (++hit!=hend);
        return n<2;
    }

    Vertex to_vertex(Halfedge h) const
    {
        return hconn_[h].vertex_;
    }

    Vertex from_vertex(Halfedge h) const
    {
        return to_vertex(opposite_halfedge(h));
    }

    void set_vertex(Halfedge h, Vertex v)
    {
        hconn_[h].vertex_ = v;
    }

    Face face(Halfedge h) const
    {
        return hconn_[h].face_;
    }

    void set_face(Halfedge h, Face f)
    {
        hconn_[h].face_ = f;
    }

    Halfedge next_halfedge(Halfedge h) const
    {
        return hconn_[h].next_halfedge_;
    }

    void set_next_halfedge(Halfedge h, Halfedge nh)
    {
        hconn_[h].next_halfedge_ = nh;
        hconn_[nh].prev_halfedge_ = h;
    }

    Halfedge prev_halfedge(Halfedge h) const
    {
        return hconn_[h].prev_halfedge_;
    }

    Halfedge opposite_halfedge(Halfedge h) const
    {
        return Halfedge((h.idx() & 1) ? h.idx()-1 : h.idx()+1);
    }

    Halfedge ccw_rotated_halfedge(Halfedge h) const
    {
        return opposite_halfedge(prev_halfedge(h));
    }

    Halfedge cw_rotated_halfedge(Halfedge h) const
    {
        return next_halfedge(opposite_halfedge(h));
    }

    Edge edge(Halfedge h) const
    {
        return Edge(h.idx() >> 1);
    }

    bool is_boundary(Halfedge h) const
    {
        return !face(h).is_valid();
    }

    bool touches_boundary(Halfedge h) const
    {
        return is_boundary(opposite_halfedge(h));
    }

    Halfedge halfedge(Edge e, unsigned int i) const
    {
        assert(i<=1);
        return Halfedge((e.idx() << 1) + i);
    }

    Vertex vertex(Edge e, unsigned int i) const
    {
        assert(i<=1);
        return to_vertex(halfedge(e, i));
    }

    Face face(Edge e, unsigned int i) const
    {
        assert(i<=1);
        return face(halfedge(e, i));
    }

    bool is_boundary(Edge e) const
    {
        return (is_boundary(halfedge(e, 0)) || is_boundary(halfedge(e, 1)));
    }

    Halfedge halfedge(Face f) const
    {
        return fconn_[f].halfedge_;
    }

    void set_halfedge(Face f, Halfedge h)
    {
        fconn_[f].halfedge_ = h;
    }

    bool is_boundary(Face f) const
    {
        Halfedge h  = halfedge(f);
        Halfedge hh = h;
        do
        {
            if (is_boundary(opposite_halfedge(h)))
                return true;
            h = next_halfedge(h);
        }
        while (h != hh);
        return false;
    }





public: //--------------------------------------------------- property handling



    template <class T> Vertex_property<T> add_vertex_property(const std::string& name, T t=T())
    {
        return Vertex_property<T>(vprops_.add<T>(name, give(t)));
    }
    template <class T> Halfedge_property<T> add_halfedge_property(const std::string& name, T t=T())
    {
        return Halfedge_property<T>(hprops_.add<T>(name, give(t)));
    }
    template <class T> Edge_property<T> add_edge_property(const std::string& name, T t=T())
    {
        return Edge_property<T>(eprops_.add<T>(name, t));
    }
    template <class T> Face_property<T> add_face_property(const std::string& name, const T t=T())
    {
        return Face_property<T>(fprops_.add<T>(name, t));
    }
    template <class T> Mesh_property<T> add_mesh_property(const std::string& name, const T t=T())
    {
        return Mesh_property<T>(mprops_.add<T>(name, t));
    }

    template <class T> Vertex_property<T> get_vertex_property(const std::string& name) const
    {
        return Vertex_property<T>(vprops_.get<T>(name));
    }
    template <class T> Halfedge_property<T> get_halfedge_property(const std::string& name) const
    {
        return Halfedge_property<T>(hprops_.get<T>(name));
    }
    template <class T> Edge_property<T> get_edge_property(const std::string& name) const
    {
        return Edge_property<T>(eprops_.get<T>(name));
    }
    template <class T> Face_property<T> get_face_property(const std::string& name) const
    {
        return Face_property<T>(fprops_.get<T>(name));
    }
    template <class T> Mesh_property<T> get_mesh_property(const std::string& name) const
    {
        return Mesh_property<T>(mprops_.get<T>(name));
    }


    template <class T> Vertex_property<T> vertex_property(const std::string& name, const T t=T())
    {
        return Vertex_property<T>(vprops_.get_or_add<T>(name, give(t)));
    }
    template <class T> Halfedge_property<T> halfedge_property(const std::string& name, T t=T())
    {
        return Halfedge_property<T>(hprops_.get_or_add<T>(name, give(t)));
    }
    template <class T> Edge_property<T> edge_property(const std::string& name, const T t=T())
    {
        return Edge_property<T>(eprops_.get_or_add<T>(name, t));
    }
    template <class T> Face_property<T> face_property(const std::string& name, const T t=T())
    {
        return Face_property<T>(fprops_.get_or_add<T>(name, t));
    }

    template <class T> Mesh_property<T> mesh_property(const std::string& name, const T t=T())
    {
        return Mesh_property<T>(mprops_.get_or_add<T>(name, t));
    }

    template <class T> void rename_vertex_property(Vertex_property<T>& p,
                                                   std::string newname)
    {
        vprops_.rename(p, give(newname));
    }

    void swap_vertex_property(const std::string & name1,
                              const std::string & name2)
    {
        vprops_.swap(name1,name2);
    }

    template <class T> void remove_vertex_property(Vertex_property<T>& p)
    {
        vprops_.remove(p);
    }
    template <class T> void remove_halfedge_property(Halfedge_property<T>& p)
    {
        hprops_.remove(p);
    }
    template <class T> void remove_edge_property(Edge_property<T>& p)
    {
        eprops_.remove(p);
    }
    template <class T> void remove_face_property(Face_property<T>& p)
    {
        fprops_.remove(p);
    }
    template <class T> void remove_mesh_property(Mesh_property<T>& p)
    {
        mprops_.remove(p);
    }


    const std::type_info& get_vertex_property_type(const std::string& name)
    {
        return vprops_.get_type(name);
    }
    const std::type_info& get_halfedge_property_type(const std::string& name)
    {
        return hprops_.get_type(name);
    }
    const std::type_info& get_edge_property_type(const std::string& name)
    {
        return eprops_.get_type(name);
    }
    const std::type_info& get_face_property_type(const std::string& name)
    {
        return fprops_.get_type(name);
    }
    const std::type_info& get_mesh_property_type(const std::string& name)
    {
        return mprops_.get_type(name);
    }


    std::vector<std::string> vertex_properties() const
    {
        return vprops_.properties();
    }
    std::vector<std::string> halfedge_properties() const
    {
        return hprops_.properties();
    }
    std::vector<std::string> edge_properties() const
    {
        return eprops_.properties();
    }
    std::vector<std::string> face_properties() const
    {
        return fprops_.properties();
    }
    std::vector<std::string> mesh_properties() const
    {
        return mprops_.properties();
    }

    void property_stats() const;





public: //--------------------------------------------- iterators & circulators



    Vertex_iterator vertices_begin() const
    {
        return Vertex_iterator(Vertex(0), this);
    }

    Vertex_iterator vertices_end() const
    {
        return Vertex_iterator(Vertex(vertices_size()), this);
    }

    Vertex_container vertices() const
    {
        return Vertex_container(vertices_begin(), vertices_end());
    }

    Halfedge_iterator halfedges_begin() const
    {
        return Halfedge_iterator(Halfedge(0), this);
    }

    Halfedge_iterator halfedges_end() const
    {
        return Halfedge_iterator(Halfedge(halfedges_size()), this);
    }

    Halfedge_container halfedges() const
    {
        return Halfedge_container(halfedges_begin(), halfedges_end());
    }

    Edge_iterator edges_begin() const
    {
        return Edge_iterator(Edge(0), this);
    }

    Edge_iterator edges_end() const
    {
        return Edge_iterator(Edge(edges_size()), this);
    }

    Edge_container edges() const
    {
        return Edge_container(edges_begin(), edges_end());
    }

    Face_iterator faces_begin() const
    {
        return Face_iterator(Face(0), this);
    }

    Face_iterator faces_end() const
    {
        return Face_iterator(Face(faces_size()), this);
    }

    Face_container faces() const
    {
        return Face_container(faces_begin(), faces_end());
    }

    Vertex_around_vertex_circulator vertices(Vertex v) const
    {
        return Vertex_around_vertex_circulator(this, v);
    }

    Halfedge_around_vertex_circulator halfedges(Vertex v) const
    {
        return Halfedge_around_vertex_circulator(this, v);
    }

    Face_around_vertex_circulator faces(Vertex v) const
    {
        return Face_around_vertex_circulator(this, v);
    }

    Vertex_around_face_circulator vertices(Face f) const
    {
        return Vertex_around_face_circulator(this, f);
    }

    Halfedge_around_face_circulator halfedges(Face f) const
    {
        return Halfedge_around_face_circulator(this, f);
    }






public: //--------------------------------------------- higher-level operations



    bool is_triangle_mesh() const;

    bool is_quad_mesh() const;

    void triangulate();

    void triangulate(Face f);


    bool is_collapse_ok(Halfedge h);

    void collapse(Halfedge h);


    Vertex split(Face f, const Point& p) { Vertex v=add_vertex(p); split(f,v); return v; }

    void split(Face f, Vertex v);


    void quad_split(Face f, Vertex v, Halfedge s);
    void quad_split();
    
    Halfedge split(Edge e, const Point& p) { return split(e, add_vertex(p)); }


    Halfedge split(Edge e, Vertex v);


    Halfedge insert_vertex(Edge e, const Point& p)
    {
        return insert_vertex(halfedge(e,0), add_vertex(p));
    }

    Halfedge insert_vertex(Edge e, Vertex v)
    {
        return insert_vertex(halfedge(e,0), v);
    }

    Halfedge insert_vertex(Halfedge h, Vertex v);


    Halfedge insert_edge(Halfedge h0, Halfedge h1);


    bool is_flip_ok(Edge e) const;

    void flip(Edge e);


    unsigned int valence(Vertex v) const;

    unsigned int valence(Face f) const;

    unsigned int face_valence_sum() const;

    Halfedge find_halfedge(Vertex start, Vertex end) const;

    Edge find_edge(Vertex a, Vertex b) const;

    void delete_vertex(Vertex v);

    void delete_edge(Edge e);

    void delete_face(Face f);



public: //------------------------------------------ geometry-related functions



    const Point& position(Vertex v) const { return vpoint_[v]; }

    Point& position(Vertex v) { return vpoint_[v]; }

    std::vector<Point>& points() { return vpoint_.vector(); }

    void update_face_normals();

    Normal compute_face_normal(Face f) const;

    void update_vertex_normals();

    Normal compute_vertex_normal(Vertex v) const;

    Scalar edge_length(Edge e) const;





private: //---------------------------------------------- allocate new elements

    Vertex new_vertex()
    {
        vprops_.push_back();
        return Vertex(vertices_size()-1);
    }

    Halfedge new_edge(Vertex start, Vertex end)
    {
        assert(start != end);

        eprops_.push_back();
        hprops_.push_back();
        hprops_.push_back();

        Halfedge h0(halfedges_size()-2);
        Halfedge h1(halfedges_size()-1);

        set_vertex(h0, end); //to_vertex*
        set_vertex(h1, start); //to_vertex

        return h0;
    }

    Face new_face()
    {
        fprops_.push_back();
        return Face(faces_size()-1);
    }

    friend std::ostream& operator<<(std::ostream& os, const gsSurfMesh & sm)
    {
        os<<"gsSurfMesh with "<<sm.n_vertices()<<" vertices, "<<sm.n_edges()<<
            " edges and "<<sm.n_faces()<<" faces.\n";
        return os;
    }

public: // Catmull-Clark functions

    void cc_subdivide();

    Vertex_property<Point> cc_limit_points(std::string label = "v:limit");

    Vertex_property<Point> cc_limit_normals(std::string label = "v:normal",
                                            bool normalize = true);

    Vertex_property<Point> cc_limit_tangent_vec(std::string label = "v:tanvec",
                                                bool normalize = true);

    gsMultiPatch<real_t> cc_acc3(bool comp_topology = false) const;

    gsMultiPatch<real_t> linear_patches() const;
    
private: //--------------------------------------------------- helper functions

    void adjust_outgoing_halfedge(Vertex v);

    void remove_edge(Halfedge h);

    void remove_loop(Halfedge h);

    bool garbage() const { return garbage_; }

private: //------------------------------------------------------- private data

    friend bool GISMO_EXPORT read_poly(gsSurfMesh& mesh, const std::string& filename);

    gsProperty_container vprops_;
    gsProperty_container hprops_;
    gsProperty_container eprops_;
    gsProperty_container fprops_;
    gsProperty_container mprops_;

    Vertex_property<Vertex_connectivity>      vconn_;
    Halfedge_property<Halfedge_connectivity>  hconn_;
    Face_property<Face_connectivity>          fconn_;

    Vertex_property<bool>  vdeleted_;
    Edge_property<bool>    edeleted_;
    Face_property<bool>    fdeleted_;

    Vertex_property<Point>   vpoint_;
    Vertex_property<Normal>  vnormal_;
    Face_property<Normal>    fnormal_;

    unsigned int deleted_vertices_;
    unsigned int deleted_edges_;
    unsigned int deleted_faces_;
    bool garbage_;

    // helper data for add_face()
    typedef std::pair<Halfedge, Halfedge>  NextCacheEntry;
    typedef std::vector<NextCacheEntry>    NextCache;
    std::vector<Vertex>      add_face_vertices_;
    std::vector<Halfedge>    add_face_halfedges_;
    std::vector<bool>        add_face_is_new_;
    std::vector<bool>        add_face_needs_adjust_;
    NextCache                add_face_next_cache_;
};


//------------------------------------------------------------ output operators


inline std::ostream& operator<<(std::ostream& os, gsSurfMesh::Vertex v)
{
    return (os << 'v' << v.idx());
}

inline std::ostream& operator<<(std::ostream& os, gsSurfMesh::Halfedge h)
{
    return (os << 'h' << h.idx());
}

inline std::ostream& operator<<(std::ostream& os, gsSurfMesh::Edge e)
{
    return (os << 'e' << e.idx());
}

inline std::ostream& operator<<(std::ostream& os, gsSurfMesh::Face f)
{
    return (os << 'f' << f.idx());
}



namespace internal
{

template<>
class GISMO_EXPORT gsXml<gsSurfMesh>
{
private:
    gsXml();
public:
    GSXML_COMMON_FUNCTIONS(gsSurfMesh)
    GSXML_GET_POINTER(gsSurfMesh)
    static std::string tag () { return "Mesh"; }
    static std::string type() { return "off"; }

    static void get_into(gsXmlNode * node, gsSurfMesh & result);
    static gsXmlNode * put (const gsSurfMesh & obj, gsXmlTree & data);
};

}//namespace internal

//=============================================================================
} // namespace gismo
//=============================================================================