inmost_mesh.h

#ifndef INMOST_MESH_H_INCLUDED
#define INMOST_MESH_H_INCLUDED

#include "inmost_common.h"

#if defined(USE_MESH)

#define __NDT 3
#define __NET 6

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 ElementType;
	static const ElementType NONE = 0x00;
	static const ElementType NODE = 0x01;
	static const ElementType EDGE = 0x02;
	static const ElementType FACE = 0x04;
	static const ElementType CELL = 0x08;
	static const ElementType ESET = 0x10;
	static const ElementType MESH = 0x20;

	typedef INMOST_DATA_ENUM_TYPE HandleType;
	static const INMOST_DATA_ENUM_TYPE handle_etype_bits  = 3;
	static const INMOST_DATA_ENUM_TYPE handle_etype_shift = sizeof(HandleType)*8-handle_etype_bits;
	static const INMOST_DATA_ENUM_TYPE handle_id_mask     = (1 << handle_etype_shift)-1;
	
	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;

	typedef array<INMOST_DATA_REAL_TYPE>    inner_real_array;
	typedef array<INMOST_DATA_INTEGER_TYPE> inner_integer_array;
	typedef array<INMOST_DATA_BULK_TYPE>    inner_bulk_array;
	typedef array<HandleType>               inner_reference_array;

	enum DataType
	{
		DATA_REAL      = 0, 
		DATA_INTEGER   = 1, 
		DATA_BULK      = 2,
		DATA_REFERENCE = 3
	};

	typedef INMOST_DATA_ENUM_TYPE MarkerType; // low 8 bits - marker mask, rest high bits - position of marker
	static const INMOST_DATA_ENUM_TYPE MarkerFields = 16;   // number of chars to hold all markers, total number (MarkerFields * bits_per_char)
	static const INMOST_DATA_ENUM_TYPE MarkerMask   = static_cast<INMOST_DATA_BULK_TYPE>(-1); // bit mask to obtain marker mask within MarkerType
	static const INMOST_DATA_ENUM_TYPE MarkerShift  = sizeof(INMOST_DATA_BULK_TYPE)*8;    // sizeof(char) * bits_per_char
	__INLINE static MarkerType InvalidMarker() {return ENUMUNDEF;}


	static const INMOST_DATA_ENUM_TYPE chunk_bits_elems       = 13;
	static const INMOST_DATA_ENUM_TYPE chunk_bits_empty       = 8;
	static const INMOST_DATA_ENUM_TYPE chunk_bits_tags        = 6;
	static const INMOST_DATA_ENUM_TYPE chunk_bits_dense_array = 6;

	//Use topolgy checking for debug purposes
	typedef INMOST_DATA_ENUM_TYPE TopologyCheck; 
	static const TopologyCheck THROW_EXCEPTION        = 0x00000001; //done//throw TopologyError exception on error
	static const TopologyCheck PRINT_NOTIFY           = 0x00000002; //done//print topology notify to std::cerr
	static const TopologyCheck DELETE_ON_ERROR        = 0x00000004; //done//element should be deleted if there is an error in EndTopologyCheck
	static const TopologyCheck MARK_ON_ERROR          = 0x00000008; //done//bitwise OR of errors is recorded to sparse integer tag of element (if not deleted), availible by TopologyErrorTag()
	static const TopologyCheck DUPLICATE_EDGE         = 0x00000010; //done//check that edge already exists, on occurance CreateEdge silently returns existant edge
	static const TopologyCheck DUPLICATE_FACE         = 0x00000020; //done//check that face already exists, on occurance CreateFace silently returns existant face
	static const TopologyCheck DUPLICATE_CELL         = 0x00000040; //done//check that cell already exists, on occurance CreateCell silently returns existant cell
	static const TopologyCheck DEGENERATE_EDGE        = 0x00000080; //done//check that edge have more then two nodes (no support for this kind of object)
	static const TopologyCheck DEGENERATE_FACE        = 0x00000100; //done//produce error if face consists of less then 3 edges
	static const TopologyCheck DEGENERATE_CELL        = 0x00000200; //done//produce error if cell consists of less then 4 faces
	static const TopologyCheck FACE_ORIENTATION       = 0x00000400; //done//produce error if face have wrong orientation of edges, for star-shaped elements only
	static const TopologyCheck FACE_PLANARITY         = 0x00000800; //done//produce error if face is non planar
	static const TopologyCheck INTERLEAVED_FACES      = 0x00001000; //done//produce error if there is another face that already uses same nodes
	static const TopologyCheck TRIPLE_SHARED_FACE     = 0x00002000; //done//check that every face have exectly two neighbours, if DUPLICATE_CELL is activated, then checks for cell duplication
	static const TopologyCheck FLATTENED_CELL         = 0x00004000; //done//produce error if one of the faces of the cell contains all the nodes of the cell
	static const TopologyCheck ADJACENT_DUPLICATE     = 0x00008000; //done//produce error if provided array of elements for construction contain duplications
	static const TopologyCheck ADJACENT_HIDDEN        = 0x00010000; //done//hidden elements should not be used when new elements are created
	static const TopologyCheck ADJACENT_VALID         = 0x00020000; //done//check that all handles are valid
	static const TopologyCheck ADJACENT_DIMENSION     = 0x00040000; //done//produce error if provided array of elements have wrong geometric dimension
	static const TopologyCheck PROHIBIT_MULTILINE     = 0x00080000; //done//(needs NEED_TEST_CLOSURE) produce error if edges of faces are not closed, or face is folded
	static const TopologyCheck PROHIBIT_POLYGON       = 0x00100000; //done//allow only known types of elements: Tet,Quad, Line
	static const TopologyCheck PROHIBIT_MULTIPOLYGON  = 0x00200000; //done//(needs NEED_TEST_CLOSURE) produce error if faces of cell are not closed, or cell is folded
	static const TopologyCheck PROHIBIT_POLYHEDRON    = 0x00400000; //done//allow only known types of elements: Tet,Hex,Prism,Pyramid
	static const TopologyCheck FACE_EDGES_ORDER       = 0x00800000; //done//edges of the face should form one closed loop
	static const TopologyCheck PROHIBIT_CONCAVE_FACE  = 0x01000000; //not implemented//don't allow concave face
	static const TopologyCheck PROHIBIT_CONCAVE_CELL  = 0x02000000; //not implemented//don't allow concave cell
	static const TopologyCheck PROHIBIT_NONSTAR_FACE  = 0x04000000; //not implemented//don't allow non-star shaped face
	static const TopologyCheck PROHIBIT_NONSTAR_CELL  = 0x08000000; //not implemented//don't allow non-star shaped concave cell
	static const TopologyCheck FACE_SELF_INTERSECTION = 0x10000000; //not implemented//edges of the face don't cross each other
	static const TopologyCheck CELL_SELF_INTERSECTION = 0x20000000; //not implemented//faces of the cell don't cross each other
	static const TopologyCheck NEED_TEST_CLOSURE      = 0x40000000; //done//silent, test's for closure in ComputeGeometricType, needed to detect MultiLine and MultiPolygon
	static const TopologyCheck DISABLE_2D             = 0x80000000; //done//don't allow 2d grids, where edges appear to be vertexes, faces are edges and cells are faces
	static const TopologyCheck GRID_CONFORMITY        = NEED_TEST_CLOSURE | PROHIBIT_MULTILINE | PROHIBIT_MULTIPOLYGON  | INTERLEAVED_FACES | TRIPLE_SHARED_FACE;
	static const TopologyCheck DEFAULT_CHECK          = THROW_EXCEPTION | DUPLICATE_EDGE | DUPLICATE_FACE | DUPLICATE_CELL | PRINT_NOTIFY;
	const char *                             TopologyCheckNotifyString(TopologyCheck c); //mesh.cpp
	const char *                             DataTypeName         (DataType t); //tag.cpp
	__INLINE bool                            OneType              (ElementType t) {return t > 0 && (t & (t-1)) == 0;}
	__INLINE ElementType                     FirstElementType     () {return NODE;}
	__INLINE ElementType                     LastElementType      () {return MESH << 1;}
	__INLINE ElementType                     NextElementType      (ElementType etype) {return etype << 1;}
	__INLINE ElementType                     PrevElementType      (ElementType etype) {return etype >> 1;}
	__INLINE ElementType                     ElementTypeFromDim   (INMOST_DATA_INTEGER_TYPE dim) {return 1 << dim;}
	const char *                             ElementTypeName      (ElementType t); //mesh.cpp
	__INLINE INMOST_DATA_INTEGER_TYPE        ElementNum           (ElementType t)
	{
		unsigned int v = static_cast<unsigned int>(t);  // 32-bit value to find the log2 of 
		//static const unsigned int b[] = {0xAAAAAAAA, 0xCCCCCCCC, 0xF0F0F0F0, 0xFF00FF00, 0xFFFF0000};
		static const unsigned int b[] = {0xAAAAAAAA, 0xCCCCCCCC, 0xF0F0F0F0};
		register unsigned int r = (v & b[0]) != 0;
		//r |= ((v & b[4]) != 0) << 4;
		//r |= ((v & b[3]) != 0) << 3;
		r |= ((v & b[2]) != 0) << 2;
		r |= ((v & b[1]) != 0) << 1;
		return static_cast<INMOST_DATA_INTEGER_TYPE>(r);
	}
	__INLINE HandleType                      InvalidHandle        () {return 0;}
	__INLINE INMOST_DATA_INTEGER_TYPE        GetHandleID          (HandleType h) {return (h & handle_id_mask)-1;}
	__INLINE INMOST_DATA_INTEGER_TYPE        GetHandleElementNum  (HandleType h) {return h >> handle_etype_shift;}
	__INLINE ElementType                     GetHandleElementType (HandleType h) {return 1 << GetHandleElementNum(h);}
	__INLINE HandleType                      ComposeHandle        (ElementType etype, INMOST_DATA_INTEGER_TYPE ID) {return ID == -1 ? InvalidHandle() : ((ElementNum(etype) << handle_etype_shift) + (1+ID));}
	__INLINE HandleType                      ComposeHandleNum     (INMOST_DATA_INTEGER_TYPE etypenum, INMOST_DATA_INTEGER_TYPE ID) {return ID == -1 ? InvalidHandle() : ((etypenum << handle_etype_shift) + (1+ID));}
	__INLINE bool                            isValidHandle        (HandleType h) {return h != 0;}
	

	

	class TagMemory //implemented in tag.cpp
	{
	public:
		~TagMemory();
		TagMemory(Mesh * m, const TagMemory & other);
		TagMemory & operator =(TagMemory const & other);
	private:
		TagMemory();
		INMOST_DATA_ENUM_TYPE  pos[__NET]; 
		DataType               dtype;                 
		std::string            tagname;            
		INMOST_MPI_Type        bulk_data_type;   
		INMOST_DATA_ENUM_TYPE  size, bytes_size;
		bool                   sparse[__NET];
		INMOST_DATA_ENUM_TYPE  record_size;
		Mesh *                 m_link;
		friend class Tag;
		friend class Storage; //for debug
	};

	class Tag //implemented in tag.cpp
	{
	private:
		TagMemory * mem;
		Tag(Mesh * m, std::string name, DataType _dtype, INMOST_DATA_ENUM_TYPE size); //Create me through TagManager::CreateTag
		__INLINE INMOST_DATA_ENUM_TYPE  GetRecordSize() const {return mem->record_size;}
		__INLINE void                   SetSize(INMOST_DATA_ENUM_TYPE size) {mem->size = size;}
		__INLINE void                   SetPosition(INMOST_DATA_ENUM_TYPE pos, ElementType type) {mem->pos[ElementNum(type)] = pos;}
		__INLINE INMOST_DATA_ENUM_TYPE  GetPosition(ElementType type) const {assert(mem != NULL); return mem->pos[ElementNum(type)];}
		__INLINE void                   SetSparse(ElementType type) {mem->sparse[ElementNum(type)] = true;}
		__INLINE INMOST_DATA_ENUM_TYPE  GetPositionByDim(INMOST_DATA_ENUM_TYPE typenum) const {return mem->pos[typenum];}
	public:
		~Tag() {mem = NULL;}
		Tag() {mem = NULL;}
		Tag(const Tag & other) {mem = other.mem;}
		bool operator <(const Tag & other) const {return mem < other.mem;}
		bool operator >(const Tag & other) const {return mem > other.mem;}
		bool operator ==(const Tag & other) const {return mem == other.mem;}
		bool operator !=(const Tag & other) const {return mem != other.mem;}
		Tag & operator =(Tag const & other) {mem = other.mem; return *this;	}
		__INLINE DataType               GetDataType() const {assert(mem!=NULL); return mem->dtype;}
		__INLINE INMOST_MPI_Type        GetBulkDataType() const {assert(mem!=NULL); return mem->bulk_data_type;}
		__INLINE INMOST_DATA_ENUM_TYPE  GetBytesSize() const {assert(mem!=NULL); return mem->bytes_size;}
		__INLINE INMOST_DATA_ENUM_TYPE  GetSize() const {assert(mem!=NULL); return mem->size;}
		__INLINE std::string            GetTagName() const {assert(mem!=NULL); return mem->tagname;}
		__INLINE bool                   isDefined(ElementType type) const {assert(mem!=NULL && OneType(type)); return GetPosition(type) != ENUMUNDEF;}
		__INLINE bool                   isSparse(ElementType type) const {assert(mem!=NULL && OneType(type)); return mem->sparse[ElementNum(type)];}
		__INLINE bool                   isValid() const {return mem != NULL;}
		__INLINE Mesh *                 GetMeshLink() const {assert(mem!=NULL); return mem->m_link;}		
		__INLINE bool                   isSparseByDim(INMOST_DATA_INTEGER_TYPE typenum) const {assert(mem!=NULL); return mem->sparse[typenum];}
		__INLINE bool                   isDefinedByDim(INMOST_DATA_INTEGER_TYPE typenum) const {assert(mem!=NULL); return GetPositionByDim(typenum) != ENUMUNDEF;}
		__INLINE void                   SetBulkDataType(INMOST_MPI_Type type) {assert(mem!=NULL && mem->dtype == DATA_BULK ); mem->bulk_data_type = type;}
		friend class TagManager;
		friend class Storage;
		friend class Mesh;
	};
	
	class TagManager //implemented in tag.cpp
	{
	protected:
		TagManager();
		TagManager(const TagManager & other);
		TagManager & operator = (TagManager const & other);
		typedef chunk_array<INMOST_DATA_ENUM_TYPE,chunk_bits_empty>    empty_data;
		typedef chunk_array<Tag, chunk_bits_tags>                      tag_array_type;
		typedef chunk_bulk_array<chunk_bits_elems>                     dense_sub_type;
		typedef chunk_array<dense_sub_type,chunk_bits_dense_array>     dense_data_array_type;
		typedef struct{void * tag, * rec;}                             sparse_sub_record;
		typedef array< sparse_sub_record >                             sparse_sub_type;
		typedef chunk_array< sparse_sub_type,chunk_bits_elems>         sparse_data_array_type;
		typedef chunk_array<INMOST_DATA_INTEGER_TYPE,chunk_bits_elems> back_links_type;
		//typedef std::vector<INMOST_DATA_INTEGER_TYPE>                  back_links_type;
	public:
		typedef tag_array_type::iterator                               iteratorTag;
	public:
		virtual ~TagManager();
		bool                             HaveTag            (std::string name) const;
		Tag                              GetTag             (std::string name) const;
		void                             ListTagNames       (std::vector<std::string> & list) const;
		Tag                              CreateTag          (Mesh * m, std::string name, DataType dtype, ElementType etype, ElementType sparse, INMOST_DATA_ENUM_TYPE size = ENUMUNDEF); 
		virtual Tag                      DeleteTag          (Tag tag, ElementType mask); 
		bool                             ElementDefined     (Tag const & tag, ElementType etype) const;
	protected:
		void                             ReallocateData     (const Tag & t, INMOST_DATA_INTEGER_TYPE etypenum,INMOST_DATA_ENUM_TYPE new_size);
		void                             ReallocateData     (INMOST_DATA_INTEGER_TYPE etypenum, INMOST_DATA_ENUM_TYPE new_size);
		__INLINE sparse_sub_type const & GetSparseData      (int etypenum, int local_id) const {return sparse_data[etypenum][local_id];}
		__INLINE sparse_sub_type &       GetSparseData      (int etypenum, int local_id) {return sparse_data[etypenum][local_id];}
		__INLINE dense_sub_type const &  GetDenseData       (int pos) const {return dense_data[pos];}
		__INLINE dense_sub_type &        GetDenseData       (int pos) {return dense_data[pos];}
		static void                      CopyData           (const Tag & t, void * adata, const void * bdata);
		static void                      DestroyVariableData(const Tag & t, void * adata);
	protected:
		typedef tag_array_type::iterator            tag_iterator;
		typedef tag_array_type::const_iterator      tag_const_iterator;
	protected:
		tag_array_type                              tags;
		empty_data                                  empty_dense_data;
		dense_data_array_type                       dense_data;
		sparse_data_array_type                      sparse_data[6];
		back_links_type                             back_links[6];
	};
	
	/// Base class for Mesh, Element, and ElementSet classes.
	/// This base class is used for the Mesh class, as well as Element classes, and ElementSet class.
	/// 
	/// Storage class is used for representing different data objects in memory.
	/// Each data object is associated with corresponding Tag.
	class Storage //implemented in storage.cpp
	{
	public:
		/// Storage type for representing real values.
		typedef INMOST_DATA_REAL_TYPE     real; 
		/// Storage type for representing integer values.
		typedef INMOST_DATA_INTEGER_TYPE  integer;
		/// Storage type for representing one byte of abstact data.
		typedef INMOST_DATA_BULK_TYPE     bulk;
		/// type for representing unsigned integer values.
		typedef INMOST_DATA_ENUM_TYPE     enumerator;
		/// Storage type for representing references to Element.
		typedef HandleType                reference;
		/// Storage type for representing arrays of real values.
		typedef shell<real>               real_array;
		/// Storage type for representing arrays of integer values.
		typedef shell<integer>            integer_array;
		/// Storage type for representing abstact data as a series of bytes.
		typedef shell<bulk>               bulk_array;
		/// Storage type for representing arrays of Element references.
		class reference_array : public shell<HandleType>
		{
			Mesh * m;
		public:
			reference_array() : m(NULL) {}
			reference_array(Mesh * m, inner_reference_array & arr) : shell<reference>(arr), m(m) {} 
			reference_array(Mesh * m, reference * arr, size_type size) : shell<reference>(arr,size), m(m) {}
			~reference_array() {}
			void push_back(const Storage & elem);
			void push_back(HandleType h) {shell<reference>::push_back(h);} //is it needed?
			Element operator[] (size_type n);
			Element operator[] (size_type n) const;
			class iterator : public shell<HandleType>::iterator
			{
				Mesh * m;
			public:
				iterator(Mesh * m, const shell<HandleType>::iterator & other) : shell<HandleType>::iterator(other), m(m) {}
				iterator(const iterator & other) : shell<HandleType>::iterator(other), m(other.m) {}
				iterator & operator =(iterator const & other) {m = other.m; shell<HandleType>::iterator::operator=(other); return *this;}
				iterator & operator ++() {shell<HandleType>::iterator::operator++(); return *this;}
				iterator   operator ++(int) {iterator ret(*this); shell<HandleType>::iterator::operator++(); return ret;}
				iterator & operator --() {shell<HandleType>::iterator::operator--(); return *this;}
				iterator   operator --(int) {iterator ret(*this); shell<HandleType>::iterator::operator--(); return ret;}
				Element operator->();
			};
			class const_iterator : public shell<HandleType>::const_iterator
			{
				Mesh * m;
			public:
				const_iterator(Mesh * m, const shell<HandleType>::const_iterator & other) : shell<HandleType>::const_iterator(other) , m(m) {}
				const_iterator(const const_iterator & other) : shell<HandleType>::const_iterator(other), m(other.m) {}
				const_iterator & operator =(const_iterator const & other) {m = other.m; shell<HandleType>::const_iterator::operator=(other); return *this;}
				const_iterator & operator ++() {shell<HandleType>::const_iterator::operator++(); return *this;}
				const_iterator   operator ++(int) {const_iterator ret(*this); shell<HandleType>::const_iterator::operator++(); return ret;}
				const_iterator & operator --() {shell<HandleType>::const_iterator::operator--(); return *this;}
				const_iterator   operator --(int) {const_iterator ret(*this); shell<HandleType>::const_iterator::operator--(); return ret;}
				Element operator->();
			};
			class reverse_iterator : public shell<HandleType>::reverse_iterator
			{
				Mesh * m;
			public:
				reverse_iterator(Mesh * m, const shell<HandleType>::reverse_iterator & other) : shell<HandleType>::reverse_iterator(other), m(m) {}
				reverse_iterator(const reverse_iterator & other) : shell<HandleType>::reverse_iterator(other), m(other.m)  {}
				reverse_iterator & operator =(reverse_iterator const & other) {m = other.m; shell<HandleType>::reverse_iterator::operator=(other); return *this;}
				reverse_iterator & operator ++() {shell<HandleType>::reverse_iterator::operator++(); return *this;}
				reverse_iterator   operator ++(int) {reverse_iterator ret(*this); shell<HandleType>::reverse_iterator::operator++(); return ret;}
				reverse_iterator & operator --() {shell<HandleType>::reverse_iterator::operator--(); return *this;}
				reverse_iterator   operator --(int) {reverse_iterator ret(*this); shell<HandleType>::reverse_iterator::operator--(); return ret;}
				Element operator->();
			};
			class const_reverse_iterator : public shell<HandleType>::const_reverse_iterator
			{
				Mesh * m;
			public:
				const_reverse_iterator(Mesh * m, const shell<HandleType>::const_reverse_iterator & other) : shell<HandleType>::const_reverse_iterator(other), m(m) {}
				const_reverse_iterator(const const_reverse_iterator & other) : shell<HandleType>::const_reverse_iterator(other), m(other.m) {}
				const_reverse_iterator & operator =(const_reverse_iterator const & other) {m = other.m; shell<HandleType>::const_reverse_iterator::operator=(other); return *this;}
				const_reverse_iterator & operator ++() {shell<HandleType>::const_reverse_iterator::operator++(); return *this;}
				const_reverse_iterator   operator ++(int) {const_reverse_iterator ret(*this); shell<HandleType>::const_reverse_iterator::operator++(); return ret;}
				const_reverse_iterator & operator --() {shell<HandleType>::const_reverse_iterator::operator--(); return *this;}
				const_reverse_iterator   operator --(int) {const_reverse_iterator ret(*this); shell<HandleType>::const_reverse_iterator::operator--(); return ret;}
				Element operator->();
			};
			iterator begin() {return iterator(m,shell<HandleType>::begin());}
			const_iterator begin() const {return const_iterator(m,shell<HandleType>::begin());}
			reverse_iterator rbegin() {return reverse_iterator(m,shell<HandleType>::rbegin());}
			const_reverse_iterator rbegin() const {return const_reverse_iterator(m,shell<HandleType>::rbegin());}
			iterator end() {return iterator(m,shell<HandleType>::end());}
			const_iterator end() const {return const_iterator(m,shell<HandleType>::end());}
			reverse_iterator rend() {return reverse_iterator(m,shell<HandleType>::rend());}
			const_reverse_iterator rend() const {return const_reverse_iterator(m,shell<HandleType>::rend());}
			void SetMeshLink(Mesh * new_m) {m = new_m;}
		};
		//typedef shell<reference>          reference_array;
	protected:
		HandleType                        handle;
		HandleType *                      handle_link;
	private:
		Mesh *                            m_link;
	public:
		Storage(const Storage & other) : handle(other.handle), handle_link(other.handle_link), m_link(other.m_link) {}
		Storage(Mesh * mesh, HandleType handle) : handle(handle), handle_link(NULL), m_link(mesh) {}
		/// This constructor allows for remote handle modification
		Storage(Mesh * mesh, HandleType * handle) : handle(*handle), handle_link(handle), m_link(mesh) {}
		/// If there is a link to handle provided (automatically by ElementArray and reference_array),
		/// then remote handle value will be modified
		Storage & operator =(Storage const & other); 
		bool      operator <(const Storage & other) const {return handle < other.handle;}
		bool      operator >(const Storage & other) const {return handle > other.handle;}
		bool      operator <=(const Storage & other) const {return handle <= other.handle;}
		bool      operator >=(const Storage & other) const {return handle >= other.handle;}
		bool      operator ==(const Storage & other) const {return handle == other.handle;} //m_link may be NULL if Invalidxxx is used
		bool      operator !=(const Storage & other) const {return handle != other.handle;} //m_link may be NULL if Invalidxxx is used
		Storage * operator ->() {return this;}
		const Storage * operator->() const {return this;}
		Storage & self() {return *this;}
		const Storage & self() const {return *this;}
		virtual ~Storage() {}
	public:
		/// Retrieve real value associated with Tag.
		real      &              Real            (const Tag & tag) const;
		/// Retrieve integer value associated with Tag.
		integer   &              Integer         (const Tag & tag) const;
		/// Retrieve one byte of abstract data associated with Tag.
		bulk      &              Bulk            (const Tag & tag) const;
		/// Retrieve Element reference associated with Tag.
		reference &              Reference       (const Tag & tag) const;
		/// Retrieve array of real values associated with Tag.
		real_array               RealArray       (const Tag & tag) const;
		/// Retrieve array of integer values associated with Tag.
		integer_array            IntegerArray    (const Tag & tag) const;
		/// Retrieve abstract data associated with Tag as a series of bytes.
		bulk_array               BulkArray       (const Tag & tag) const;
		/// Retrieve array of Element references associated with Tag.
		reference_array          ReferenceArray  (const Tag & tag) const;
		
		//optimized data requests for dense data with fixed size
		real_array               RealArrayDF     (const Tag & tag) const;
		integer_array            IntegerArrayDF  (const Tag & tag) const;
		bulk_array               BulkArrayDF     (const Tag & tag) const;
		reference_array          ReferenceArrayDF(const Tag & tag) const;
		real      &              RealDF          (const Tag & tag) const;
		integer   &              IntegerDF       (const Tag & tag) const;
		bulk      &              BulkDF          (const Tag & tag) const;
		reference &              ReferenceDF     (const Tag & tag) const;
		
		//optimized data requests for dense data with variable size
		real_array               RealArrayDV     (const Tag & tag) const;
		integer_array            IntegerArrayDV  (const Tag & tag) const;
		bulk_array               BulkArrayDV     (const Tag & tag) const;
		reference_array          ReferenceArrayDV(const Tag & tag) const;
		real      &              RealDV          (const Tag & tag) const;
		integer   &              IntegerDV       (const Tag & tag) const;
		bulk      &              BulkDV          (const Tag & tag) const;
		reference &              ReferenceDV     (const Tag & tag) const;
		
		/// Return the data length associated with Tag.
		/// For abstract data return the number of bytes, otherwise return the length of associated array. 
		/// @see Storage::SetDataSize
		INMOST_DATA_ENUM_TYPE    GetDataSize     (const Tag & tag) const; //For DATA_BULK return number of bytes, otherwise return the length of array
		/// Set the length of  data associated with Tag.
		/// @param tag Identifying Tag.
		/// @param new_size The number of bytes for abstract data, otherwise the length of the array.
		/// @see Storage::GetDataSize
		void                     SetDataSize     (const Tag & tag,INMOST_DATA_ENUM_TYPE new_size) const;
		/// Extract part of the data associated with Tag.
		/// Copy part of the associated array or data to the destination memory.
		/// @param tag Identifying Tag.
		/// @param shift Starting position of the copied data.
		/// For abstact data – number of bytes to skip, otherwise number of values to skip.
		/// @param size Number of elements to copy.
		/// For abstact data – number of bytes to copy, otherwise number of values to copy.
		/// @param data Destination position to copy data to.
		/// @see Storage::SetData
		void                     GetData         (const Tag & tag, INMOST_DATA_ENUM_TYPE shift, INMOST_DATA_ENUM_TYPE size, void * data) const;
		void                     SetData         (const Tag & tag, INMOST_DATA_ENUM_TYPE shift, INMOST_DATA_ENUM_TYPE size, const void * data) const;
		void                     DelData         (const Tag & tag) const;
		/// Deallocates space allocated for sparse data, frees variable array if necessary
		void                     DelSparseData   (const Tag & tag) const;
		/// Frees variable array or fills field with zeroes
		void                     DelDenseData    (const Tag & tag) const;
		/// Check if any data is associated with Tag.
		bool                     HaveData        (const Tag & tag) const;
		__INLINE ElementType     GetElementType  () const {return GetHandleElementType(handle);}
		__INLINE integer         GetElementNum   () const {return GetHandleElementNum(handle);}
		void                     SetMarker       (MarkerType n) const;
		bool                     GetMarker       (MarkerType n) const;
		void                     RemMarker       (MarkerType n) const;
		void                     ClearMarkerSpace() const;
		void                     GetMarkerSpace  (bulk copy[MarkerFields]) const;
		void                     SetMarkerSpace  (bulk source[MarkerFields]) const;
		__INLINE integer         LocalID         () const {return GetHandleID(handle);}
		/// This number is guaranteed to be between 0 and Mesh::NumberOf(type of element)
		/// after Mesh::ReorderEmpty
		integer                  DataLocalID     () const;
		bool                     isValid         () const;
		__INLINE Mesh *          GetMeshLink     () const {return m_link;}
		__INLINE HandleType      GetHandle       () const {return handle;}
		friend class Mesh;
	};
	
	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);}
			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);}
			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;
	};
	
			
	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
		Node                        getAsNode               () const;
		Edge                        getAsEdge               () const;
		Face                        getAsFace               () const;
		Cell                        getAsCell               () const;
		ElementSet                  getAsSet                () const;
		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
		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.