inmost_mesh.h

#ifndef INMOST_MESH_H_INCLUDED
#define INMOST_MESH_H_INCLUDED

#include "inmost_common.h"
#include "inmost_data.h"

#if defined(USE_MESH)



namespace INMOST
{


	class Mesh;
	class Storage;
	class Element;
	class TagManager;
	class Node;
	class Edge;
	class Face;
	class Cell;
	class ElementSet;


	
	typedef INMOST_DATA_BULK_TYPE GeometricData;
	static const GeometricData CENTROID    = 0;
	static const GeometricData NORMAL      = 1;
	static const GeometricData ORIENTATION = 2;
	static const GeometricData MEASURE     = 3;
	static const GeometricData BARYCENTER  = 4;
	
	typedef INMOST_DATA_BULK_TYPE SyncBitOp; //< This type is used for marker synchronization
	static const SyncBitOp SYNC_BIT_SET = 0;
	static const SyncBitOp SYNC_BIT_OR  = 1;
	static const SyncBitOp SYNC_BIT_XOR = 2;
	static const SyncBitOp SYNC_BIT_AND = 3;


	///Definition for data type of topology error or event. This type may be extended later to 64 bits
  /// to accomodate more topological errors or events.
	typedef INMOST_DATA_ENUM_TYPE         TopologyCheck; 
  ///Throw TopologyError exception on error.
	static const TopologyCheck            THROW_EXCEPTION        = 0x00000001; 
  ///Print topology notify to std::cerr.
	static const TopologyCheck            PRINT_NOTIFY           = 0x00000002; 
  ///Element should be deleted if there is an error in EndTopologyCheck.
	static const TopologyCheck            DELETE_ON_ERROR        = 0x00000004; 
  ///Bitwise OR of errors is recorded to sparse integer tag of element (if not deleted), availible 
  /// by TopologyErrorTag().
	static const TopologyCheck            MARK_ON_ERROR          = 0x00000008; 
  ///Check that edge already exists, on occurance CreateEdge silently returns existant edge.
	static const TopologyCheck            DUPLICATE_EDGE         = 0x00000010; 
  ///Check that face already exists, on occurance CreateFace silently returns existant face.
	static const TopologyCheck            DUPLICATE_FACE         = 0x00000020; 
  ///Check that cell already exists, on occurance CreateCell silently returns existant cell.
	static const TopologyCheck            DUPLICATE_CELL         = 0x00000040; 
  ///Check that edge have more then two nodes (no support for this kind of object).
	static const TopologyCheck            DEGENERATE_EDGE        = 0x00000080; 
  ///Produce error if face consists of less then 3 edges.
	static const TopologyCheck            DEGENERATE_FACE        = 0x00000100; 
  ///Produce error if cell consists of less then 4 faces.
	static const TopologyCheck            DEGENERATE_CELL        = 0x00000200; 
  ///Produce error if face have wrong orientation of edges, for star-shaped elements only.
	static const TopologyCheck            FACE_ORIENTATION       = 0x00000400; 
  ///Produce error if face is non planar.
	static const TopologyCheck            FACE_PLANARITY         = 0x00000800; 
  ///Produce error if there is another face that already uses same nodes.
	static const TopologyCheck            INTERLEAVED_FACES      = 0x00001000; 
  ///Check that every face have exectly two neighbours, if DUPLICATE_CELL is activated, then checks 
  /// for cell duplication.
	static const TopologyCheck            TRIPLE_SHARED_FACE     = 0x00002000;
  ///Produce error if one of the faces of the cell contains all the nodes of the cell.
	static const TopologyCheck            FLATTENED_CELL         = 0x00004000;
  ///Produce error if provided array of elements for construction contain duplications.
	static const TopologyCheck            ADJACENT_DUPLICATE     = 0x00008000;
  ///Hidden elements should not be used when new elements are created.
	static const TopologyCheck            ADJACENT_HIDDEN        = 0x00010000; 
  ///Check that all handles are valid.
	static const TopologyCheck            ADJACENT_VALID         = 0x00020000; 
  ///Produce error if provided array of elements have wrong geometric dimension.
	static const TopologyCheck            ADJACENT_DIMENSION     = 0x00040000; 
  ///Produce error if edges of faces are not closed, or face is folded.
  /// \warning This check needs NEED_TEST_CLOSURE to be activated.
	static const TopologyCheck            PROHIBIT_MULTILINE     = 0x00080000;
  /// Allow only known types of elements: Tet, Quad, Line.
	static const TopologyCheck            PROHIBIT_POLYGON       = 0x00100000;
  ///Produce error if faces of cell are not closed, or cell is folded.
  /// \warning This check needs NEED_TEST_CLOSURE to be activated.
	static const TopologyCheck            PROHIBIT_MULTIPOLYGON  = 0x00200000; 
  ///Allow only known types of elements: Tet,Hex,Prism,Pyramid.
	static const TopologyCheck            PROHIBIT_POLYHEDRON    = 0x00400000; 
  ///Edges of the face should form one closed loop.
	static const TopologyCheck            FACE_EDGES_ORDER       = 0x00800000; 
  ///Don't allow concave face.
  /// \warning Not implemented.
	static const TopologyCheck            PROHIBIT_CONCAVE_FACE  = 0x01000000; 
  ///Don't allow concave cell.
  /// \warning Not implemented.
	static const TopologyCheck            PROHIBIT_CONCAVE_CELL  = 0x02000000;
  ///Don't allow non-star shaped face.
  /// \warning Not implemented.
	static const TopologyCheck            PROHIBIT_NONSTAR_FACE  = 0x04000000;
  ///Don't allow non-star shaped concave cell.
  /// \warning Not implemented.
	static const TopologyCheck            PROHIBIT_NONSTAR_CELL  = 0x08000000;
  ///Edges of the face don't cross each other.
  /// \warning Not implemented.
	static const TopologyCheck            FACE_SELF_INTERSECTION = 0x10000000;
  ///Faces of the cell don't cross each other.
  /// \warning Not implemented.
	static const TopologyCheck            CELL_SELF_INTERSECTION = 0x20000000; 
  ///Silent, test's for closure in ComputeGeometricType, needed to detect MultiLine and 
  /// MultiPolygon.
	static const TopologyCheck            NEED_TEST_CLOSURE      = 0x40000000; 
  ///Don't allow 2d grids, where edges appear to be vertexes, faces are edges and cells are faces
	static const TopologyCheck            DISABLE_2D             = 0x80000000;
  ///A shortcut to test for grid conformity.
	static const TopologyCheck            GRID_CONFORMITY        = NEED_TEST_CLOSURE 
                                                               | PROHIBIT_MULTILINE 
                                                               | PROHIBIT_MULTIPOLYGON  
                                                               | INTERLEAVED_FACES 
                                                               | TRIPLE_SHARED_FACE;
  ///Default set of options.
	static const TopologyCheck            DEFAULT_CHECK          = THROW_EXCEPTION 
                                                               | DUPLICATE_EDGE 
                                                               | DUPLICATE_FACE 
                                                               | DUPLICATE_CELL 
                                                               | PRINT_NOTIFY;
  ///Returns a string explaining each topology check. Implemented in mesh.cpp.
  /// @param c Single topology check.
	const char * TopologyCheckNotifyString(TopologyCheck c);
  
 
	
	template <typename StorageType>
	class ElementArray
	{
	public:
		typedef dynarray<HandleType,64>      cont_t;
		typedef cont_t::size_type            size_type;
	private:
		Mesh *                               m_link;
		cont_t                               container;
	public:
		const cont_t & GetContainer() {return container;}
		~ElementArray() {}
		ElementArray() : m_link(NULL) {}
		ElementArray(Mesh * m_link) : m_link(m_link) {}
		ElementArray(Mesh * m_link, size_type n, HandleType h = InvalidHandle())  : m_link(m_link), container(n,h) {}
		ElementArray(Mesh * m_link, const cont_t & c)  : m_link(m_link), container(c) {}
		ElementArray(const ElementArray & other) : m_link(other.m_link), container(other.container) {}
		template<class InputIterator>
		ElementArray(Mesh * m_link, InputIterator first, InputIterator last) :m_link(m_link), container(first,last) {}
		ElementArray & operator=(ElementArray const & other) {m_link = other.m_link; container = other.container; return *this;}
	public:
		class iterator : public cont_t::iterator
		{
			Mesh * m_link;
		public:
			iterator(Mesh * m, const cont_t::iterator & other ) : cont_t::iterator(other) , m_link(m){assert(m_link);}
			iterator(const iterator & other ) : cont_t::iterator(other), m_link(other.m_link) {assert(m_link);}
      ptrdiff_t     operator -(const iterator & other) const {return static_cast<const cont_t::iterator>(*this)-static_cast<const cont_t::iterator>(other);}
      iterator      operator +(size_t n){iterator ret(*this); ret.cont_t::iterator::operator+(n); return ret;}
      iterator      operator -(size_t n){iterator ret(*this); ret.cont_t::iterator::operator-(n); return ret;}
			iterator &    operator ++() {cont_t::iterator::operator++(); return *this;}
			iterator      operator ++(int) {iterator ret(*this); cont_t::iterator::operator++(); return ret;}
			iterator &    operator --() {cont_t::iterator::operator--(); return *this;}
			iterator      operator --(int) {iterator ret(*this); cont_t::iterator::operator--(); return ret;}
			iterator &    operator =(iterator const & other) {m_link = other.m_link; cont_t::iterator::operator=(static_cast<cont_t::iterator const &>(other)); return *this; }
			HandleType &  operator *() { return cont_t::iterator::operator *(); }
			StorageType   operator->() { return StorageType(m_link,&cont_t::iterator::operator *()); }
		};
		class reverse_iterator : public cont_t::reverse_iterator
		{
			Mesh * m_link;
		public:
			reverse_iterator(Mesh * m, const cont_t::reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(m) {assert(m_link);}
			reverse_iterator(const reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(other.m_link) {assert(m_link);}
      ptrdiff_t             operator -(const reverse_iterator & other) const {return static_cast<const cont_t::reverse_iterator>(other)-static_cast<const cont_t::reverse_iterator>(*this);}
      reverse_iterator      operator +(size_t n){reverse_iterator ret(*this); ret.cont_t::reverse_iterator::operator+(n); return ret;}
      reverse_iterator      operator -(size_t n){reverse_iterator ret(*this); ret.cont_t::reverse_iterator::operator-(n); return ret;}
			reverse_iterator &    operator ++() {cont_t::reverse_iterator::operator++(); return *this;}
			reverse_iterator      operator ++(int) {reverse_iterator ret(*this); cont_t::reverse_iterator::operator++(); return ret;}
			reverse_iterator &    operator --() {cont_t::reverse_iterator::operator--(); return *this;}
			reverse_iterator      operator --(int) {reverse_iterator ret(*this); cont_t::reverse_iterator::operator--(); return ret;}
			reverse_iterator & operator =(reverse_iterator const & other) {m_link = other.m_link; cont_t::reverse_iterator::operator=(static_cast<cont_t::reverse_iterator const &>(other)); return *this; }
			HandleType &       operator *() { return cont_t::reverse_iterator::operator *(); }
			StorageType        operator->() { return StorageType(m_link,&cont_t::reverse_iterator::operator *()); }
		};
		class const_iterator : public cont_t::const_iterator
		{
			Mesh * m_link;
		public:
			const_iterator(Mesh * m, const cont_t::const_iterator & other ) : cont_t::const_iterator(other), m_link(m) {assert(m_link);}
			const_iterator(const const_iterator & other ) : cont_t::const_iterator(other), m_link(other.m_link) {assert(m_link);}
			const_iterator &    operator ++() {cont_t::const_iterator::operator++(); return *this;}
			const_iterator      operator ++(int) {const_iterator ret(*this); cont_t::const_iterator::operator++(); return ret;}
			const_iterator &    operator --() {cont_t::const_iterator::operator--(); return *this;}
			const_iterator      operator --(int) {const_iterator ret(*this); cont_t::const_iterator::operator--(); return ret;}
			const_iterator &    operator =(const_iterator const & other) {m_link = other.m_link; cont_t::const_iterator::operator=(static_cast<cont_t::const_iterator const &>(other)); return *this; }
			const HandleType &  operator *() { return cont_t::const_iterator::operator *(); }
			StorageType         operator->() { return StorageType(m_link,cont_t::const_iterator::operator *()); }
		};
		class const_reverse_iterator : public cont_t::const_reverse_iterator
		{
			Mesh * m_link;
		public:
			const_reverse_iterator(Mesh * m, const cont_t::const_reverse_iterator & other) : cont_t::const_reverse_iterator(other), m_link(m) {assert(m_link);}
			const_reverse_iterator(const const_reverse_iterator & other ) : cont_t::const_reverse_iterator(other), m_link(other.m_link) {assert(m_link);}
			const_reverse_iterator &    operator ++() {cont_t::const_reverse_iterator::operator++(); return *this;}
			const_reverse_iterator      operator ++(int) {const_reverse_iterator ret(*this); cont_t::const_reverse_iterator::operator++(); return ret;}
			const_reverse_iterator &    operator --() {cont_t::const_reverse_iterator::operator--(); return *this;}
			const_reverse_iterator      operator --(int) {const_reverse_iterator ret(*this); cont_t::const_reverse_iterator::operator--(); return ret;}
			const_reverse_iterator & operator =(const_reverse_iterator const & other) { cont_t::const_reverse_iterator::operator=(static_cast<cont_t::const_reverse_iterator const &>(other)); return *this; }
			const HandleType &       operator *() { return cont_t::const_reverse_iterator::operator *(); }
			StorageType              operator->() { return StorageType(m_link,cont_t::const_reverse_iterator::operator *()); }
		};
	public:
		template<class InputIterator>
		__INLINE void             insert      (iterator pos,InputIterator pbeg, InputIterator pend) {container.insert(static_cast<cont_t::iterator>(pos),pbeg,pend);}
		__INLINE iterator         erase       (iterator pos) {return iterator(m_link,container.erase(static_cast<cont_t::iterator>(pos)));}

		__INLINE iterator         begin       () { return iterator(m_link,container.begin()); }
		__INLINE iterator         end         () { return iterator(m_link,container.end()); }
		__INLINE reverse_iterator rbegin      () { return reverse_iterator(m_link,container.rbegin()); }
		__INLINE reverse_iterator rend        () { return reverse_iterator(m_link,container.rend()); }
		__INLINE const_iterator         begin () const { return const_iterator(m_link,container.begin()); }
		__INLINE const_iterator         end   () const { return const_iterator(m_link,container.end()); }
		__INLINE const_reverse_iterator rbegin() const { return const_reverse_iterator(m_link,container.rbegin()); }
		__INLINE const_reverse_iterator rend  () const { return const_reverse_iterator(m_link,container.rend()); }
		
		__INLINE StorageType      operator [] (size_type n) {assert(m_link); return StorageType(m_link,&container[n]);}
		__INLINE StorageType      operator [] (size_type n) const {assert(m_link); return StorageType(m_link,container[n]);}
		__INLINE StorageType      front       () {assert(m_link); return StorageType(m_link,&container.front()); }
		__INLINE StorageType      front       () const {assert(m_link); return StorageType(m_link,container.front()); }
		__INLINE StorageType      back        ()  {assert(m_link); return StorageType(m_link,&container.back()); }
		__INLINE StorageType      back        () const {assert(m_link); return StorageType(m_link,container.back()); }
		__INLINE HandleType       atfront     () const { return container.front(); }
		__INLINE HandleType       atback      () const { return container.back(); }
		__INLINE HandleType &     atfront     () { return container.front(); }
		__INLINE HandleType &     atback      () { return container.back(); }
		__INLINE HandleType &     at          (size_type n) { return container.at(n); }
		__INLINE HandleType       at          (size_type n) const { return container.at(n); }
		__INLINE void             swap        (ElementArray<StorageType> & other) {Mesh * t = m_link; m_link = other.m_link; other.m_link = t; container.swap(other.container);}
		__INLINE void             push_back   (const Storage & x) {if( m_link == NULL ) m_link = x.GetMeshLink(); assert(x.GetMeshLink() == m_link); container.push_back(x.GetHandle());}
		//__INLINE void             push_back   (const StorageType & x) {container.push_back(x.GetHandle());}
		__INLINE void             push_back   (HandleType x) {assert(m_link != NULL); container.push_back(x);}
		__INLINE void             pop_back    () {container.pop_back();}
		__INLINE void             resize      (size_type n) {container.resize(n);}
		__INLINE bool             empty       () const {return container.empty();}
		__INLINE void             clear       () {container.clear();}
		__INLINE void             reserve     (size_type n) {container.reserve(n);}
		__INLINE size_type        size        () const {return container.size(); }
		__INLINE HandleType *     data        () {return container.data();}
		__INLINE const HandleType*data        () const {return container.data();}
		__INLINE Mesh *           GetMeshLink () const {assert(m_link); return m_link;}
		__INLINE void             SetMeshLink (Mesh * m) {m_link = m;}
		//implemented in mesh_array.cpp
		void                      Unite       (const HandleType * h, INMOST_DATA_ENUM_TYPE num);
		void                      Subtract    (const HandleType * h, INMOST_DATA_ENUM_TYPE num);
		void                      Intersect   (const HandleType * h, INMOST_DATA_ENUM_TYPE num);
		template<typename EType>
		void                      Unite       (const ElementArray<EType> & other) {Unite(other.data(),static_cast<INMOST_DATA_ENUM_TYPE>(other.size()));}
		template<typename EType>
		void                      Subtract    (const ElementArray<EType> & other) {Subtract(other.data(),static_cast<INMOST_DATA_ENUM_TYPE>(other.size()));}
		template<typename EType>
		void                      Intersect   (const ElementArray<EType> & other) {Intersect(other.data(),static_cast<INMOST_DATA_ENUM_TYPE>(other.size()));}
		void                      SetMarker   (MarkerType n) const;
		void                      RemMarker   (MarkerType n) const;
    void                      SetPrivateMarker   (MarkerType n) const;
		void                      RemPrivateMarker   (MarkerType n) const;
    template<typename Etype>
    ElementArray<Etype>       Convert() {return ElementArray<Etype>(m_link,container);}
	};
	
			
	class Element : public Storage //implemented in element.cpp
	{
	public:
		typedef INMOST_DATA_BULK_TYPE GeometricType;
		static const GeometricType Unset        = 0;
		static const GeometricType Vertex       = 1;
		static const GeometricType Line         = 2;
		static const GeometricType MultiLine    = 3;
		static const GeometricType Tri          = 4;
		static const GeometricType Quad         = 5;
		static const GeometricType Polygon      = 6;
		static const GeometricType MultiPolygon = 7;
		static const GeometricType Tet          = 8;
		static const GeometricType Hex          = 9;
		static const GeometricType Prism        = 10;
		static const GeometricType Pyramid      = 11;
		static const GeometricType Polyhedron   = 12;
		static const GeometricType Set          = 253;
		static const GeometricType MeshPart     = 254;
		static const GeometricType MaxType      = 255;
		//enum GeometricType {Unset,Vertex,Line,MultiLine,Tri,Quad,Polygon,MultiPolygon,Tet,Hex,Prism,Pyramid,Polyhedron,Set};
		static const char *       GeometricTypeName(GeometricType t);
		static integer            GetGeometricDimension(GeometricType m_type);
		typedef INMOST_DATA_BULK_TYPE Status;
		static const Status Owned  = 1;
		static const Status Shared = 2;
		static const Status Ghost  = 4;
		static const Status Any    = 0;
		static const char * StatusName(Status s);
	public:
		typedef inner_reference_array            adj_type;
		typedef adj_type::iterator               adj_iterator;
		typedef adj_type::const_iterator         const_adj_iterator;
		typedef adj_type::reverse_iterator       adj_reverse_iterator;
		typedef adj_type::const_reverse_iterator const_adj_reverse_iterator;
	public:
		//adj_type &                  HighConn                () const;
		//adj_type &                  LowConn                 () const;
	protected:
		void                        SetGeometricType        (GeometricType t);
	public:
		Element() : Storage(NULL,InvalidHandle()) {}
		Element(Mesh * m, HandleType h) : Storage(m,h) {}
		Element(Mesh * m, HandleType * h) : Storage(m,h) {}
		Element(const Element & other) : Storage(other) {}
		Element & operator =(Element const & other) {Storage::operator =(other); return *this;}
		Element * operator ->() {return this;}
		const Element * operator ->() const {return this;}
		Element & self() {return *this;}
		const Element & self() const {return *this;}
		virtual ~Element() {}
	public:
		/// Retrieve number of adjacent elements.
		/// For etype you can either pass one type as CELL,
		/// or several types as bitwise mask: NODE | CELL.
		/// @param etype bitwise mask of element types
		/// @return number of adjacent elements.
		virtual enumerator                nbAdjElements           (ElementType etype) const;
		/// Retrieve unordered array of adjacent elements.
		/// If you care about orderness of nodes for face you should you Face::getNodes() instead.
		/// If you want faster access you may use direct access to handles stored in memory
		/// through Mesh::HighConn for upper adjacencies and Mesh::LowConn for lower adjacencies.
		/// @param etype bitwise mask of element types
		/// @return array of elements
		virtual ElementArray<Element>     getAdjElements          (ElementType etype) const;  //unordered
		/// Retrieve number of adjacent elements with marker.
		/// As etype you can either pass one type as CELL,
		/// or several types as bitwise mask: NODE | CELL
		/// @param etype bitwise mask of element types
		/// @param mask marker to be set 
		/// @param invert_mask if true then those are selected on wich marker is not set
		/// @return number of adjacent elements.
		virtual enumerator                nbAdjElements           (ElementType etype, MarkerType mask, bool invert_mask = false) const;
		/// Retrieve unordered array of adjacent elements with marker.
		/// @param etype bitwise mask of element types
		/// @param mask marker to be set 
		/// @param invert_mask if true then those are selected on wich marker is not set
		/// @return array of elements
		virtual ElementArray<Element>     getAdjElements          (ElementType etype, MarkerType mask, bool invert_mask = false) const;  //unordered
		ElementArray<Element>       BridgeAdjacencies       (ElementType Bridge, ElementType Dest, MarkerType mask = 0, bool invert_mask = false) const;
		ElementArray<Node>          BridgeAdjacencies2Node  (ElementType Bridge, MarkerType mask = 0, bool invert_mask = false) const;
		ElementArray<Edge>          BridgeAdjacencies2Edge  (ElementType Bridge, MarkerType mask = 0, bool invert_mask = false) const;
		ElementArray<Face>          BridgeAdjacencies2Face  (ElementType Bridge, MarkerType mask = 0, bool invert_mask = false) const;
		ElementArray<Cell>          BridgeAdjacencies2Cell  (ElementType Bridge, MarkerType mask = 0, bool invert_mask = false) const;
		/// Retrieve all the nodes of the element.
		///
		/// For a node returns itself.
		///
		/// For an edge returns ordered pair of nodes. The order of nodes in the edge is preserved from the first construction.
		///
		/// For a face returns ordered set of nodes. In the case Face::CheckNormalOrientation returns true the
		/// order of nodes will follow right hand side rule with respect to normal vector, otherwise it follows
		/// left hand side rule with respect to normal vector.
		///
		/// For a cell returns the same order that was provided through suggest_nodes_oreder in Mesh::CreateCell.
		/// In the case suggest_nodes_order was not provided, the order of nodes follows VTK format for known types
		/// of elements such as Element::Tet, Element::Prism, Element::Hex, Element::Pyramid. For a general polyhedron
		/// the order is unspecified.
		///
		/// @return array of elements
		/// @see Face::CheckNormalOrientation
		/// @see Face::FaceOrientedOutside
		virtual ElementArray<Node>  getNodes                () const; //unordered
		/// Retrieve all the edges of the element.
		///
		/// For a node returns unordered set of edges.
		///
		/// For an edge returns itself.
		///
		/// For a face returns ordered set of edges.
		///
		/// For a cell returns unordered set of edges.
		///
		/// @return array of elements
		virtual ElementArray<Edge>  getEdges                () const; //unordered
		/// Retrieve all the faces of the element.
		///
		/// For a node returns unordered set of faces.
		///
		/// For an edge returns unordered set of faces.
		///
		/// For a face returns itself.
		///
		/// For a cell return ordered set of faces. The order of faces in the cell is preserved from the first construction.
		///
		/// @return array of elements
		virtual ElementArray<Face>  getFaces                () const; //unordered
		/// Return all the cells of the element.
		///
		/// For a node returns unordered set of cells.
		///
		/// For an edge returns unordered set of cells.
		///
		/// For a face returns a pair of cells. In the case Face::CheckNormalOrientation returns true
		/// then the normal points from the first cell to the second and in oppisite direction otherwise.
		///
		/// For a cell returns itself.
		///
		/// @return array of elements
		/// @see Face::FaceOrientedOutside
		virtual ElementArray<Cell>  getCells                () const; //unordered
		virtual ElementArray<Node>  getNodes                (MarkerType mask,bool invert_mask = false) const; //unordered
		virtual ElementArray<Edge>  getEdges                (MarkerType mask,bool invert_mask = false) const; //unordered
		virtual ElementArray<Face>  getFaces                (MarkerType mask,bool invert_mask = false) const; //unordered
		virtual ElementArray<Cell>  getCells                (MarkerType mask,bool invert_mask = false) const; //unordered
		GeometricType               GetGeometricType        () const;
		unsigned int                GetElementDimension     () const {return GetGeometricDimension(GetGeometricType());}
		Status                      GetStatus               () const;
		void                        SetStatus               (Status status) const;
		Storage::integer &          GlobalID                () const;
		bool                        CheckElementConnectivity() const;
		void                        PrintElementConnectivity() const;
		static bool                 CheckConnectivity       (Mesh * m);
		//implemented in geometry.cpp
		void                        CastRay                 (real * pos, real * dir, dynarray< std::pair<Element, real> , 16 > & hits) const;
		void                        ComputeGeometricType    () const;
		void                        Centroid                (real * cnt) const;
		void                        Barycenter              (real * cnt) const;
		Storage::real               Mean                    (real (*func)(real* x,real t),real time) const;
		bool                        Boundary                () const;
		bool                        Planarity               () const; // check that all nodes lay on one plane
		//implemented in modify.cpp
		bool                        Hide                    () const; // if true then element was hidden, works only inside BeginModification and EndModification, on EndModification all Hidden elements are deleted
		bool                        Show                    () const; // if true then element was recovered
    /// \warning This function will not resolve an ierarchical strucutre of ElementSet, use ElemnetSet::DeleteSet instead
		bool                        Delete                  (); // if true then element was deleted, otherwise it was hidden
		bool                        Hidden                  () const;
		bool                        New                     () const;
		void                        Disconnect              (bool delete_upper_adjacent) const; //disconnect all elements, delete upper dependent
		/// Disconnect element.
		/// Disconnects nodes from this edge, edges from this face, faces from this cell, cannot disconnect cells from this node; \n
		/// Disconnects edges from this node, faces from this edge, cells from this face, cannot disconnect nodes from this cell; \n
		/// Updates geometric data and cell nodes automatically.
		void                        Disconnect              (const HandleType * adjacent, INMOST_DATA_ENUM_TYPE num) const;
		/// \brief Connects lower adjacencies to current element.
		/// Geometric data and cell nodes are updated automatically.
		///
		/// \todo
		///	  1. Asserts in this function should be replaced by Topography checks.
		///	  2. This function should be used for construction of elements instead of current implementation.
		///	  3. Should correctly account for order of edges (may be implemented through CheckEdgeOrder, FixEdgeOrder).
		void                        Connect                 (const HandleType * adjacent, INMOST_DATA_ENUM_TYPE num) const; 
		/// Update geometric data for element, calls RecomputeGeometricData from Mesh.
		void                        UpdateGeometricData     () const; 
	};
	
	__INLINE const Element & InvalidElement() {static Element ret(NULL,InvalidHandle()); return ret;}
	
	class Node : public Element //implemented in node.cpp
	{
	public:
		Node() : Element() {}
		Node(const Node & other) : Element(other) {}
		Node(Mesh * m, HandleType h) : Element(m,h) {}
		Node(Mesh * m, HandleType * h) : Element(m,h) {}
		Node & operator =(Node const & other) {Element::operator =(other); return *this;}
		Node * operator ->() {return this;}
		const Node * operator ->() const {return this;}
		Node & self() {return *this;}
		const Node & self() const {return *this;}
	public:
		ElementArray<Edge>          getEdges                () const; //unordered
		ElementArray<Face>          getFaces                () const; //unordered
		ElementArray<Cell>          getCells                () const; //unordered

		ElementArray<Edge>          getEdges                (MarkerType mask,bool invert_mask = false) const; //unordered
		ElementArray<Face>          getFaces                (MarkerType mask,bool invert_mask = false) const; //unordered
		ElementArray<Cell>          getCells                (MarkerType mask,bool invert_mask = false) const; //unordered

		Storage::real_array         Coords                  () const; 
	};

	__INLINE const Node & InvalidNode() {static Node ret(NULL,InvalidHandle()); return ret;}
	
	class Edge : public Element //implemented in edge.cpp
	{
	public:
		Edge() : Element() {}
		Edge(const Edge & other) : Element(other) {}
		Edge(Mesh * m, HandleType h) : Element(m,h) {}
		Edge(Mesh * m, HandleType * h) : Element(m,h) {}
		Edge & operator =(Edge const & other) {Element::operator =(other); return *this;}
		Edge * operator ->() {return this;}
		const Edge * operator ->() const {return this;}
		Edge & self() {return *this;}
		const Edge & self() const {return *this;}
	public:
		
		ElementArray<Node>          getNodes                () const; //ordered
		ElementArray<Face>          getFaces                () const; //unordered
		ElementArray<Cell>          getCells                () const; //unordered

		ElementArray<Node>          getNodes                (MarkerType mask,bool invert_mask = false) const; //ordered
		ElementArray<Face>          getFaces                (MarkerType mask,bool invert_mask = false) const; //unordered
		ElementArray<Cell>          getCells                (MarkerType mask,bool invert_mask = false) const; //unordered

		Node                        getBeg                  () const;
		Node                        getEnd                  () const;
		//implemented in modify.cpp
		static Edge                 UniteEdges              (ElementArray<Edge> & edges, MarkerType del_protect);
		static bool                 TestUniteEdges          (const ElementArray<Edge> & edges, MarkerType del_protect);
		static ElementArray<Edge>   SplitEdge               (Edge e, const ElementArray<Node> & nodes, MarkerType del_protect); //provide ordered array of nodes, that lay between former nodes of the edge
		static bool                 TestSplitEdge           (Edge e, const ElementArray<Node> & nodes, MarkerType del_protect);
		//implemented in geometry.cpp
		Storage::real               Length                  () const;
	};

	__INLINE const Edge & InvalidEdge() {static Edge ret(NULL,InvalidHandle()); return ret;}
	

	/// \brief An interface for elements of type FACE.
	///
	/// This interface carry the link to the mesh and the element's handle that
	/// represents position of the element's data in the mesh. 
	///
	/// Interface provides some operations that can be done uniquely on the object of class 
	/// Element for which Element::GetElementType retruns FACE.
	///
	/// For the basic set of operations on all of the elements check class Element.
	///
	/// For the basic set of operations on the data of the element check class Storage.
	///
	/// You can obtain object of class Face from Mesh::iteratorFace,
	/// in this case obtained object is always valid.
	/// Also you can get it through Mesh::FaceByLocalID, check with Element::isValid to see whether you
	/// have obtained a valid object.
	/// You can convert an object of class Element into an object of class Face by Element::getAsFace. In debug
	/// mode it will internally check that the element is of type FACE.
	/// You may compose an object of class Face by using constructor and specifing pointer to the mesh
	/// and providing element's handle. You can make a handle with ComposeHandle(FACE,local_id) function.
	class Face : public Element //implemented in face.cpp
	{
	public:
		Face() : Element() {}
		Face(const Face & other) : Element(other) {}
		Face(Mesh * m, HandleType h) : Element(m,h) {}
		Face(Mesh * m, HandleType * h) : Element(m,h) {}
		Face & operator =(Face const & other) {Element::operator =(other); return *this;}
		Face * operator ->() {return this;}
		const Face * operator ->() const {return this;}
		Face & self() {return *this;}
		const Face & self() const {return *this;}
	public:
		
		ElementArray<Node>          getNodes                () const; //ordered
		ElementArray<Edge>          getEdges                () const; //ordered
		ElementArray<Cell>          getCells                () const; //unordered

		ElementArray<Node>          getNodes                (MarkerType mask,bool invert_mask = false) const; //ordered
		ElementArray<Edge>          getEdges                (MarkerType mask,bool invert_mask = false) const; //ordered
		ElementArray<Cell>          getCells                (MarkerType mask,bool invert_mask = false) const; //unordered

		//this is for 2d case when the face is represented by segment
		Node                        getBeg                  () const;
		Node                        getEnd                  () const;
		/// \brief Retrieve the cell for which the normal points outwards.
		///
		/// \warning Depending on the grid construction the normal may incorrectly point inwards.
		/// You can resolve this situation by Face::FixNormalOrientation.
		///
		/// @return the cell for which normal points outwards.
		/// @see Face::FixNormalOrientation
		Cell                        BackCell                () const;
		/// \brief Retrieve the cell for which the normal points inwards.
		///
		/// \warning Depending on the grid construction the normal may incorrectly point outwards.
		/// You can resolve this situation by Face::FixNormalOrientation.
		///
		/// @return the cell for which normal points inwards.
		/// @see Face::FixNormalOrientation
		Cell                        FrontCell               () const;
		bool                        FaceOrientedOutside     (Cell c) const;
		void                        ReorderEdges            () const;
		bool                        CheckEdgeOrder          () const; //not implemented// returns true if edges of face form an ordered closed loop
		bool                        FixEdgeOrder            () const; //not implemented// returns true if edges were successfully reordered to form a closed loop
		//implemented in modify.cpp
		static Face                 UniteFaces              (ElementArray<Face> & faces, MarkerType del_protect);
		static bool                 TestUniteFaces          (const ElementArray<Face> & faces,  MarkerType del_protect);
		static ElementArray<Face>   SplitFace               (Face face, const ElementArray<Edge> & edges, MarkerType del_protect); //provide all edges that lay inside face
		static bool                 TestSplitFace           (Face face, const ElementArray<Edge> & edges, MarkerType del_protect);	
		void                        SwapCells               (); //swap back cell and front cell
		//implemented in geometry.cpp
		Storage::real               Area                    () const;
		void                        Normal                  (real * nrm) const;
		void                        UnitNormal              (real * nrm) const;
		void                        OrientedNormal          (Cell c, real * nrm) const;
		void                        OrientedUnitNormal      (Cell c, real * nrm) const;
		bool                        FixNormalOrientation    () const;  //returns true if orientation was corrected, otherwise returns false
		bool                        CheckNormalOrientation  () const; //returns true if orientation is correct, otherwise returns false
		bool                        Closure                 () const; // test integrity of polygon
	};

	__INLINE const Face & InvalidFace() {static Face ret(NULL,InvalidHandle()); return ret;}
	
	//implemented in cell.cpp
	/// \brief An interface for elements of type CELL.
	///
	/// This interface carry the link to the mesh and the element's handle that
	/// represents position of the element's data in the mesh. 
	///
	/// Interface provides some operations that can be done uniquely on the object of class 
	/// Element for which Element::GetElementType retruns CELL.
	///
	/// For the basic set of operations on all of the elements check class Element.
	///
	/// For the basic set of operations on the data of the element check class Storage.
	///
	/// You can obtain object of class Cell from Mesh::iteratorCell,
	/// in this case obtained object is always valid.
	/// Also you can get it through Mesh::CellByLocalID, check with Element::isValid to see whether you
	/// have obtained valid object.
	/// You can convert an object of class Element into an object of class Cell by Element::getAsCell. In debug
	/// mode it will internally check that the element is of type CELL.
	/// You may compose an object of class Cell by using constructor and specifing pointer to the mesh
	/// and providing element's handle. You can make a handle with ComposeHandle(CELL,local_id) function.
	class Cell : public Element
	{
	public:
		/// \brief Basic constructor.
		///
		/// Constructs invalid cell.
		Cell() : Element() {} 
		/// \brief Basic constructor with fixed handle.
		///
		/// When constructed this way, only handle within object may change.
		///
		/// @param m Pointer to the mesh to which the element belongs.
		/// @param h Unique handle that describes position of the element within the mesh.
		Cell(Mesh * m, HandleType h) : Element(m,h) {}
		/// \brief Basic constructor with an assignable handle.
		///
		/// When constructed this way, the provided memory location for handle
		/// will be modified on assignment.
		///
		/// The purpose of this function is to be used in various non-constant iterators
		/// of containers that allow underlaying contents to be changed.
		///
		/// @param m Pointer to the mesh to which the element belongs.
		/// @param h Pointer to the handle that describes position of the element within mesh.
		Cell(Mesh * m, HandleType * h) : Element(m,h) {}
		/// \brief Copy constructor.
		///
		/// New object will inherit the assignment behavior of the initial object.
		///
		/// @param other Object of type Cell to be duplicated.
		Cell(const Cell & other) : Element(other) {}
		/// \brief Assignment operator
		///
		/// Assigned object will inherit the asignment behavior of the initial object only
		/// in case it was constructed without pointer to a handle. Otherwise it will
		/// only modify the memory location to which it points.
		///
		/// @param other Object of type Cell to be duplicated.
		/// @return Reference to the current object.
		Cell & operator =(Cell const & other) {Element::operator =(other); return *this;} ///< Assignment Operator.
		/// \brief Operator of dereference to pointer.
		///
		/// This is needed for iterators to work properly.
		///
		/// @return Pointer to the current object.
		Cell * operator ->() {return this;} 
		/// \brief Operator of dereference to constant pointer. 
		///
		/// This is needed for const_iterators to work properly.
		///
		/// @return Constant pointer to the current object.
		const Cell * operator ->() const {return this;}
		/// \brief Get self-reference.
		///
		/// Main purpose is to convert iterators into elements and to pass them as arguments of functions.
		///
		/// @return Constant reference to the current object.
		Cell & self() {return *this;} 
		/// \brief Get constant self-reference.
		///
		/// Main purpose is to convert iterators into elements and to pass them as arguments of functions.
		///
		/// @return Reference to the current object.
		const Cell & self() const {return *this;}
	public:
		/// \brief Get all the nodes of the current cell.
		///
		/// This function traverses up the adjacency graph by one level.
		///
		/// The order of nodes will be preserved from the moment of the construction of the cell.
		/// When suggest_nodes array is not supplied into the Mesh::CreateCell functions, then
		/// for known types the order of nodes follows VTK convention. For polyhedral cells
		/// the order is unspecified. When suggest_nodes_order was provided into Mesh::CreateCell
		/// then the order follows provided order.
		///
		/// @return Set of nodes that compose current cell.
		ElementArray<Node>          getNodes                () const;
		/// \brief Get all the edges of the current cell.
		///
		/// This function traverses down the adjacency graph by two levels.
		///
		/// The order of the edges is unspecified.
		///
		/// @return Set of edges that compose current cell.
		///
		/// \todo
		/// One should thoroughly check three scenarios of function execution in shared parallel environment
		/// for different types of cells (simple tet/hex cells as well as complex polyhedral cells) and
		/// draw a conclusion on the best scenario for each condition. One of the development versions
		/// contains all the algorithms, ask for the files.
		///  1. Use of markers (current wariant).
		///  2. Put all elements into array with duplications, then run std::sort and std::unique.
		///  3. Put all elements into array, check for duplication by running through array.
		///
		/// \warning Note that this function uses markers to check for duplication of edges
		/// in output array. This allows for faster algorithm, but will lead to penalties
		/// in shared parallel execution, since operations on markers are declared atomic.
		/// For atomic declaration to work you have to define USE_OMP during CMake configuration
		/// or in inmost_common.h. If you attempt to use this function in shared parallel execution
		/// without USE_OMP you may expect side effects.
		ElementArray<Edge>          getEdges                () const;
		/// \brief Get all the faces of the current cell.
		///
		/// This function traverses down the adjacency graph by one level.
		///
		/// The order of the faces returned by this function is preserved from 
		/// the moment of the construction of the cell.
		///
		/// @return Set of faces that compose current cell.
		ElementArray<Face>          getFaces                () const;
		/// \brief Get the subset of the nodes of the current cell that are (not) marked by provided marker.
		///
		/// This function traverses up the adjacency graph by one level.
		///
		/// The order of nodes will be preserved from the moment of the construction of the cell.
		/// When suggest_nodes array is not supplied into the Mesh::CreateCell functions, then
		/// for known types the order of nodes follows VTK convention. For polyhedral cells
		/// the order is unspecified. When suggest_nodes_order was provided into Mesh::CreateCell
		/// then the order follows provided order.
		///
		/// @param mask Marker that should be used to filter elements.
		/// @param invert_mask If false then those elements that are marked will be taken, otherwise
		///                    elements that are not marked will be taken.
		/// @return Set of the nodes that compose current cell and (not) marked by marker.
		///
		/// \warning To work correctly in shared parallel environment this function requires
		/// USE_OMP to be defined in CMake or in inmost_common.h.
		ElementArray<Node>          getNodes                (MarkerType mask,bool invert_mask = false) const;
		/// \brief Get the subset of the edges of the current cell that are (not) marked by provided marker.
		///
		/// This function traverses down the adjacency graph by two levels.
		///
		/// The order of the edges is unspecified.
		///
		/// @param mask Marker that should be used to filter elements.
		/// @param invert_mask If false then those elements that are marked will be taken, otherwise
		///                    elements that are not marked will be taken.
		/// @return Set of the edges that compose current cell and (not) marked by marker.
		///
		/// \warning To work correctly in shared parallel environment this function requires
		/// USE_OMP to be defined in CMake or in inmost_common.h.
		ElementArray<Edge>          getEdges                (MarkerType mask,bool invert_mask = false) const;
		/// \brief Get the subset of the faces of the current cell that are (not) marked by provided marker.
		///
		/// This function traverses down the adjacency graph by one level.
		///
		/// The order of the faces returned by this function is preserved from 
		/// the moment of the construction of the cell.
		///
		/// @param mask Marker that should be used to filter elements.
		/// @param invert_mask If false then those elements that are marked will be taken, otherwise
		///                    elements that are not marked will be taken.
		/// @return Set of the faces that compose current cell and (not) marked by marker.
		///
		/// \warning To work correctly in shared parallel environment this function requires
		/// USE_OMP to be defined in CMake or in inmost_common.h.
		ElementArray<Face>          getFaces                (MarkerType mask,bool invert_mask = false) const;
		/// \brief Check that sequence of edges form a closed loop and each edge have a node that matches one of the nodes at the next edge.
		///
		/// This function works for cells for which Element::GetElementDimension returns 2.
		///
		/// @return True if edges form the correct closed loop.
		///
		/// \todo
		///   1. Implement.
		///   2. Use in topology check algorithms.
		///   3. Use in Element::Connect.
		bool                        CheckEdgeOrder          () const; //not implemented//2D only, returns true if edges of face form an ordered closed loop
		/// \brief Repair the sequence of edges so that each edge have node that matches one of the nodes at the next edge.
		///
		/// This function works for cells for which Element::GetGeometricDimension returns 2.
		///
		/// This function may fail if all the edges does not form a closed loop.
		///
		/// @return True if edges were successfully reordered to form a closed loop.
		///
		/// \todo
		///   1. Implement.
		///   2. Use in topology check algorithms.
		///   3. Use in Element::Connect.
		bool                        FixEdgeOrder            () const; //not implemented//2D only, returns true if edges were successfully reordered to form a closed loop
		//implemented in modify.cpp
		/// \brief Unite a set of given cells into one cell.
		///
		/// This will create a cell whose faces are formed by symmetric difference
		/// of faces of given cells.
		/// If you specify a nonzero marker then the procedure will fail 
		/// if any marked element have to be deleted during union.
		/// @param cells A set of cells to be united.
		/// @param del_protect A marker that protects elements from deletion. Zero means no check.
		/// @return A new cell.
		static Cell                 UniteCells              (ElementArray<Cell> & cells, MarkerType del_protect);
		/// \brief Test that no marked element will be deleted during union of given cells.
		/// @param cells A set of cells to be united.
		/// @param del_protect A marker that protects elements from deletion. Zero means no check.
		/// @return True if no element deleted, false otherwise.
		static bool                 TestUniteCells          (const ElementArray<Cell> & cells, MarkerType del_protect);
		/// \brief Separate a cell according to the set of provided faces.
		///
		/// You should first separate all edges with new nodes by Edge::SplitEdge, then faces with new edges by Face::SplitFace
		/// and then create faces using new edges from faces of initial cell and optionally edges from inside of the cell. 
		/// All faces are supposed to be internal with respect to the current cell. 
		/// Internally this function will resolve geometry of new cells inside of the 
		/// current cell by running through adjacency graph.
		///
		/// @param cell A cell to be split.
		/// @param faces A set of faces, internal for current cell, that will be used for construction of new cells.
		/// @param del_protect Marker that may be used to protect some elements from deletion by algorithm. 
		///        Zero means no check.
		/// @return A set of new cells.
		/// @see Edge::SplitEdge
		/// @see Face::SplitFace
		///
		/// \todo
		///   1. The algorithm inside is minimizing the size of the adjacency graph for each new cell.
		///      The correct behavior is to calculate volume of the cell for each adjacency graph and choose the graph with minimal volume.
		///      This requires calculation of volume for non-convex cells. For correct calculation of volume on non-convex cells one should 
		///      find one face for which normal orientation can be clearly determined and then orient all edges of the cell with respect to
		///      the orientation of edges of this face and establish normals for all faces. Once the algorithm is implemented here it should 
		///      be implemented in geometrical services or vice verse.
		///   2. Probably the algorithm should minimize the volume and adjacency graph size altogether. Between the cells with smallest volume
		///      within some tolerance select those that have smallest adjacency graph.
		static ElementArray<Cell>   SplitCell               (Cell cell, const ElementArray<Face> & faces, MarkerType del_protect); //provide all faces, that lay inside cell
		/// \brief This functions checks is it possible to split the cell by the given set of faces without deleting marked elements.
		///
		/// @param cell A cell to be split.
		/// @param faces Set of faces, internal for current cell, that will be used for new cells.
		/// @param del_protect Marker that may be used to protect some elements from deletion by algorithm. Zero means no check.
		/// @return True if possible, otherwise false.
		static bool                 TestSplitCell           (Cell cell, const ElementArray<Face> & faces, MarkerType del_protect);
		//implemented in geometry.cpp
		/// \brief Get a cell that share a face with the current cell.
		///
		/// Don't forget to check that the returned cell is valid.
		/// It would return invalid cell for boundary face.
		///
		/// @param face A face of current cell for which neighbouring cell is determined.
		/// @return A cell that shares a given face with the current cell.
		Cell                        Neighbour               (Face face) const;
		/// \brief Get all cells that share the face with the current cell.
		/// @return Set of cells that share a face with the current cell.
		ElementArray<Cell>          NeighbouringCells       () const; // get all cells that share any face with current
		/// \brief Determine, if point lies inside element.
		///
		/// Now it works only for 3-dimensional elements, at future it will be
		/// extended to support polygons and segments.
		/// @param point coordinates of the point, it is assumed that number of the
		/// coordinates is the same as the number of dimensions in the mesh.
		/// @return returns true if points inside element, false otherwise
		/// @see Mesh::GetDimensions
		///
		/// \todo
		///   1. Should be checked or extended for 2d cells. (done, testing)
		bool                        Inside                  (const real * point) const; //is point inside cell, check for 2d case
		/// \brief Return volume of the cell.
		///
		/// Note that currently the volume for non-convex cells may be calculated incorrectly.
		/// \todo
		///   1. Geometric services should correctly resolve volume for non-convex cells.
		real                        Volume                  () const;
		/// \brief Test that faces of the cell form the closed set.
		///
		/// This is automatically checked for if you activate NEED_TEST_CLOSURE
		/// in Mesh::SetTopologyCheck.
		/// @return True if faces of the cell form the closed set, false otherwise.
		bool                        Closure                 () const; // test integrity of cell
	};

	///
	__INLINE const Cell & InvalidCell() {static Cell ret(NULL,InvalidHandle()); return ret;}

	class ElementSet : public Element //implemented in eset.cpp
	{
	public:
    ///Number of reserved positions in HighConn array.
    ///The first position is the handle to parent set.
    ///The second position is handle to sibling set.
    ///The third position is handle to child set.
    ///The fourth position is number of sorted elements in the set.
    ///All the rest are positions of deleted elements.
		static const enumerator             high_conn_reserved  = 4;
    __INLINE static HandleType &        hParent             (Element::adj_type & arr) {return arr[0];}
	  __INLINE static HandleType &        hSibling            (Element::adj_type & arr) {return arr[1];}
	  __INLINE static HandleType &        hChild              (Element::adj_type & arr) {return arr[2];}
	  __INLINE static HandleType &        hSorted             (Element::adj_type & arr) {return arr[3];}
		typedef INMOST_DATA_BULK_TYPE       ComparatorType;
		static const ComparatorType         UNSORTED_COMPARATOR = 0;
		static const ComparatorType         GLOBALID_COMPARATOR = 1;
		static const ComparatorType         CENTROID_COMPARATOR = 2;
		static const ComparatorType         HIERARCHY_COMPARATOR  = 3;
		static const ComparatorType         HANDLE_COMPARATOR   = 4;
		                                    ElementSet          () : Element() {}
		                                    ElementSet          (Mesh * m, HandleType h) : Element(m,h) {}
		                                    ElementSet          (Mesh * m, HandleType * h) : Element(m,h) {}
		                                    ElementSet          (const ElementSet & other) : Element(other) {}
		__INLINE ElementSet &               operator =          (ElementSet const & other) {Element::operator =(other); return *this;}
		__INLINE ElementSet *               operator ->         () {return this;}
		__INLINE const ElementSet *         operator ->         () const {return this;}
		__INLINE ElementSet &               self                ()              {return *this;}
		__INLINE const ElementSet &         self                () const        {return *this;}
	public:
		/// Get name of the set
		std::string                 GetName() const;
		/// Retrieve parent of the set
		ElementSet                  GetParent() const;
		/// Retrieve sibling set of the set, this will be next child for the parent
		ElementSet                  GetSibling() const;
		/// Retrieve child set of the set.
		/// Following example shows how to iterate through children of current set:
		/// for(ElementSet it = set->GetChild(); it->isValid(); it = it->GetSibling()) ...
		ElementSet                  GetChild() const;
		
		/// This will create new child for the parent
		void                        AddSibling(const ElementSet & sibling) const;
		/// Add child to current set
		void                        AddChild(const ElementSet & child) const;

		/// This will erase sibling or parent's child
		void                        RemSibling(const ElementSet & sibling) const;
		/// This will erase my child
		void                        RemChild(const ElementSet & child) const;
		/// How many there are siblings to the right of me including me
		enumerator                  CountSiblings() const;
		/// How many children I have 
		enumerator                  CountChildren() const;

		bool                        HaveSibling() const;
		bool                        HaveParent() const;
		bool                        HaveChild() const;

		/// Direct raw access to stored elements, no copy involved.
		/// The actual pointer may change when you add new elements due to reallocation of memory.
		/// Modify at your own risk.
		HandleType *                getHandles() const;
		/// Retrieve number of stored handles, including invalid.
		/// If you want to get number of valid elements use ElementSet::Size
		enumerator                  nbHandles() const;
		/// Retrieve position after last sorted element. This does not correctly
		/// represent total number of sorted elements, since some of them may be deleted
		enumerator                  nbSorted() const;
		/// Retrieve all elements by type
		enumerator                  nbAdjElements(ElementType etype) const;
		enumerator                  nbAdjElements(ElementType etype, MarkerType select, bool invert = false) const;
		/// Retrieve all elements by type
		ElementArray<Element>       getAdjElements(ElementType etype) const;
		ElementArray<Element>       getAdjElements(ElementType etype, MarkerType select, bool invert = false) const;
		
		/// Retrieve only nodes
		ElementArray<Node>          getNodes() const;
		ElementArray<Node>          getNodes(MarkerType select, bool invert = false) const;
		/// Retrieve only edges
		ElementArray<Edge>          getEdges() const;
		ElementArray<Edge>          getEdges(MarkerType select, bool invert = false) const;
		/// Retrieve only faces
		ElementArray<Face>          getFaces() const;
		ElementArray<Face>          getFaces(MarkerType select, bool invert = false) const;
		/// Retrieve only cells
		ElementArray<Cell>          getCells() const;
		ElementArray<Cell>          getCells(MarkerType select, bool invert = false) const;
		/// Put one element without checking of the existance of duplicate.
		/// For sorted set new element will appear at unsorted part.
		void                        PutElement(HandleType e) const;
		/// Put one element without checking of the existance of duplicate.
		/// For sorted set new element will appear at unsorted part
		void                        PutElement(const Storage & e) const {PutElement(e->GetHandle());}
		/// Put multiple handles without checking of the existance of duplicate
		void                        PutElements(const HandleType * handles, enumerator num) const;
		/// Put multiple handles of the other set without checking of the existance of duplicate
		void                        PutElements(const ElementSet & other) const {PutElements(other->getHandles(),other->nbHandles());}
		/// Put multiple handles without checking
		template<typename EType>
		void                        PutElements(const ElementArray<EType> & elems) const {PutElements(elems.data(),static_cast<enumerator>(elems.size()));}
		/// Put one element with checking of the existance of duplicate.
		/// Preserves order for sorted set, thus may be expensive
		void                        AddElement(HandleType e) const;
		/// Put one element with checking of the existance of duplicate.
		/// Preserves order for sorted set, thus may be expensive
		void                        AddElement(const Storage & e) const {AddElement(e->GetHandle());}
		/// Add multiple elements with checking of the existance of duplicate.
		/// Preserves order for sorted set, thus may be expensive.
		/// This will also remove any duplicates in unsorted part of the set.
		/// If you inserted duplicated elements through PutElements into previously sorted array
		/// then this operation does not guarantee that those duplications will be stored.
    /// \todo Recheck usage of markers.
		void                        AddElements(const HandleType * handles, enumerator num) const;
		/// Add elements of other set
		void                        AddElements(const ElementSet & other) {Unite(other);}
		/// Add multiple elements with checking of the existance of duplicate
		template<typename EType>
		void                        AddElements(const ElementArray<EType> & elems) const {AddElements(elems.data(),static_cast<enumerator>(elems.size()));}

		void                        RemoveElement(const Storage & e) const;
		void                        RemoveElements(const HandleType * handles, enumerator num) const;
		/// Remove multiple elements from the set
		template<typename EType>
		void                        RemoveElements(const ElementArray<EType> & elems) const {RemoveElements(elems.data(),static_cast<enumerator>(elems.size()));}

		/// Compute and return union with other set.
		/// Result is unordered.
		/// All elements will be unique.
		ElementArray<Element> Union(const ElementSet & other) const;
		/// Compute and return union with raw handles.
		/// Result is unordered.
		/// All elements will be unique.
		ElementArray<Element> Union(const HandleType * handles, enumerator num) const;
		/// Compute and return union with elements
		template<typename EType>
		ElementArray<Element> Union(const ElementArray<EType> & elems) const {return Union(elems.data(),static_cast<enumerator>(elems.size()));}

		ElementArray<Element> Difference(const ElementSet & other) const;
		/// Compute and return difference with raw handles.
		/// If initial set was ordered, result will preserve the order.
		/// All elements will be unique.
		ElementArray<Element> Difference(const HandleType * handles, enumerator num) const;
		/// Compute and return difference with elements
		template<typename EType>
		ElementArray<Element> Difference(const ElementArray<EType> & elems) const {return Difference(elems.data(),static_cast<enumerator>(elems.size()));}

		ElementArray<Element> Intersection(const ElementSet & other) const;
		/// Compute and return intersection with raw handles.
		/// If initial set was ordered, result will preserve the order.
		/// All elements will be unique.
		ElementArray<Element> Intersection(const HandleType * handles, enumerator num) const;
		/// Compute and return intersection with elements.
		template<typename EType>
		ElementArray<Element> Intersection(const ElementArray<EType> & elems) const {return Intersection(elems.data(),static_cast<enumerator>(elems.size()));}
		
		/// Compute and store union with raw handles.
		void Unite(const ElementSet & other) const;
		/// Compute and store union with raw handles.
		void Unite(const HandleType * handles, enumerator num) const;
		/// Compute and store union with elements.
		template<typename EType>
		void Unite(const ElementArray<EType> & elems) const {Unite(elems.data(),static_cast<enumerator>(elems.size()));}

		/// Compute and store difference with raw handles.
		/// \todo
		///   If other and current sets are sorted in same way, may perform narrowing traversal by retriving
		///   mutual lower_bound/higher_bound O(log(n)) operations for detecting common subsets in sorted sets.
		///   May work good when deleting handles by small chunks, ApplyModification may greatly benefit.
		void Subtract(const ElementSet & other) const;
		/// Compute and store difference with raw handles
		void Subtract(const HandleType * handles, enumerator num) const;
		/// Compute and store difference with elements
		template<typename EType>
		void Subtract(const ElementArray<EType> & elems) const {Subtract(elems.data(),static_cast<enumerator>(elems.size()));}


		/// Compute and store intersection with raw handles
		void Intersect(const ElementSet & other) const;
		/// Compute and store intersection with raw handles
		void Intersect(const HandleType * handles, enumerator num) const;
		/// Compute and store intersection with elements
		template<typename EType>
		void Intersect(const ElementArray<EType> & elems) const {Intersect(elems.data(),static_cast<enumerator>(elems.size()));}
		/// Performs sort of the set of elements. If the set was previously sorted but
		/// have unsorted part, then unsorted part will be sorted and two parts will be merged.
		/// If you need all the set to be resorted (for example in the case global ids were changed)
		/// then invoke SortSet with UNSORTED_COMPARATOR first and then with needed comparator.
		///
		/// Internally it uses:
		/// - std::sort with Mesh::CentroidComparator for CENTROID_COMPARATOR
		/// - std::sort with Mesh::IerarhyComparator  for  HIERARCHY_COMPARATOR
		/// - radix sort starting from certain size for   GLOBALID_COMPARATOR
		/// - radix sort starting from certain size for     HANDLE_COMPARATOR
		/// - only changes state, no sort performed for   UNSORTED_COMPARATOR 
		///
		/// After the set was sorted all the invalid handles should appear at the end of the set
		/// and then removed from array, so it will implicitly work as ReorderEmpty function.
		/// No checks that elements are hidden performed (Maybe this checks should be done in comparators)
		/// In the case you formed the set by running over all mesh elements from NODE to ESET in 
		/// increasing order then your set will be automatically sorted by handles, in this case
		/// you are encouraged to override Mesh::SetComparatorTag with HANDLE_COMPARATOR on the set
		/// without invoking SortSet, so that SortSet does not redundantly perform any work.
		/// You are encouraged to do so even if you are not going to use this information -
		/// some internal algorithms may benefit from it.
		///
		/// \todo
		/// !TODO 52 - check radix sort on big endian computer
		/// @param comp one of the comparators from description
		/// @see Mesh::SetComparatorTag
		void SortSet(ComparatorType comp) const;
		/// Perform binary search by global id. In set sorted with GLOBALID_COMPARATOR in O(log(n)) time
		/// otherwise search needs O(n) comparisons
		Element FindElementByGlobalID(integer global_id) const;
		/// Perform binary search by centroid. In set sorted with CENTROID_COMPARATOR in O(log(n)) time
		/// otherwise search needs O(n) comparisons
		Element FindElementByCentroid(real * centroid) const;
		/// Performs linear search in unsorted set. 
		/// - binary search by handle in set sorted with      HANDLE_COMPARATOR
		/// - binary search by centroid in set sorted with  CENTROID_COMPARATOR
		/// - binary search by global id in set sorted with GLOBALID_COMPARATOR
		/// - binary search by ierarhy in set sorted with    HIERARCHY_COMPARATOR (may be too expensive)
		///
		/// If you have a lot of elements to test against set, then you'd better first put marker
		/// on elements of the set by SetMarkerElements then test which of your elements have marker
		/// and then remove markers by RemMarkerElements. This will consume only O(n+m) operations
		/// where n is the number of my elements and m is the size of the test set.
		/// @param h handle of element
		/// @param use_comparator use information about current comparator or perform linear search
		/// @return true if element is present
		/// @see ElementSet::SetMarkerElements
		/// @see ElementSet::RemMarkerElements
		bool FindHandle(HandleType h, bool use_comparator) const;
		/// Set markers on all the elements of given type
		void SetMarkerElements(MarkerType m, ElementType etype = ESET|CELL|FACE|EDGE|NODE) const;
    void SetPrivateMarkerElements(MarkerType m, ElementType etype = ESET|CELL|FACE|EDGE|NODE) const;
		/// Remove markers from all the elements of given type
		void RemMarkerElements(MarkerType m, ElementType etype = ESET|CELL|FACE|EDGE