mesh.cpp

#include "inmost.h"
#include <sstream>
#if defined(USE_MESH)
#define WAITNL 	{char c;scanf("%c",&c);}

namespace INMOST
{
  static std::vector<Mesh *> allocated_meshes;


#if defined(USE_PARALLEL_WRITE_TIME)
	void Mesh::AtExit(void)
	{
		while(!allocated_meshes.empty())
		{
			if( allocated_meshes.back() != NULL )
			{
				allocated_meshes.back()->FinalizeFile();
			}
			allocated_meshes.pop_back();
		}
	}
#endif //USE_PARALLEL_WRITE_TIME

  std::string Mesh::GetMeshName()
  {
    return name;
  }
  void Mesh::SetMeshName(std::string new_name)
  {
    name = new_name;
  }
  Mesh * Mesh::GetMesh(std::string name)
  {
    for(int q = 0; q < (int)allocated_meshes.size(); ++q)
      if( allocated_meshes[q]->GetMeshName() == name )
        return allocated_meshes[q];
    return NULL;
  }

	const char * TopologyCheckNotifyString(TopologyCheck c)
	{
		switch(c)
		{
			case THROW_EXCEPTION:        return "exception thrown";
			case PRINT_NOTIFY:           return "print notify";
			case DELETE_ON_ERROR:        return "element should be deleted on error";
			case MARK_ON_ERROR:          return "element is marked on error";
			case DUPLICATE_EDGE:         return "search for duplicate edge";
			case DUPLICATE_FACE:         return "search for duplicate face";
			case DUPLICATE_CELL:         return "search for duplicate cell";
			case DEGENERATE_EDGE:        return "TOPOLOGY ERROR: curvilinear edge found"; 
			case DEGENERATE_FACE:        return "TOPOLOGY ERROR: degenerate face found";
			case DEGENERATE_CELL:        return "TOPOLOGY ERROR: degenerate cell found";
			case FACE_ORIENTATION:       return "TOPOLOGY ERROR: bad face orientation";
			case FACE_PLANARITY:         return "TOPOLOGY ERROR: non-planar face found";
			case INTERLEAVED_FACES:      return "TOPOLOGY ERROR: interleaving faces found";
			case TRIPLE_SHARED_FACE:     return "TOPOLOGY ERROR: face have more then two neighbours"; 
			case FLATTENED_CELL:         return "TOPOLOGY ERROR: flattened cell found"; 
			case ADJACENT_DUPLICATE:     return "TOPOLOGY ERROR: duplicates in adjacent elements";
			case ADJACENT_HIDDEN:        return "TOPOLOGY ERROR: hidden element is used as adjacent"; 
			case ADJACENT_VALID:         return "TOPOLOGY ERROR: invalid handle is used as adjacent"; 
			case ADJACENT_DIMENSION:     return "TOPOLOGY ERROR: wrong dimension of adjacent elements";
			case PROHIBIT_MULTILINE:     return "TOPOLOGY ERROR: multiline is prohibited"; 
			case PROHIBIT_POLYGON:       return "TOPOLOGY ERROR: polygon is prohibited"; 
			case PROHIBIT_MULTIPOLYGON:  return "TOPOLOGY ERROR: multipolygon is prohibited"; 
			case PROHIBIT_POLYHEDRON:    return "TOPOLOGY ERROR: polyhedron is prohibited"; 
			case FACE_EDGES_ORDER:       return "TOPOLOGY ERROR: no order in face edges"; 
			case PROHIBIT_CONCAVE_FACE:  return "TOPOLOGY ERROR: concave faces are prohibited"; 
			case PROHIBIT_CONCAVE_CELL:  return "TOPOLOGY ERROR: concave cells are prohibited"; 
			case PROHIBIT_NONSTAR_FACE:  return "TOPOLOGY ERROR: non star-shaped faces are prohibited"; 
			case PROHIBIT_NONSTAR_CELL:  return "TOPOLOGY ERROR: non star-shpaed cells are prohibited"; 
			case FACE_SELF_INTERSECTION: return "TOPOLOGY ERROR: self intersection of face edges detected"; 
			case CELL_SELF_INTERSECTION: return "TOPOLOGY ERROR: self intersection of cell faces detected"; 
			case DISABLE_2D:             return "TOPOLOGY ERROR: 2d mesh support is disabled"; 
			default: return "unknown";
		}
	}

	const char * ElementTypeName(ElementType t)
	{
		switch(t)
		{
			case NONE: return "NONE";
			case NODE: return "NODE";
			case EDGE: return "EDGE";
			case FACE: return "FACE";
			case CELL: return "CELL";
			case ESET: return "ESET";
			case MESH: return "MESH";
		}
		return "UNKNOWN";
	}

  void Mesh::Init(std::string name)
  {
    m_link = this;
		integer selfid = TieElement(5);
		assert(selfid == 0);
		dim = 3;
		have_global_id = NONE;
		checkset = DEFAULT_CHECK;
		errorset = 0;
		new_element = hide_element = 0;

		memset(remember,0,sizeof(remember));
		tag_coords        = CreateTag("PROTECTED_COORD",DATA_REAL, NODE,NONE,dim);
		tag_high_conn     = CreateTag("PROTECTED_HIGH_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_low_conn      = CreateTag("PROTECTED_LOW_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_markers       = CreateTag("PROTECTED_MARKERS",DATA_BULK,CELL|FACE|EDGE|NODE|ESET|MESH,NONE,MarkerFields);
		tag_geom_type     = CreateTag("PROTECTED_GEOM_TYPE",DATA_BULK,CELL|FACE|EDGE|NODE,NONE,1);
		tag_setname       = CreateTag("PROTECTED_SET_NAME",DATA_BULK,ESET,NONE);
		tag_setcomparator = CreateTag("PROTECTED_SET_COMPARATOR",DATA_BULK,ESET,NONE,1);
    AllocatePrivateMarkers();
		for(ElementType etype = NODE; etype <= MESH; etype = etype << 1)
			ReallocateData(ElementNum(etype),GetArrayCapacity(ElementNum(etype)));
		epsilon = 1.0e-8;
		m_state = Mesh::Serial;

#if defined(USE_MPI)
		{
			int test;
			MPI_Initialized(&test);
			if( test == 0 ) MPI_Init(NULL,NULL);
			comm = INMOST_MPI_COMM_WORLD;
		}
#endif

#if defined(USE_PARALLEL_WRITE_TIME)
		num_exchanges = 0;
		std::stringstream temp;
		temp << "time_" << GetProcessorRank() << ".xml";
		out_time.open(temp.str().c_str(),std::ios::out);
		out_time << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>" << std::endl;
		out_time << "<?xml-stylesheet type=\"text/xsl\" href=\"style.xsl\"?>" << std::endl;
		out_time << "<Debug>" << std::endl;
		tab = 1;
		func_id = 0;
#endif
    allocated_meshes.push_back(this);
  }
		
	Mesh::Mesh()
	:TagManager(), Storage(NULL,ComposeHandle(MESH,0))
	{
    std::stringstream tmp;
    tmp << "Mesh" << allocated_meshes.size();
    name = tmp.str();
    Init(name);
	}

  Mesh::Mesh(std::string name)
	:TagManager(), Storage(NULL,ComposeHandle(MESH,0))
	{
    Init(name);
	}

  bool Mesh::GetPrivateMarker(HandleType h, MarkerType n) const
  {
    assert(isPrivate(n));
    n &= ~MarkerPrivateBit;
    int thread = 0;
#if defined(USE_OMP)
    thread = omp_get_thread_num();
#endif
    const bulk * mem = static_cast<const bulk *>(MGetDenseLink(h,tag_private_markers[thread]));
    return (mem[n >> MarkerShift] & static_cast<bulk>(n & MarkerMask)) != 0;
  }

  void Mesh::SetPrivateMarker(HandleType h,MarkerType n)  
  {
    assert(isPrivate(n));
    n &= ~MarkerPrivateBit;
    int thread = 0;
#if defined(USE_OMP)
    thread = omp_get_thread_num();
#endif
    bulk * mem = static_cast<bulk *>(MGetDenseLink(h,tag_private_markers[thread]));
    mem[n >> MarkerShift] |= static_cast<bulk>(n & MarkerMask);
  }

  void Mesh::RemPrivateMarker(HandleType h,MarkerType n) 
  {
    assert(isPrivate(n));
    n &= ~MarkerPrivateBit;
    int thread = 0;
#if defined(USE_OMP)
    thread = omp_get_thread_num();
#endif
    bulk * mem = static_cast<bulk *>(MGetDenseLink(h,tag_private_markers[thread]));
    mem[n >> MarkerShift] &= ~static_cast<bulk>(n & MarkerMask);
  }

  void Mesh::AllocatePrivateMarkers()
  {
#if defined(USE_OMP)
#pragma omp parallel
    {
#pragma omp single
      {
        tag_private_markers = new Tag[omp_get_num_threads()];
      }
//#pragma omp ordered
      {
        std::stringstream name;
        name << "PROTECTED_PRIVATE_MARKERS_" << omp_get_thread_num();
        tag_private_markers[omp_get_thread_num()] = CreateTag(name.str(),DATA_BULK,CELL|FACE|EDGE|NODE|ESET|MESH,NONE,MarkerFieldsPrivate);
      }
    }
#else
    tag_private_markers = new Tag;
    tag_private_markers[0] = CreateTag("PROTECTED_PRIVATE_MARKERS_0",DATA_BULK,CELL|FACE|EDGE|NODE|ESET|MESH,NONE,MarkerFieldsPrivate);
#endif
  }

  void Mesh::DeallocatePrivateMarkers()
  {
#if defined(USE_OMP)
#pragma omp parallel
    {
      //retrive tag before it was erased
      Tag del = tag_private_markers[omp_get_thread_num()];
#pragma omp barrier
      //delete tag, it will erase the tag
      DeleteTag(del);
#pragma omp barrier
      //deallocate space
#pragma omp single
      {
        delete [] tag_private_markers;
      }
    }
#else
    DeleteTag(tag_private_markers[0]);
    delete tag_private_markers;
#endif
  }

	Storage::enumerator Mesh::MemoryUsage(HandleType h)
	{
		if( isValidHandle(h) )
		{
			integer etypenum = GetHandleElementNum(h);
			enumerator ret = 2*sizeof(integer); //link and address occupied
			for(Mesh::iteratorTag t = BeginTag(); t != EndTag(); ++t)
			{
				if( t->isDefinedByDim(etypenum) )
				{
					bool have_data = true;
					if( t->isSparseByDim(etypenum) )
					{
						have_data = HaveData(h,*t);
						if( have_data ) 
							ret += sizeof(sparse_sub_record); //size occupied for storage of data link
					}
					if( have_data )
					{
						ret += t->GetRecordSize(); //all utilized data for fixed data, size of support structure for variable data
						if( t->GetSize() == ENUMUNDEF )
							ret += GetDataCapacity(h,*t); //actually occupied size for sparse data
					}
				}
				if( !sparse_data[etypenum].empty() ) ret += sizeof(sparse_sub_type); //size needed to support sparse data
			}
			//Any additional impact of supporting huge structures over all elements may be added later
			return ret;
		}
		else return 0;
	}
	
	Mesh::Mesh(const Mesh & other)
	:TagManager(other),Storage(NULL,ComposeHandle(MESH,0))
	{
    {
      std::stringstream tmp;
      tmp << other.name << "_copy";
      name = tmp.str();
    }
		m_link = this;
		integer selfid = TieElement(5);
		assert(selfid == 0);
		//TagManager constuctor copied only tags
		//copy links:
		for(int i = 0; i < 5; i++)
		{
			links[i] = other.links[i];
			empty_links[i] = other.empty_links[i];
			empty_space[i] = other.empty_space[i];
		}
		//this should alocate space for data and copy it including markers and connections
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			ReallocateData(ElementNum(etype),GetArrayCapacity(ElementNum(etype)));
			for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
			{
				if( tags[i].isDefined(etype) )
				{
					if( tags[i].isSparse(etype) )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) if( isValidElement(etype,lid) )
						{
							HandleType h = ComposeHandle(etype,lid);
							if( other.HaveData(h,tags[i]) )
								TagManager::CopyData(tags[i],MGetLink(h,tags[i]),other.MGetLink(h,other.tags[i]));
						}
					}
					else
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) if( isValidElement(etype,lid) )
						{
							HandleType h = ComposeHandle(etype,lid);
							TagManager::CopyData(tags[i],MGetLink(h,tags[i]),other.MGetLink(h,other.tags[i]));
						}
					}
				}
			}
		}
		//setup system tags shortcuts
		dim = other.dim;
		tag_coords        = CreateTag("PROTECTED_COORD",DATA_REAL, NODE,NONE,dim);
		tag_high_conn     = CreateTag("PROTECTED_HIGH_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_low_conn      = CreateTag("PROTECTED_LOW_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_markers       = CreateTag("PROTECTED_MARKERS",DATA_BULK,CELL|FACE|EDGE|NODE|ESET|MESH,NONE,MarkerFields);
		tag_geom_type     = CreateTag("PROTECTED_GEOM_TYPE",DATA_BULK,CELL|FACE|EDGE|NODE,NONE,1);
		tag_setname       = CreateTag("PROTECTED_SET_NAME",DATA_BULK,ESET,NONE);
		tag_setcomparator = CreateTag("PROTECTED_SET_COMPARATOR",DATA_BULK,ESET,NONE,1);
    AllocatePrivateMarkers();
		//copy supplimentary values
		m_state = other.m_state;
		checkset = other.checkset;
		errorset = other.errorset;
		new_element = other.new_element;
		hide_element = other.hide_element;
		epsilon = other.epsilon;
		have_global_id = other.have_global_id;
		// copy communicator
		if( m_state == Mesh::Parallel ) SetCommunicator(other.comm); else comm = INMOST_MPI_COMM_WORLD;
		// reestablish geometric tags and table
		memcpy(remember,other.remember,sizeof(remember));
		RestoreGeometricTags();
		//this is not needed as it was copied with all the other data
		//recompute global ids
		//AssignGlobalID(other.have_global_id);
    allocated_meshes.push_back(this);
	}
	
	Mesh & Mesh::operator =(Mesh const & other)
	{
		if( this == &other ) return *this; //don't do anything
    {
      std::stringstream tmp;
      tmp << other.name << "_copy";
      name = tmp.str();
    }
		//first delete everything
		//delete parallel vars
#if defined(USE_MPI)
#if defined(USE_MPI_P2P)
		if( m_state == Mesh::Parallel )
		{
			MPI_Free_mem(shared_space);
			MPI_Win_free(&window);
		}
#endif
#endif
		//clear all data fields
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
			{
				if( tags[i].isDefined(etype) )
				{
					if( tags[i].isSparse(etype) )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelSparseData(ComposeHandle(etype,lid),tags[i]);
					}
					else if( tags[i].GetSize() == ENUMUNDEF )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelDenseData(ComposeHandle(etype,lid),tags[i]);
					}
				}
			}
		}
		//clear links
		for(int i = 0; i < 5; i++)
		{
			links[i].clear();
			empty_links[i].clear();
			empty_space[i].clear();
		}
    DeallocatePrivateMarkers();
		//this should copy tags, clear sparse data, set up dense links
		TagManager::operator =(other);
		//set up new links
		for(int i = 0; i < 5; i++)
		{
			links[i] = other.links[i];
			empty_links[i] = other.empty_links[i];
			empty_space[i] = other.empty_space[i];
		}
		//this should alocate space for data and copy it including markers and connections
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			ReallocateData(ElementNum(etype),GetArrayCapacity(ElementNum(etype)));
			for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
			{
				if( tags[i].isDefined(etype) )
				{
					if( tags[i].isSparse(etype) )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) if( isValidElement(etype,lid) )
						{
							HandleType h = ComposeHandle(etype,lid);
							if( other.HaveData(h,tags[i]) )
								TagManager::CopyData(tags[i],MGetLink(h,tags[i]),other.MGetLink(h,other.tags[i]));
						}
					}
					else
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) if( isValidElement(etype,lid) )
						{
							HandleType h = ComposeHandle(etype,lid);
							TagManager::CopyData(tags[i],MGetLink(h,tags[i]),other.MGetLink(h,other.tags[i]));
						}
					}
				}
			}
		}
		//setup system tags shortcuts
		dim = other.dim;
		tag_coords        = CreateTag("PROTECTED_COORD",DATA_REAL, NODE,NONE,dim);
		tag_high_conn     = CreateTag("PROTECTED_HIGH_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_low_conn      = CreateTag("PROTECTED_LOW_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
		tag_markers       = CreateTag("PROTECTED_MARKERS",DATA_BULK,CELL|FACE|EDGE|NODE|ESET|MESH,NONE,MarkerFields);
		tag_geom_type     = CreateTag("PROTECTED_GEOM_TYPE",DATA_BULK,CELL|FACE|EDGE|NODE,NONE,1);
		tag_setname       = CreateTag("PROTECTED_SET_NAME",DATA_BULK,ESET,NONE);
		tag_setcomparator = CreateTag("PROTECTED_SET_COMPARATOR",DATA_BULK,ESET,NONE,1);
    AllocatePrivateMarkers();
		//copy supplimentary values
		m_state = other.m_state;
		checkset = other.checkset;
		errorset = other.errorset;
		new_element = other.new_element;
		hide_element = other.hide_element;
		epsilon = other.epsilon;
		have_global_id = other.have_global_id;
		// copy communicator
		if( m_state == Mesh::Parallel ) SetCommunicator(other.comm); else comm = INMOST_MPI_COMM_WORLD;
		// reestablish geometric tags and table
		memcpy(remember,other.remember,sizeof(remember));
		RestoreGeometricTags();
		//this is not needed as it was copied with all the other data
		//recompute global ids
		//AssignGlobalID(other.have_global_id);
		return *this;
	}

	void Mesh::Clear()
	{
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
			{
				if( tags[i].isDefined(etype) )
				{
					if( tags[i].isSparse(etype) )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelSparseData(ComposeHandle(etype,lid),tags[i]);
					}
					else if( tags[i].GetSize() == ENUMUNDEF )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelDenseData(ComposeHandle(etype,lid),tags[i]);
					}
				}
			}
		}
    DeallocatePrivateMarkers();
		memset(remember,0,sizeof(bool)*15);
		tags.clear();
		//clear links
		dense_data.clear();
		for(int i = 0; i < 5; i++)
		{
			links[i].clear();
			empty_links[i].clear();
			empty_space[i].clear();
		}
		for(int i = 0; i < 6; i++)
		{
			sparse_data[i].clear();
			back_links[i].clear();
		}
		RemTopologyCheck(ENUMUNDEF);
		SetTopologyCheck(DEFAULT_CHECK);

	}
	
	Mesh::~Mesh()
	{
		//clear all data fields
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
			{
				if( tags[i].isDefined(etype) )
				{
					if( tags[i].isSparse(etype) )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelSparseData(ComposeHandle(etype,lid),tags[i]);
					}
					else if( tags[i].GetSize() == ENUMUNDEF )
					{
						for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
							if( isValidElement(etype,lid) )
								DelDenseData(ComposeHandle(etype,lid),tags[i]);
					}
				}
			}
		}
    DeallocatePrivateMarkers();
		//clear links
		for(int i = 0; i < 5; i++)
		{
			links[i].clear();
			empty_links[i].clear();
			empty_space[i].clear();
		}
#if defined(USE_MPI)
#if defined(USE_MPI_P2P)
		if( m_state == Mesh::Parallel )
		{
			int test = 0;
			MPI_Finalized(&test);
			if( !test )
			{
				MPI_Free_mem(shared_space);
				MPI_Win_free(&window);
			}
			else
			{
				std::cout << "Cannot release memory and window allocated by MPI" << std::endl;
				std::cout << "since MPI was already finalized. Most likely this" << std::endl;
				std::cout << "happens when you define class Mesh in main() function" << std::endl;
				std::cout << "so that destructor get called after Mesh::Finalize()" << std::endl;
				std::cout << "Please enclose area where you use Mesh class" << std::endl;
				std::cout << "with scopes, so that destructor get called when" << std::endl;
				std::cout << "execution goes out of the scope or dynamically" << std::endl;
				std::cout << "allocate and explicitly delete Mesh class." << std::endl;
				std::cout << "Thank you!" << std::endl;
			}
			//~ MPI_Comm_free(&comm);
		}
#endif //USE_MPI_P2P
#endif //USE_MPI
#if defined(USE_PARALLEL_WRITE_TIME)
		FinalizeFile();
#endif //USE_PARALLEL_WRITE_TIME
    {
      for(size_t q = 0; q < allocated_meshes.size(); ++q)
			  if (allocated_meshes[q] == this)
				  allocated_meshes[q] = NULL;
      std::sort(allocated_meshes.rbegin(), allocated_meshes.rend());
      while(!allocated_meshes.empty() && allocated_meshes.back() == NULL) allocated_meshes.pop_back();
    }
		//arrays for data are deallocated inside ~TagManager()
	}
	
	
	
	Tag Mesh::CreateTag(std::string name, DataType dtype, ElementType etype,ElementType sparse, INMOST_DATA_ENUM_TYPE size)
	{
		Tag ret = TagManager::CreateTag(this,name,dtype,etype,sparse,size);
		return ret;
	}
	Tag Mesh::DeleteTag(Tag tag, ElementType type_mask)
	{
		//deallocate data on elements
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
		{
			if( (etype & type_mask) && tag.isDefined(etype) )
			{
				if( tag.isSparse(etype) )
				{
					for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
						if( isValidElement(etype,lid) )
							DelSparseData(ComposeHandle(etype,lid),tag);
				}
				else if( tag.GetSize() == ENUMUNDEF )
				{
					for(integer lid = 0; lid < LastLocalID(etype); ++lid) 
						if( isValidElement(etype,lid) )
							DelDenseData(ComposeHandle(etype,lid),tag);
				}
			}
		}
#if defined(USE_OMP)
#pragma omp critical
#endif
    {
		  tag = TagManager::DeleteTag(tag,type_mask);
    }
		return tag;
	}
	
	
	
	HandleType Mesh::FindSharedAdjacency(const HandleType * arr, enumerator s) const
	{
		if( s == 0 ) return InvalidHandle();
		if( !HideMarker() )
		{
			{
				enumerator flag0, flag1, i;
				dynarray<Element::adj_type const *, 64> hcarr(s);
				hcarr[0] = &HighConn(arr[0]);
				if( !hcarr[0]->empty() ) 
				{
					for(i = 1; i < s; i++) hcarr[i] = &HighConn(arr[i]);
				}
				else return InvalidHandle();
				enumerator it, iend = static_cast<enumerator>(hcarr[0]->size()), jt, jend;
				for(it = 0; it < iend; ++it)
				{
					flag0 = 0;
					for(i = 1; i < s; ++i)
					{
						jend = static_cast<enumerator>(hcarr[i]->size());
						flag1 = 0;
						for(jt = 0; jt < jend; ++jt)
						{
							if( hcarr[0]->at(it) == hcarr[i]->at(jt) )
							{
								flag0++;
								flag1 = 1;
								break;
							}
						}
						if( flag1 == 0 ) break;
					}
					if( flag0 == s-1 ) 
					{
						return hcarr[0]->at(it);
					}
				}
			}
		}
		else
		{
			enumerator ss = Count(arr,s,HideMarker()), nextk;
			enumerator k = ENUMUNDEF, i;
			k = getNext(arr,s,k,HideMarker());
			nextk = getNext(arr,s,k,HideMarker());
			Element::adj_type const & hc0 = HighConn(arr[k]);
			Element::adj_type::size_type it, iend = hc0.size(), jt, jend, flag0, flag1;
			for(it = 0; it < iend; ++it) if( !GetMarker(hc0[it],HideMarker()) )
			{
				flag0 = 0;
				i = nextk;
				while( i < s )
				{
					flag1 = 0;
					Element::adj_type const & ihc = HighConn(arr[i]);
					jend = ihc.size();
					for(jt = 0; jt < jend; ++jt) if( !GetMarker(ihc[jt],HideMarker()) )
					{
						if( hc0[it] == ihc[jt] )
						{
							flag0++;
							flag1 = 1;
							break;
						}
					}
					if( flag1 == 0 ) break;
					i = getNext(arr,s,i,HideMarker());
				}
				if( flag0 == ss-1 ) return hc0[it];
			}
		}
		return InvalidHandle();
	}
	
	
	TopologyCheck Mesh::BeginTopologyCheck(ElementType etype, const HandleType * adj, enumerator s)
	{
		enumerator i,j,d = ENUMUNDEF;
		TopologyCheck chk = 0;
		if( GetTopologyCheck(ADJACENT_VALID) )
		{
			for(i = 0; i < s; i++)
				if( !isValidHandle(adj[i]) )
				{
					chk |= ADJACENT_VALID;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_VALID) << std::endl;
				}
		}
		if( GetTopologyCheck(ADJACENT_HIDDEN) )
		{
			for(i = 0; i < s; i++)
				if( GetMarker(adj[i],HideMarker()) )
				{
					chk |= ADJACENT_HIDDEN;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_HIDDEN) << std::endl;
				}
		}
		if( GetTopologyCheck(ADJACENT_DUPLICATE) )
		{
			bool have_dup = false;
			MarkerType dup = CreateMarker();
			for(i = 0; i < s; i++)
			{
				if( GetMarker(adj[i],dup) ) 
				{
					have_dup = true; // duplication of element
					break;
				}
				else SetMarker(adj[i],dup);
			}
			for(i = 0; i < s; i++) RemMarker(adj[i],dup);
			ReleaseMarker(dup);	
			
			if( have_dup )
			{
				chk |= ADJACENT_DUPLICATE;
				if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DUPLICATE) << std::endl;
			}
		}
		
		if( GetTopologyCheck(ADJACENT_DIMENSION | DEGENERATE_CELL | DEGENERATE_FACE | DEGENERATE_EDGE | DISABLE_2D) )
		{
			bool happen = false;
			for(i = 0; i < s; i++)
			{
				j = ElementByHandle(adj[i])->GetElementDimension();
				if( j == ENUMUNDEF ) 
				{
					chk |= ADJACENT_DIMENSION;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DIMENSION) << std::endl;
					happen = true;
					break;
				}
				if( d == ENUMUNDEF ) d = j;
				if( d != j ) 
				{
					chk |= ADJACENT_DIMENSION;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DIMENSION) << std::endl;
					happen = true;
					break;
				}
			}
			if(!happen) switch(etype) //if happen = true, we cannot be sure about value in d
			{
				case EDGE: 
					if( GetTopologyCheck(DISABLE_2D) )
					{
						if( s == 1 && d == 0 ) 
						{
							chk |= DISABLE_2D;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DISABLE_2D) << std::endl;
						}
					}
					if( GetTopologyCheck(ADJACENT_DIMENSION) )
					{
						if( d > 0 ) 
						{
							chk |= ADJACENT_DIMENSION;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DIMENSION) << std::endl;
						}
					}
					if( GetTopologyCheck(DEGENERATE_EDGE) )
					{
						if( s > 2 )
						{
							chk |= DEGENERATE_EDGE;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DEGENERATE_EDGE) << std::endl;
						}
					}
				break;
				case FACE:
					if( GetTopologyCheck(DISABLE_2D) )
					{
						if( s == 2 && d == 0 ) 
						{
							chk |= DISABLE_2D;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DISABLE_2D) << std::endl;
						}
					}
					if( GetTopologyCheck(ADJACENT_DIMENSION) )
					{
						if( d > 1 ) 
						{
							chk |= ADJACENT_DIMENSION;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DIMENSION) << std::endl;
						}
					}
					if( GetTopologyCheck(DEGENERATE_EDGE) )
					{
						if(d == 0 && s != 2) //This should be Edge (2D)
						{
							chk |= DEGENERATE_EDGE;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DEGENERATE_EDGE) << std::endl;
						}
					}
					if( GetTopologyCheck() & DEGENERATE_FACE )
					{
						if(d == 1 && s < 3) // Should not be less then 3 Line or Curve
						{
							chk |= DEGENERATE_FACE;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DEGENERATE_FACE) << std::endl;
						}					
					}
				break;
				case CELL:
					if( GetTopologyCheck(DISABLE_2D) )
					{
						if( d == 1 ) 
						{
							chk |= DISABLE_2D;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DISABLE_2D) << std::endl;
						}
					}
					if( GetTopologyCheck(ADJACENT_DIMENSION) )
					{
						if( d > 2 || d == 0 ) 
						{
							chk |= ADJACENT_DIMENSION;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(ADJACENT_DIMENSION) << std::endl;
						}
					}
					if( GetTopologyCheck(DEGENERATE_FACE) )
					{
						if( d == 1 && s < 3 ) //Should not be less then 3 Line
						{
							chk |= DEGENERATE_FACE;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DEGENERATE_FACE) << std::endl;
						} 
					}
					if( GetTopologyCheck(DEGENERATE_CELL) )
					{
						if( d == 2 && s < 4 ) //Should not be less then 4 Tri, Quad, Polygon
						{
							chk |= DEGENERATE_CELL;
							if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(DEGENERATE_CELL) << std::endl;
						} 
					}
				break;
			}
		}
		if( etype == CELL )
		{
			if( GetTopologyCheck(TRIPLE_SHARED_FACE) )
			{
				bool check = false;
				for(i = 0; i < s; i++)
				{
					if( ElementByHandle(adj[i])->nbAdjElements(CELL) == 2 )
					{
						check = true;
						break;
					}
				}
				if( check )
				{
					bool happen = true;
					if( (GetTopologyCheck(DUPLICATE_CELL)) && (FindSharedAdjacency(adj,s) != InvalidHandle()) ) happen = false;
					if( happen )
					{
						chk |= TRIPLE_SHARED_FACE;
						if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(TRIPLE_SHARED_FACE) << std::endl;
					}
				}
			}
		}
		else if( etype == FACE )
		{
			if( GetTopologyCheck(INTERLEAVED_FACES) )
			{
				ElementArray<Node> allnodes(this);
				for(i = 0; i < s; i++) allnodes.Unite(ElementByHandle(adj[i])->getNodes());
				ElementArray<Face> faces = allnodes[0].getFaces();
				for(i = 1; i < static_cast<enumerator>(allnodes.size()) && !faces.empty(); i++) faces.Intersect(allnodes[i].getFaces());
				for(i = 0; i < static_cast<enumerator>(faces.size()); i++)
				{
					if( faces[i]->nbAdjElements(NODE) != allnodes.size() )
					{
						chk |= INTERLEAVED_FACES;
						if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(INTERLEAVED_FACES) << std::endl;
					}
				}
			}
		}
		errorset |= chk;
		return chk;
	}
	
	TopologyCheck Mesh::EndTopologyCheck(HandleType he)
	{
		TopologyCheck chk = 0;
		Element e = ElementByHandle(he);
		switch(e->GetGeometricDimension(e->GetGeometricType()))
		{
		case 3:
			if( GetTopologyCheck(FLATTENED_CELL) )
			{
				ElementArray<Face> faces = e->getFaces();
				enumerator num = e->nbAdjElements(NODE);
				for(enumerator i = 0; i < static_cast<enumerator>(faces.size()); i++)
				{
					if( faces[i].nbAdjElements(NODE) == num )
					{
						chk |= FLATTENED_CELL;
						if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(FLATTENED_CELL) << std::endl;
					}
				}
			}
			if( GetTopologyCheck(FACE_ORIENTATION) )
			{
				ElementArray<Face> faces = e->getFaces();
				for(enumerator i = 0; i < static_cast<enumerator>(faces.size()); i++)
				{
					if( faces[i].nbAdjElements(CELL) == 1 && !faces[i].CheckNormalOrientation() )
					{
						chk |= FACE_ORIENTATION;
						if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(FACE_ORIENTATION) << std::endl;
						break;
					}
				}
			}
			break;
		case 2:
			if( GetTopologyCheck(FACE_PLANARITY) )
			{
				if( (e->GetElementType() == FACE && !e->getAsFace()->Planarity()) || (e->GetElementType() == CELL && !e->getAsCell()->Planarity()) )
				{
					chk |= FACE_PLANARITY;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(FACE_PLANARITY) << std::endl;
				}
			}
			if( GetTopologyCheck(FACE_EDGES_ORDER) )
			{
				if( (e->GetElementType() == FACE && !e->getAsFace()->CheckEdgeOrder()) || (e->GetElementType() == CELL && !e->getAsCell()->CheckEdgeOrder()) )
				{
					chk |= FACE_EDGES_ORDER;
					if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(FACE_EDGES_ORDER) << std::endl;
				}
			}
			break;
		}
		if( GetTopologyCheck(PROHIBIT_MULTILINE) )
		{
			if( e->GetGeometricType() == Element::MultiLine )
			{
				chk |= PROHIBIT_MULTILINE;
				if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(PROHIBIT_MULTILINE) << std::endl;
			}
		}
		if( GetTopologyCheck(PROHIBIT_MULTIPOLYGON) )
		{
			if( e->GetGeometricType() == Element::MultiPolygon )
			{
				chk |= PROHIBIT_MULTIPOLYGON;
				if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(PROHIBIT_MULTIPOLYGON) << std::endl;
			}
		}
		if( GetTopologyCheck(PROHIBIT_POLYGON) )
		{
			if( e->GetGeometricType() == Element::Polygon )
			{
				chk |= PROHIBIT_POLYGON;
				if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(PROHIBIT_POLYGON) << std::endl;
			}
		}
		if( GetTopologyCheck(PROHIBIT_POLYHEDRON) )
		{
			if( e->GetGeometricType() == Element::Polyhedron )
			{
				chk |= PROHIBIT_POLYHEDRON;
				if( GetTopologyCheck(PRINT_NOTIFY) ) std::cerr << TopologyCheckNotifyString(PROHIBIT_POLYHEDRON) << std::endl;
			}
		}
		errorset |= chk;
		return chk;
	}
	
	Node Mesh::CreateNode(const real * coords)
	{
		integer id = TieElement(0);
		HandleType h = ComposeHandleNum(0,id);
		SetGeometricType(h,Element::Vertex);
		real * v = static_cast<Storage::real *>(MGetDenseLink(h,CoordsTag()));
		for(integer i = 0; i < dim; i++) v[i] = coords[i];
		SetMarker(h,NewMarker());
		return Node(this,h);
	}

	void PrintHandle(HandleType h)
	{
		std::cout << ElementTypeName(GetHandleElementType(h)) << " " << GetHandleID(h);
	}

	void PrintAdjElem(Element::adj_type const & arr)
	{
		for(Element::adj_type::size_type it = 0; it < arr.size(); ++it)
		{
			PrintHandle(arr[it]);
			std::cout << " ";
		}
		std::cout << std::endl;
	}

	void PrintHandles(const HandleType * beg, const HandleType * end)
	{
		for(const HandleType * it = beg; it != end; ++it)
		{
			PrintHandle(*it);
			std::cout << " ";
		}
		std::cout << std::endl;
	}

	std::pair<Edge,bool> Mesh::CreateEdge(const ElementArray<Node> & nodes)
	{
		HandleType he = InvalidHandle();
		if( !nodes.empty() )
		{
			TopologyCheck chk = 0;
			chk |= BeginTopologyCheck(EDGE,nodes.data(),static_cast<enumerator>(nodes.size()));
			if( chk != 0 )
			{
				if( GetTopologyCheck(THROW_EXCEPTION) ) throw TopologyCheckError;
			}
			if( GetTopologyCheck() & DUPLICATE_EDGE )
			{
				HandleType test = FindSharedAdjacency(nodes.data(),static_cast<enumerator>(nodes.size()));
				if (test != InvalidHandle()) return std::make_pair(Edge(this,test),false);
			}
			integer id = TieElement(1);
			he = ComposeHandleNum(1,id);
			for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++)
			{
				Element::adj_type & hc = HighConn(nodes.at(i));

				//PrintHandle(nodes.at(i)); std::cout << " current: "; PrintAdjElem(hc); 

				hc.push_back(he);

				//std::cout << "add "; PrintHandle(he); std::cout << " result: "; PrintAdjElem(hc);
			}
			Element::adj_type & lc = LowConn(he);

			//PrintHandle(he); std::cout << " current: "; PrintAdjElem(lc);
			
			lc.insert(lc.end(),nodes.data(),nodes.data()+nodes.size());

			//std::cout << "add "; PrintHandles(nodes.data(),nodes.data()+nodes.size());
			//std::cout << "result: "; PrintAdjElem(lc);
			//DEBUG
			/*
			bool halt = false;
			if(!Element(this,he)->CheckElementConnectivity()) 
			{
				Element(this,he)->PrintElementConnectivity();
				halt = true;
			}
			for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++) 
				if(!nodes[i]->CheckElementConnectivity()) 
				{
					nodes[i]->PrintElementConnectivity();
					halt = true;
				}
			assert(!halt);
			*/
			//DEBUG END
			ComputeGeometricType(he);
			SetMarker(he,NewMarker());
			RecomputeGeometricData(he);
			chk |= EndTopologyCheck(he);
			if( chk != 0 )
			{
				if( GetTopologyCheck(MARK_ON_ERROR) ) Integer(he,TopologyErrorTag()) = chk;
				if( GetTopologyCheck(DELETE_ON_ERROR) ) { Destroy(he); he = InvalidHandle();}
				if( GetTopologyCheck(THROW_EXCEPTION) ) throw TopologyCheckError;
			}
		}
		return std::make_pair(Edge(this,he),true);
	}
	
	
	std::pair<Face,bool> Mesh::CreateFace(const ElementArray<Node> & f_nodes)
	{
		ElementArray<Edge> f_edges(this,f_nodes.size());
		ElementArray<Node> e_nodes(this,2);
		if( f_nodes.size() == 2 ) //This is an edge for 2d!
		{
			e_nodes.resize(1);
			e_nodes.at(0) = f_nodes.at(0);
			f_edges.at(0) = CreateEdge(e_nodes).first->GetHandle();
			e_nodes.at(0) = f_nodes.at(1);
			f_edges.at(1) = CreateEdge(e_nodes).first->GetHandle();
		}
		else
		{
			for(ElementArray<Node>::size_type i = 0; i < f_nodes.size(); i++)
			{
				e_nodes.at(0) = f_nodes.at(i);
				e_nodes.at(1) = f_nodes.at((i+1)%f_nodes.size());
				f_edges.at(i) = CreateEdge(e_nodes).first->GetHandle();
			}
		}
		return CreateFace(f_edges);
	}
	
	std::pair<Face,bool> Mesh::CreateFace(const ElementArray<Edge> & f_edges)
	{
		HandleType he = InvalidHandle();
		if( !f_edges.empty() )
		{
			TopologyCheck chk = 0;
			chk |= BeginTopologyCheck(FACE,f_edges.data(),static_cast<enumerator>(f_edges.size()));
			if( chk != 0 )
			{
				if( GetTopologyCheck(THROW_EXCEPTION) ) throw TopologyCheckError;
			}
			if( GetTopologyCheck(DUPLICATE_FACE) )
			{
				HandleType test = FindSharedAdjacency(f_edges.data(),static_cast<enumerator>(f_edges.size()));
				if (test != InvalidHandle()) return std::make_pair(Face(this,test),false);
			}
			integer id = TieElement(2);
			he = ComposeHandleNum(2,id);
			for(ElementArray<Edge>::size_type i = 0; i < f_edges.size(); i++)
			{
				Element::adj_type & hc = HighConn(f_edges.at(i));

				//PrintHandle(f_edges.at(i)); std::cout << " current: "; PrintAdjElem(hc); 

				hc.push_back(he);

				//std::cout << "add "; PrintHandle(he); std::cout << " result: "; PrintAdjElem(hc);
			}
			Element::adj_type & lc = LowConn(he);

			//PrintHandle(he); std::cout << "current: "; PrintAdjElem(lc);

			lc.insert(lc.end(),f_edges.data(),f_edges.data()+f_edges.size());

			//std::cout << "add "; PrintHandles(f_edges.data(),f_edges.data()+f_edges.size());
			//std::cout << "result: "; PrintAdjElem(lc);
			//DEBUG
			/*
			bool halt = false;
			if(!Element(this,he)->CheckElementConnectivity())
			{
				Element(this,he)->PrintElementConnectivity();
				halt = true;
			}
			for(ElementArray<Edge>::size_type i = 0; i < f_edges.size(); i++) 
			{
				if(!f_edges[i]->CheckElementConnectivity())
				{
					f_edges[i]->PrintElementConnectivity();
					halt = true;
				}
			}
			assert(!halt);
			*/
			//DEBUG END
			ComputeGeometricType(he);
			SetMarker(he,NewMarker());
			RecomputeGeometricData(he);
			chk |= EndTopologyCheck(he);
			if( chk != 0 )
			{
				if( GetTopologyCheck(MARK_ON_ERROR)   ) Integer(he,TopologyErrorTag()) = chk;
				if( GetTopologyCheck(DELETE_ON_ERROR) ) { Destroy(he); he = InvalidHandle();}
				if( GetTopologyCheck(THROW_EXCEPTION) ) throw TopologyCheckError;
			}
		}
		return std::make_pair(Face(this,he),true);
	}
	
	
	std::pair<Cell,bool> Mesh::CreateCell(const ElementArray<Node> & c_f_nodes, const integer * c_f_sizes, integer s, const ElementArray<Node> & suggest_nodes_order)
	{
		ElementArray<Face> c_faces(this,s);
		ElementArray<Node>::size_type j = 0;
		for(integer i = 0; i < s; i++)
		{
			c_faces.at(i) = CreateFace(ElementArray<Node>(this, c_f_nodes.data()+j, c_f_nodes.data()+j + c_f_sizes[i])).first->GetHandle();
			j += c_f_sizes[i];
		}
		return CreateCell(c_faces,suggest_nodes_order);
	}
	

	std::pair<Cell, bool> Mesh::CreateCell(const ElementArray<Node> & c_f_nodes, const integer * c_f_nodeinds, const integer * c_f_numnodes, integer s, const ElementArray<Node> & suggest_nodes_order)
	{
		integer j = 0;
		ElementArray<Node> temp(this);
		ElementArray<Face> c_faces(this,s);
		for(integer i = 0; i < s; i++)
		{
			temp.resize(c_f_numnodes[i]);
			for(integer k = j; k < j+c_f_numnodes[i]; k++)
				temp.at(k-j) = c_f_nodes.at(c_f_nodeinds[k]);
			c_faces.at(i) = CreateFace(temp).first->GetHandle();
			j += c_f_numnodes[i];
			temp.clear();
		}
		return CreateCell(c_faces,suggest_nodes_order);
	}
	
	
	void Mesh::RestoreCellNodes(HandleType hc, ElementArray<Node> & ret)
	{
		ret.clear();
		Cell c(this,hc);
		switch(GetGeometricType(hc))
		{
			case Element::Vertex:
			case Element::Line:
			{
				assert( false );
			}
			case Element::MultiLine:
			case Element::Tri:
			case Element::Quad:
			case Element::Polygon:
			{
				ElementArray<Edge> edges = c->getEdges();
				ret.reserve(edges.size());
				for(ElementArray<Edge>::iterator it = edges.begin(); it != edges.end(); it++)
				{
					ElementArray<Node> nodes = it->getNodes();
					assert( nodes.size() == 1 );
					ret.push_back(nodes.data()[0]);
				}
				break;
			}
				/*
				 6 7
				 4 5
				 2 3
				 0 1
				 */
			case Element::Hex:
			{
				MarkerType mrk = CreateMarker();
				MarkerType cemrk = CreateMarker();
				MarkerType femrk = CreateMarker();
				//printf("%lx %lx %lx\n",mrk,cemrk,femrk);
				ElementArray<Face> faces = c->getFaces();
				Face face = faces[0];
				ret.reserve(8);
				ElementArray<Node> verts = face->getNodes();
				if( face->BackCell() == c )
					for(ElementArray<Node>::iterator it = verts.begin(); it != verts.end(); it++)
					{
						ret.push_back(*it);
						it->SetMarker(mrk);
					}
				else
					for(ElementArray<Node>::reverse_iterator it = verts.rbegin(); it != verts.rend(); it++)
					{
						ret.push_back(*it);
						it->SetMarker(mrk);
					}
				ElementArray<Edge> c_edges = c->getEdges();
				for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++)
					it->SetMarker(cemrk);
				ElementArray<Edge> f_edges = face->getEdges();
				for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++)
					it->SetMarker(femrk);
				for(unsigned int k = 0; k < 4; k++)
				{
					ElementArray<Edge> v_edges = ret[k]->getEdges();
					for(ElementArray<Edge>::iterator it = v_edges.begin(); it != v_edges.end(); it++)
					{
						if( it->GetMarker(cemrk) && !it->GetMarker(femrk) )
						{
							ElementArray<Node> nn = it->getNodes();
							if( nn[0].GetMarker(mrk) )
								ret.push_back(nn[1]);
							else
								ret.push_back(nn[0]);
							break;
						}
					}
					ret[k]->RemMarker(mrk);
				}
				for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->RemMarker(cemrk);
				for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->RemMarker(femrk);
				ReleaseMarker(mrk);
				ReleaseMarker(cemrk);
				ReleaseMarker(femrk);
				break;
			}
				/*
				 5
				 3 4
				 2 
				 0 1
				 */
			case Element::Prism:
			{
				MarkerType mrk = CreateMarker();
				MarkerType cemrk = CreateMarker();
				MarkerType femrk = CreateMarker();
				ret.reserve(6);
				Face face;
				ElementArray<Face> faces = c->getFaces();
				for(ElementArray<Face>::size_type i = 0; i < faces.size(); i++) //iterate over faces
					if( faces[i].nbAdjElements(EDGE) == 3 ) //number of edges in i-th face
					{
						face = faces[i];
						break;
					}
				ElementArray<Node> verts = face->getNodes();
				if( face->BackCell() == c )
					for(ElementArray<Node>::iterator it = verts.begin(); it != verts.end(); it++)
					{
						ret.push_back(*it);
						it->SetMarker(mrk);
					}
				else
					for(ElementArray<Node>::reverse_iterator it = verts.rbegin(); it != verts.rend(); it++)
					{
						ret.push_back(*it);
						it->SetMarker(mrk);
					}
				ElementArray<Edge> c_edges = c->getEdges();
				for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->SetMarker(cemrk);
				ElementArray<Edge> f_edges = face->getEdges();
				for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->SetMarker(femrk);
				for(ElementArray<Cell>::size_type k = 0; k < 3; k++)
				{
					ElementArray<Edge> v_edges = ret[k]->getEdges();
					for(ElementArray<Edge>::iterator it = v_edges.begin(); it != v_edges.end(); it++)
						if( it->GetMarker(cemrk) && !it->GetMarker(femrk) )
						{
							ElementArray<Node> nn = it->getNodes();
							if( nn[0].GetMarker(mrk) )
								ret.push_back(nn[1]);
							else
								ret.push_back(nn[0]);
							break;
						}
					ret[k]->RemMarker(mrk);
				}
				for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->RemMarker(cemrk);
				for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->RemMarker(femrk);
				ReleaseMarker(mrk);
				ReleaseMarker(cemrk);
				ReleaseMarker(femrk);
				break;
			}
				/*
				  4
				 2 3
				 0 1
				 */
			case Element::Pyramid:
			{
				ret.reserve(5);