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
//=============================================================================