eset.cpp

//#define _HAS_ITERATOR_DEBUGGING 0

#include "inmost.h"
#if defined(USE_MESH)

namespace INMOST
{
	//shortcuts
	

	ElementSet::iterator & ElementSet::iterator::operator ++() 
	{
		++pos;
		Element::adj_type::size_type end = ptr->size();
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			while(pos < end && (ptr->at(pos) == InvalidHandle() || m->GetMarker(ptr->at(pos),hm)) ) ++pos; 
		}
		else 
		{
			while(pos < end && ptr->at(pos) == InvalidHandle()) ++pos; 
		}
		return *this;
	}


	Storage::enumerator ElementSet::nbSorted() const
	{
		return hSorted(GetMeshLink()->HighConn(GetHandle()));
	}
	ElementSet ElementSet::GetParent() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return ElementSet(GetMeshLink(),hParent(hc));
	}

	ElementSet ElementSet::GetSibling() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return ElementSet(GetMeshLink(),hSibling(hc));
	}

	ElementSet ElementSet::GetChild() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return ElementSet(GetMeshLink(),hChild(hc));
	}


	bool ElementSet::HaveParent() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return hParent(hc) != InvalidHandle();
	}

	bool ElementSet::HaveSibling() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return hSibling(hc) != InvalidHandle();
	}

	bool ElementSet::HaveChild() const
	{
		assert(GetElementType() == ESET);
		Element::adj_type & hc = GetMeshLink()->HighConn(GetHandle());
		return hChild(hc) != InvalidHandle();
	}

	void ElementSet::AddChild(const ElementSet & child) const
	{
		assert(GetElementType() == ESET);
		assert(child->GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & parent_hc = m->HighConn(GetHandle());
		assert(isValid());
		if(hChild(parent_hc) == InvalidHandle()) 
		{
			Element::adj_type & child_hc = m->HighConn(child->GetHandle());
			assert(hParent(child_hc) == InvalidHandle()); // child is not connected to another parent
			hParent(child_hc) = GetHandle(); //connect child to me
			hChild(parent_hc) = child->GetHandle(); //connect me to child
		}
		else 
			ElementSet(m,hChild(parent_hc))->AddSibling(child); //connect as sibling to my children
	}

	void ElementSet::AddSibling(const ElementSet & sibling) const
	{
		assert(GetElementType() == ESET);
		assert(sibling->GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		//perform checks
		Element::adj_type & current_hc = m->HighConn(GetHandle());
		assert(isValid());
		assert( hParent(current_hc) != InvalidHandle() ); //no parent for current set
		if( hSibling(current_hc) == InvalidHandle() ) //no sibling here
		{
			Element::adj_type & sibling_hc = m->HighConn(sibling->GetHandle());
			assert( hParent(sibling_hc) == InvalidHandle() );
			hSibling(current_hc) = sibling->GetHandle(); //connect sibling to me
			hParent(sibling_hc) = hParent(current_hc); //connect sibling to my parent
		}
		else 
			ElementSet(m,hSibling(current_hc))->AddSibling(sibling);
	}

	void ElementSet::RemChild(const ElementSet & child) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & parent_hc = m->HighConn(GetHandle());
		assert(isValid());
		assert(hChild(parent_hc) != InvalidHandle());
		if( hChild(parent_hc) == child->GetHandle() )
		{
			Element::adj_type & child_hc = m->HighConn(child->GetHandle());
			assert(hParent(child_hc) == GetHandle()); //was connected to me
			hParent(child_hc) = InvalidHandle();
			hChild(parent_hc) = hSibling(child_hc); //get next sibling as my child
		}
		else ElementSet(m,hChild(parent_hc))->RemSibling(child);
	}

	void ElementSet::RemSibling(const ElementSet & sibling) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & current_hc = m->HighConn(GetHandle());
		assert(isValid());
		assert(hSibling(current_hc) != InvalidHandle());
		if( hSibling(current_hc) == sibling->GetHandle())
		{
			Element::adj_type & sibling_hc = m->HighConn(sibling->GetHandle());
			assert(hParent(sibling_hc) == hParent(current_hc));
			hParent(sibling_hc) = InvalidHandle(); //remove it's parent
			hSibling(current_hc) = hSibling(sibling_hc); //get next sibling as my sibling
		}
		else ElementSet(m,hSibling(current_hc))->RemSibling(sibling);
	}

	Storage::enumerator ElementSet::CountChildren() const
	{
		assert(GetElementType() == ESET);
		ElementSet child = GetChild();
		return child->isValid() ? child->CountSiblings() : 0;
	}

	Storage::enumerator ElementSet::CountSiblings() const
	{
		assert(GetElementType() == ESET);
		ElementSet next = GetSibling();
		return next->isValid() ? next->CountSiblings()+1 : 1;
	}

	std::string ElementSet::GetName() const
	{
		assert(GetElementType() == ESET);
		std::string ret;
		bulk_array arr = BulkArrayDV(GetMeshLink()->SetNameTag());
		ret.resize(arr.size());
		for(bulk_array::size_type it = 0; it < arr.size(); ++it)
			ret[it] = arr[it];
		return ret;
	}

	HandleType * ElementSet::getHandles() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		return lc.data();
	}

	Storage::enumerator ElementSet::nbHandles() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		return static_cast<enumerator>(lc.size());
	}

	void ElementSet::SetMarkerElements(MarkerType mrk, ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( (GetHandleElementType(*it) & etype) && !m->GetMarker(*it,hm) ) //check for type filters invalid handles
					m->SetMarker(*it,mrk);
		}
		else
		{
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(*it) & etype )
					m->SetMarker(*it,mrk);
		}
	}


  void ElementSet::SetPrivateMarkerElements(MarkerType mrk, ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( (GetHandleElementType(*it) & etype) && !m->GetMarker(*it,hm) ) //check for type filters invalid handles
					m->SetPrivateMarker(*it,mrk);
		}
		else
		{
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(*it) & etype )
					m->SetPrivateMarker(*it,mrk);
		}
	}
	void ElementSet::RemMarkerElements(MarkerType mrk, ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( (GetHandleElementType(*it) & etype) && !m->GetMarker(*it,hm) )//check for type filters invalid handles
					m->RemMarker(*it,mrk);
		}
		else
		{
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(*it) & etype )
					m->RemMarker(*it,mrk);
		}
	}
  void ElementSet::RemPrivateMarkerElements(MarkerType mrk, ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( (GetHandleElementType(*it) & etype) && !m->GetMarker(*it,hm) )//check for type filters invalid handles
					m->RemPrivateMarker(*it,mrk);
		}
		else
		{
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(*it) & etype )
					m->RemPrivateMarker(*it,mrk);
		}
	}

	
	ElementArray<Element> ElementSet::getAdjElements(ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & etype )
					ret.push_back(lc[it]);
		}
		return ret;
	}

	Storage::enumerator ElementSet::nbAdjElements(ElementType etype) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type const & lc = m->LowConn(GetHandle());
		enumerator ret = 0;
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					++ret;
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & etype )
					++ret;
		}
		return ret;
	}

	ElementArray<Element> ElementSet::getAdjElements(ElementType etype, MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		return ret;
	}

	Storage::enumerator ElementSet::nbAdjElements(ElementType etype, MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		enumerator ret = 0;
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					++ret;
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & etype) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					++ret;
		}
		return ret;
	}

	ElementArray<Node> ElementSet::getNodes() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Node> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & NODE) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & NODE )
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Node> ElementSet::getNodes(MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Node> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & NODE) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & NODE) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Edge> ElementSet::getEdges() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Edge> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & EDGE) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & EDGE )
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Edge> ElementSet::getEdges(MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Edge> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & EDGE) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & EDGE) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Face> ElementSet::getFaces() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Face> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & FACE) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & FACE )
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Face> ElementSet::getFaces(MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Face> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & FACE) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & FACE) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Cell> ElementSet::getCells() const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Cell> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & CELL) && !m->GetMarker(lc[it],hm) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)//check for type filters invalid handles
				if( GetHandleElementType(lc[it]) & CELL )
					ret.push_back(lc[it]);
		}
		return ret;
	}

	ElementArray<Cell> ElementSet::getCells(MarkerType select, bool invert) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Cell> ret(m);
		Element::adj_type const & lc = m->LowConn(GetHandle());
		if( m->isMeshModified() )
		{
			MarkerType hm = m->HideMarker();
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & CELL) && !m->GetMarker(lc[it],hm) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		else
		{
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( (GetHandleElementType(lc[it]) & CELL) && (invert ^ m->GetMarker(lc[it],select)) )//check for type filters invalid handles
					ret.push_back(lc[it]);
		}
		return ret;
	}

	void ElementSet::PutElement(HandleType e) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		switch(BulkDF(m->SetComparatorTag()))
		{
		case UNSORTED_COMPARATOR:
			{
				Element::adj_type & hc = m->HighConn(GetHandle());
				if( hc.size() > high_conn_reserved )
				{
					lc[static_cast<enumerator>(hc.back())] = e;
					hc.pop_back();
				}
				else lc.push_back(e);
			}
			break;
		case   HANDLE_COMPARATOR:
		case  HIERARCHY_COMPARATOR:
		case GLOBALID_COMPARATOR:
		case CENTROID_COMPARATOR:
			lc.push_back(e);
			break;
		}
	}

	void ElementSet::PutElements(const HandleType * h, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		switch(BulkDF(m->SetComparatorTag()))
		{
		case UNSORTED_COMPARATOR:
			{
				const HandleType * h2 = h;
				Element::adj_type & hc = m->HighConn(GetHandle());
				while( hc.size() > high_conn_reserved )
				{
					lc[static_cast<enumerator>(hc.back())] = *h2;
					hc.pop_back();
					h2++;
				}
				lc.insert(lc.end(),h2, h+num);
			}
			break;
		case   HANDLE_COMPARATOR:
		case HIERARCHY_COMPARATOR:
		case GLOBALID_COMPARATOR:
		case CENTROID_COMPARATOR:
			lc.insert(lc.end(),h,h+num);
			break;
		}
	}

	void ElementSet::AddElement(HandleType e) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->HighConn(GetHandle());
		Element::adj_type::iterator find;
		Element::adj_type::size_type sorted = static_cast<Element::adj_type::size_type>(hSorted(hc));
		switch(BulkDF(m->SetComparatorTag()))
		{
		case UNSORTED_COMPARATOR:
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( *it == e ) return;
			if( hc.size() > high_conn_reserved )
			{
				lc[static_cast<enumerator>(hc.back())] = e;
				hc.pop_back();
			}
			else lc.push_back(e);
			return;
		case   HANDLE_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e);
			break;
		case  HIERARCHY_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e,Mesh::HierarchyComparator(m));
			break;
		case GLOBALID_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e,Mesh::GlobalIDComparator(m));
			break;
		case CENTROID_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e,Mesh::CentroidComparator(m));
			break;
		}
		if( find != lc.begin()+sorted )
		{
			if( *find == e )
				return;
		}
		lc.insert(find,1,e);
		hSorted(hc)++;
	}

	void ElementSet::AddElements(const HandleType * h, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->HighConn(GetHandle());
		switch(BulkDF(m->SetComparatorTag()))
		{
		case UNSORTED_COMPARATOR:
			{
				Element::adj_type & hc = m->HighConn(GetHandle());
				MarkerType mrk = m->CreatePrivateMarker();
				SetPrivateMarkerElements(mrk);
				enumerator it = 0;
				while( it != num )
				{
					if( !m->GetPrivateMarker(h[it],mrk) )
					{
						if( hc.size() > high_conn_reserved )
						{
							lc[static_cast<enumerator>(hc.back())] = h[it];
							hc.pop_back();
						}
						else lc.push_back(hc[it]);
					}
					++it;
				}
				RemPrivateMarkerElements(mrk);
				m->ReleasePrivateMarker(mrk);
			}
			break;
		case HANDLE_COMPARATOR:
			{
				// sort 
				array<HandleType> sorted(h,h+num);
				m->SortHandles(sorted.data(),sorted.size());
				//now merge with sorted part
				Element::adj_type merge(lc.size()+sorted.size());
				Element::adj_type::iterator it = lc.begin(), iend = lc.begin()+hSorted(hc);
				Element::adj_type::iterator jt = sorted.begin(), jend = sorted.end();
				Element::adj_type::iterator mit = merge.begin(), mitprev = merge.begin();
				//insert first element
				while(it != iend && jt != jend)
				{
					if( *it == InvalidHandle() ) ++it;
					else if( *jt == InvalidHandle() ) ++jt;
					else if( *it < *jt ) {*mit++ = *it++; break;}
					else {*mit++ = *jt++; break;}
				}
				//now insertion checks for duplicates
				while(it != iend && jt != jend)
				{
					if( *it == InvalidHandle() || *it == *mitprev ) ++it;
					else if( *jt == InvalidHandle() || *jt == *mitprev ) ++jt;
					else if( *it < *jt ) {mitprev = mit; *mit = *it; ++mit; ++it;}
					else {mitprev = mit; *mit = *jt; ++mit; ++jt;}
				}
				//attach rest
				while(it != iend && *it != InvalidHandle() && *it != *mitprev )
				{
					mitprev = mit;
					*mit = *it;
					++mit;
					++it;
				}
				while(jt != jend && *jt != InvalidHandle() && *jt != *mitprev )
				{
					mitprev = mit;
					*mit = *jt;
					++mit;
					++jt;
				}
				//indicate end of sorted part
				hSorted(hc) = static_cast<HandleType>(mit-merge.begin());
				//attach unsorted part with control for duplicates
				MarkerType mrk = m->CreatePrivateMarker();
				m->SetPrivateMarkerArray(merge.data(),hSorted(hc),mrk);
				while(it != lc.end() && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					m->SetPrivateMarker(*it,mrk);
					*mit++ = *it++;
				}
				merge.resize(static_cast<Element::adj_type::size_type>(mit-merge.begin()));
				m->RemPrivateMarkerArray(merge.data(),static_cast<enumerator>(merge.size()),mrk);
				m->ReleasePrivateMarker(mrk);
				// set merged elemetns as current
				lc.swap(merge);
				//indicate that there are no empty spaces
				hc.resize(high_conn_reserved);
			}
			break;
		case GLOBALID_COMPARATOR:
			{
				// sort 
				array<HandleType> sorted(h,h+num);
				m->SortByGlobalID(sorted.data(),sorted.size());
				//now merge with sorted part
				Element::adj_type merge(lc.size()+sorted.size());
				Element::adj_type::iterator it = lc.begin(), iend = lc.begin()+hSorted(hc);
				Element::adj_type::iterator jt = sorted.begin(), jend = sorted.end();
				Element::adj_type::iterator mit = merge.begin();
				MarkerType mrk = m->CreatePrivateMarker();
        //now insertion checks for duplicates
				while(it != iend && jt != jend)
				{
					if( *it == InvalidHandle() || m->GetPrivateMarker(*it,mrk) ) ++it;
					else if( *jt == InvalidHandle() || m->GetPrivateMarker(*jt,mrk) ) ++jt;
					else if( Mesh::GlobalIDComparator(m)(*it,*jt) ) 
          {
            *mit = *it; 
            m->SetPrivateMarker(*mit,mrk);
            ++mit; 
            ++it;
          }
					else 
          {
            *mit = *jt;
            m->SetPrivateMarker(*mit,mrk);
            ++mit; 
            ++jt;
          }
				}
				//attach rest
				while(it != iend && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					*mit = *it;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++it;
				}
				while(jt != jend && *jt != InvalidHandle() && !m->GetPrivateMarker(*jt,mrk) )
				{
					*mit = *jt;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++jt;
				}
				//indicate end of sorted part
				hSorted(hc) = static_cast<HandleType>(mit-merge.begin());
				//attach unsorted part
				while(it != lc.end() && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					m->SetPrivateMarker(*it,mrk);
					*mit++ = *it++;
				}
				merge.resize(static_cast<Element::adj_type::size_type>(mit-merge.begin()));
				m->RemPrivateMarkerArray(merge.data(),static_cast<enumerator>(merge.size()),mrk);
				m->ReleasePrivateMarker(mrk);
				// set merged elemetns as current
				lc.swap(merge);
				//indicate that there are no empty spaces
				hc.resize(high_conn_reserved);
			}
			break;
		case HIERARCHY_COMPARATOR:
			{
				// sort 
				array<HandleType> sorted(h,h+num);
				std::sort(sorted.begin(),sorted.end(),Mesh::HierarchyComparator(m));
				//now merge with sorted part
				Element::adj_type merge(lc.size()+sorted.size());
				Element::adj_type::iterator it = lc.begin(), iend = lc.begin()+hSorted(hc);
				Element::adj_type::iterator jt = sorted.begin(), jend = sorted.end();
				Element::adj_type::iterator mit = merge.begin(), mitprev = merge.begin();
				MarkerType mrk = m->CreatePrivateMarker();
				//now insertion checks for duplicates
				while(it != iend && jt != jend)
				{
					if( *it == InvalidHandle() || m->GetPrivateMarker(*it,mrk) ) ++it;
					else if( *jt == InvalidHandle() || m->GetPrivateMarker(*jt,mrk) ) ++jt;
					else if( Mesh::HierarchyComparator(m)(*it,*jt) ) 
          {
            *mit = *it; 
            m->SetPrivateMarker(*mit,mrk);
            ++mit; ++it;
          }
					else 
          {
            *mit = *jt; 
            m->SetPrivateMarker(*mit,mrk);
            ++mit; ++jt;
          }
				}
				//attach rest
				while(it != iend && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					*mit = *it;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++it;
				}
				while(jt != jend && *jt != InvalidHandle() && !m->GetPrivateMarker(*jt,mrk) )
				{
					*mit = *jt;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++jt;
				}
				//indicate end of sorted part
				hSorted(hc) = static_cast<HandleType>(mit-merge.begin());
				//attach unsorted part
				while(it != lc.end() && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					m->SetPrivateMarker(*it,mrk);
					*mit++ = *it++;
				}
				merge.resize(static_cast<Element::adj_type::size_type>(mit-merge.begin()));
				m->RemPrivateMarkerArray(merge.data(),static_cast<enumerator>(merge.size()),mrk);
				m->ReleasePrivateMarker(mrk);
				// set merged elemetns as current
				lc.swap(merge);
				//indicate that there are no empty spaces
				hc.resize(high_conn_reserved);
			}
			break;
		case CENTROID_COMPARATOR:
			{
				// sort 
				array<HandleType> sorted(h,h+num);
				std::sort(sorted.begin(),sorted.end(),Mesh::CentroidComparator(m));
				//now merge with sorted part
				Element::adj_type merge(lc.size()+sorted.size());
				Element::adj_type::iterator it = lc.begin(), iend = lc.begin()+hSorted(hc);
				Element::adj_type::iterator jt = sorted.begin(), jend = sorted.end();
				Element::adj_type::iterator mit = merge.begin(), mitprev = merge.begin();
				MarkerType mrk = m->CreatePrivateMarker();
				//now insertion checks for duplicates
				while(it != iend && jt != jend)
				{
					if( *it == InvalidHandle() || m->GetPrivateMarker(*it,mrk) ) ++it;
					else if( *jt == InvalidHandle() || m->GetPrivateMarker(*jt,mrk) ) ++jt;
					else if( Mesh::CentroidComparator(m)(*it,*jt) ) 
          {
            *mit = *it; 
            m->SetPrivateMarker(*mit,mrk);
            ++mit; ++it;
          }
					else 
          {
            *mit = *jt; 
            m->SetPrivateMarker(*mit,mrk);
            ++mit; ++it;
          }
				}
				//attach rest
				while(it != iend && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					*mit = *it;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++it;
				}
				while(jt != jend && *jt != InvalidHandle() && !m->GetPrivateMarker(*jt,mrk) )
				{
					*mit = *jt;
          m->SetPrivateMarker(*mit,mrk);
					++mit;
					++jt;
				}
				//indicate end of sorted part
				hSorted(hc) = static_cast<HandleType>(mit-merge.begin());
				//attach unsorted part
				while(it != lc.end() && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
				{
					m->SetPrivateMarker(*it,mrk);
					*mit++ = *it++;
				}
				merge.resize(static_cast<Element::adj_type::size_type>(mit-merge.begin()));
				m->RemPrivateMarkerArray(merge.data(),static_cast<enumerator>(merge.size()),mrk);
				m->ReleasePrivateMarker(mrk);
				// set merged elemetns as current
				lc.swap(merge);
				//indicate that there are no empty spaces
				hc.resize(high_conn_reserved);
			}
			break;
		}
	}

	void ElementSet::RemoveElement(const Storage & e) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->HighConn(GetHandle());
		Element::adj_type::iterator find;
		Element::adj_type::size_type i, sorted = static_cast<Element::adj_type::size_type>(hSorted(hc));
		switch(BulkDF(m->SetComparatorTag()))
		{
		case UNSORTED_COMPARATOR:
			for(Element::adj_type::size_type it = 0; it < lc.size(); ++it)
				if( lc[it] == e->GetHandle() )
				{
					lc[it] = InvalidHandle();
					hc.push_back(static_cast<HandleType>(it));
					break;
				}
			return;
		case   HANDLE_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e->GetHandle());
			break;
		case  HIERARCHY_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e->GetHandle(),Mesh::HierarchyComparator(m));
			break;
		case GLOBALID_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e->GetHandle(),Mesh::GlobalIDComparator(m));
			break;
		case CENTROID_COMPARATOR:
			find = std::lower_bound(lc.begin(),lc.begin()+sorted,e->GetHandle(),Mesh::CentroidComparator(m));
			break;
		}
		if( find != lc.begin()+sorted )
		{
			if( *find == e->GetHandle() )
			{
				*find = InvalidHandle();
				hc.push_back(static_cast<HandleType>(find-lc.begin()));
			}
		}
		else for(i = sorted; i < lc.size(); ++i) //unsorted part
			if( lc[i] == e->GetHandle() )
			{
				lc[i] = InvalidHandle();
				hc.push_back(static_cast<HandleType>(i));
			}
	}

	void ElementSet::RemoveElements(const HandleType * h, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->HighConn(GetHandle());
		Element::adj_type::size_type i;
		MarkerType mrk = m->CreatePrivateMarker();
		m->SetPrivateMarkerArray(h,num,mrk);
		for(i = 0; i < lc.size(); i++)
		{
			if( m->GetPrivateMarker(lc[i],mrk) )
			{
				lc[i] = InvalidHandle();
				hc.push_back(i);
			}
		}
		m->RemPrivateMarkerArray(h,num,mrk);
		m->ReleasePrivateMarker(mrk);
	}

	ElementArray<Element> ElementSet::Union(const ElementSet & other) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && !m->GetPrivateMarker(lc[i],mrk))
			{
				ret.push_back(lc[i]);
				m->SetPrivateMarker(ret.atback(),mrk);
			}
		}
		Element::adj_type & other_lc = m->LowConn(other->GetHandle());
		for(const HandleType * h = other_lc.data(); h != other_lc.data() + other_lc.size(); ++h)
		{
			if( *h != InvalidHandle() && !m->GetPrivateMarker(*h,mrk) )
			{
				ret.push_back(*h);
				m->SetPrivateMarker(ret.atback(),mrk);
			}
		}
		ret.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}

	ElementArray<Element> ElementSet::Union(const HandleType * handles, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && !m->GetPrivateMarker(lc[i],mrk))
			{
				ret.push_back(lc[i]);
				m->SetPrivateMarker(ret.atback(),mrk);
			}
		}
		for(const HandleType * h = handles; h != handles+num; ++h)
		{
			if( *h != InvalidHandle() && !m->GetPrivateMarker(*h,mrk) )
			{
				ret.push_back(*h);
				m->SetPrivateMarker(ret.atback(),mrk);
			}
		}
		ret.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}

	ElementArray<Element> ElementSet::Difference(const ElementSet & other) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		other->SetPrivateMarkerElements(mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && !m->GetPrivateMarker(lc[i],mrk))
			{
				ret.push_back(lc[i]);
				m->SetMarker(ret.atback(),mrk);
			}
		}
		other->RemPrivateMarkerElements(mrk);
		ret.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}

	ElementArray<Element> ElementSet::Difference(const HandleType * handles, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		m->SetPrivateMarkerArray(handles,num,mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && !m->GetPrivateMarker(lc[i],mrk))
			{
				ret.push_back(lc[i]);
				m->SetPrivateMarker(ret.atback(),mrk);
			}
		}
		m->RemPrivateMarkerArray(handles,num,mrk);
		ret.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}

	ElementArray<Element> ElementSet::Intersection(const ElementSet & other) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker(), mrk2 = m->CreatePrivateMarker();
		other->SetPrivateMarkerElements(mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && m->GetPrivateMarker(lc[i],mrk) && !m->GetPrivateMarker(lc[i],mrk2))
			{
				ret.push_back(lc[i]);
				m->SetPrivateMarker(lc[i],mrk2);
			}
		}
		other->RemPrivateMarkerElements(mrk);
		ret.RemPrivateMarker(mrk2);
		m->ReleasePrivateMarker(mrk2);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}

	ElementArray<Element> ElementSet::Intersection(const HandleType * handles, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		ElementArray<Element> ret(m);
		Element::adj_type & lc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker(), mrk2 = m->CreatePrivateMarker();
		m->SetPrivateMarkerArray(handles,num,mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && m->GetPrivateMarker(lc[i],mrk) && !m->GetPrivateMarker(lc[i],mrk2))
			{
				ret.push_back(lc[i]);
				m->SetPrivateMarker(lc[i],mrk2);
			}
		}
		m->RemPrivateMarkerArray(handles,num,mrk);
		ret.RemPrivateMarker(mrk2);
		m->ReleasePrivateMarker(mrk2);
		m->ReleasePrivateMarker(mrk);
		return ret;
	}


	void ElementSet::Unite(const ElementSet & other) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		bulk cmp = BulkDF(m->SetComparatorTag());
		Element::adj_type & other_lc = m->LowConn(other->GetHandle());
		if( cmp == UNSORTED_COMPARATOR || other->BulkDF(m->SetComparatorTag()) != cmp )
			Unite(other_lc.data(),static_cast<enumerator>(other_lc.size()));
		else //other set sorted by the same comparator, just merge them
		{
			Element::adj_type & other_hc = m->HighConn(other->GetHandle());
			Element::adj_type & hc = m->HighConn(GetHandle());
			Element::adj_type & lc = m->LowConn(GetHandle());
			Element::adj_type merge(lc.size()+other_lc.size());
			Element::adj_type::iterator it = lc.begin(), iend = lc.begin()+hSorted(hc);
			Element::adj_type::iterator jt = other_lc.begin(), jend = other_lc.end()+hSorted(other_hc);
			Element::adj_type::iterator mit = merge.begin(), mitprev = merge.begin();
			MarkerType mrk = m->CreatePrivateMarker();
			//now insertion checks for duplicates
			while(it != iend && jt != jend)
			{
				if( *it == InvalidHandle() || m->GetPrivateMarker(*it,mrk) ) ++it;
				else if( *jt == InvalidHandle() || m->GetPrivateMarker(*jt,mrk) ) ++jt;
				else 
				{
					bool cmp_res = false;
					switch(cmp)
					{
					case HANDLE_COMPARATOR: cmp_res = *it < *jt; break;
					case HIERARCHY_COMPARATOR: cmp_res = Mesh::HierarchyComparator(m)(*it,*jt); break;
					case CENTROID_COMPARATOR: cmp_res = Mesh::CentroidComparator(m)(*it,*jt); break;
					case GLOBALID_COMPARATOR: cmp_res = Mesh::GlobalIDComparator(m)(*it,*jt); break;
					default: assert(false);
					}
					if( cmp_res ) 
					{
						*mit = *it;
            m->SetPrivateMarker(*mit,mrk);
						++mit; 
            ++it;
					}
					else 
					{
						*mit = *jt; 
            m->SetPrivateMarker(*mit,mrk);
						++mit; 
						++jt;
					}
				}
			}
			//attach rest one of the sets not reached it's end
			while(it != iend && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
			{
				*mit = *it;
        m->SetPrivateMarker(*mit,mrk);
				++mit;
				++it;
			}
			while(jt != jend && *jt != InvalidHandle() && !m->GetPrivateMarker(*jt,mrk) )
			{
				*mit = *jt;
        m->SetPrivateMarker(*mit,mrk);
				++mit;
				++jt;
			}
			//indicate end of sorted part
			hSorted(hc) = static_cast<HandleType>(mit-merge.begin());
			//attach unsorted part
			while(it != lc.end() && *it != InvalidHandle() && !m->GetPrivateMarker(*it,mrk) )
			{
				m->SetPrivateMarker(*it,mrk);
				*mit++ = *it++;
			}
			//attach unsorted part
			while(jt != other_lc.end() && *jt != InvalidHandle() && !m->GetPrivateMarker(*jt,mrk) )
			{
				m->SetPrivateMarker(*jt,mrk);
				*mit++ = *jt++;
			}
			merge.resize(static_cast<Element::adj_type::size_type>(mit-merge.begin()));
			m->RemPrivateMarkerArray(merge.data(),static_cast<enumerator>(merge.size()),mrk);
			m->ReleasePrivateMarker(mrk);
			// set merged elemetns as current
			lc.swap(merge);
			//indicate that there are no empty spaces
			hc.resize(high_conn_reserved);
		}
	}

	void ElementSet::Unite(const HandleType * handles, enumerator num) const
	{
		assert(GetElementType() == ESET);
		AddElements(handles,num); //These algorithms are identical
	}


	void ElementSet::Subtract(const ElementSet & other) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->LowConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		other->SetPrivateMarkerElements(mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && m->GetPrivateMarker(lc[i],mrk))
			{
				lc[i] = InvalidHandle();
				hc.push_back(i);
			}
		}
		other->RemPrivateMarkerElements(mrk);
		m->ReleasePrivateMarker(mrk);
	}

	void ElementSet::Subtract(const HandleType * handles, enumerator num) const
	{
		assert(GetElementType() == ESET);
		Mesh * m = GetMeshLink();
		Element::adj_type & lc = m->LowConn(GetHandle());
		Element::adj_type & hc = m->HighConn(GetHandle());
		MarkerType mrk = m->CreatePrivateMarker();
		m->SetPrivateMarkerArray(handles,num,mrk);
		for(Element::adj_type::size_type i = 0; i < lc.size(); ++i)
		{
			if( lc[i] != InvalidHandle() && m->GetPrivateMarker(lc[i],mrk))
			{
				lc[i] = InvalidHandle();
				hc.push_back(i);
			}
		}