mesh_parallel.cpp

#ifdef _MSC_VER //kill some warnings
#define _SCL_SECURE_NO_WARNINGS
#endif

#include "inmost.h"


#if defined(USE_MESH)
#include <iostream>
#include <fstream>
#include <sstream>


#if defined(USE_MPI)
static INMOST_DATA_BIG_ENUM_TYPE pmid = 0;
#endif

#if defined(USE_PARALLEL_WRITE_TIME)
#define REPORT_MPI(x) {WriteTab(out_time) << "<MPI><![CDATA[" << #x << "]]></MPI>" << std::endl; x;}
#define REPORT_STR(x) {WriteTab(out_time) << "<TEXT><![CDATA[" << x << "]]></TEXT>" << std::endl;}
#define REPORT_VAL(str,x) {WriteTab(out_time) << "<VALUE name=\"" << str << "\"> <CONTENT><![CDATA[" << x << "]]></CONTENT> <CODE><![CDATA[" << #x << "]]></CODE></VALUE>" << std::endl;}
#define ENTER_FUNC() double all_time = Timer(); {WriteTab(out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << func_id++ << "\">" << std::endl; Enter();}
#define EXIT_FUNC() {WriteTab(out_time) << "<TIME>" << Timer() - all_time << "</TIME>" << std::endl; Exit(); WriteTab(out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC_DIE() {WriteTab(out_time) << "<TIME>" << -1 << "</TIME>" << std::endl; Exit(); WriteTab(out_time) << "</FUNCTION>" << std::endl;}
#else
#define REPORT_MPI(x) x
#define REPORT_STR(x) {}
#define REPORT_VAL(str,x) {}
#define ENTER_FUNC() {}
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
#endif

namespace INMOST
{
	//////////////////////////////
	/// REDUCTION FUNCTIONS    ///
	//////////////////////////////
	void DefaultUnpack(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size == 0 )
		{
			if( tag.isDefined(element->GetElementType()) ) 
			{
				if( tag.isSparse(element->GetElementType()) )
				{
					if( element->HaveData(tag) )
						element->DelData(tag); 
				}
				else if( tag.GetSize() == ENUMUNDEF )
					element->SetDataSize(tag,size);
				else
				{
					std::cerr << "received zero size for dense nonzero tag" << std::endl;
					assert(false);
				}
			}
			return;
		}
		if( !element->HaveData(tag) )
			element->SetDataSize(tag,size);
		else if( size != element->GetDataSize(tag) )
		{
			if( tag.GetSize() == ENUMUNDEF )
				element->SetDataSize(tag,size);
			else
			{
				assert(false);
			}
		}
		element->SetData(tag,0,size,data);
	}


	void UnpackOnSkin(const Tag & tag, const Element & e, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size )
		{
      const Storage::integer * recv = static_cast<const Storage::integer *>(static_cast<const void *>(data));
			Storage::integer_array arr = e->IntegerArray(tag);
			arr.push_back(recv[0]);
		}
	}


	void UnpackSkin(const Tag & tag, const Element & e, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size )
		{
      for(INMOST_DATA_ENUM_TYPE k = 0; k < size; k+=2 )
      {
			  bool flag = true;
			  const Storage::integer * recv = static_cast<const Storage::integer *>(static_cast<const void *>(data))+k;
			  Storage::integer_array arr = e->IntegerArray(tag);
      
        for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it+=2)
        {
          if( *it == recv[0] )
				  {
					  flag = false;
					  break;
				  }
        }
			  if( flag ) 
			  {
				  arr.push_back(recv[0]);
				  arr.push_back(recv[1]);
			  }
      }
		}
	}
	


	void DeleteUnpack(const Tag & tag, const Element & e, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size ) 
		{
			int old_size;
			if( e->HaveData(tag) ) old_size = e->GetDataSize(tag);
			else old_size = 0;
			e->SetDataSize(tag,old_size+size);
			e->SetData(tag,old_size,size,data);
		}
	}
	
	void RedistUnpack(const Tag & tag,const Element & e, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size ) 
		{
			Storage::integer_array p1 = e->IntegerArray(tag);
			const Storage::integer * p2 = static_cast<const Storage::integer *>(static_cast<const void *>(data));
			dynarray<Storage::integer,64> result(p1.size()+size);
			dynarray<Storage::integer,64>::iterator end;
			end = std::set_union(p1.begin(),p1.end(),p2,p2+size,result.begin());
			result.resize(end-result.begin());
			p1.clear();
			p1.insert(p1.end(),result.begin(),result.end());
		}
	}

	void UnpackLayersMarker(const Tag & tag, const Element & e, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		if( size )
		{
			int old_size;
			if( e->HaveData(tag) ) old_size = e->GetDataSize(tag);
			else old_size = 0;
			e->SetDataSize(tag,old_size+size);
			e->SetData(tag,old_size,size,data);
		}
	}
	

	void UnpackSyncMarkerOR(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Bulk(tag) |= *data;
	}

	void UnpackSyncMarkerXOR(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Bulk(tag) ^= *data;
	}

	
	void UnpackSyncMarkerAND(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Bulk(tag) &= *data;
	}

	//////////////////////////////
	/// REDUCTION FUNCTIONS END///
	//////////////////////////////



	void Mesh::SynchronizeMarker(MarkerType marker, ElementType mask, SyncBitOp op)
	{
#if defined(USE_MPI)
		if( m_state == Mesh::Parallel )
		{
			Tag t = CreateTag("TEMP_SYNC_MARKER",DATA_BULK,mask,mask,1);

			//workaround for old gcc compiler
			const Element::Status SGhost = Element::Ghost;
			const Element::Status SAny = Element::Any;
			Element::Status Expr = (Element::Shared | ((op != SYNC_BIT_SET) ? SGhost : SAny));

			for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
				if( it->GetMarker(marker) && (it->GetStatus() & Expr) )
					it->Bulk(t) = 1;

			
			
			switch(op)
			{
			case SYNC_BIT_SET:
				ExchangeData(t,mask,0);
				for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
				{
					if( it->GetStatus() == Element::Ghost )
					{
						if( it->HaveData(t) ) it->SetMarker(marker); else it->RemMarker(marker);
					}
				}
				break;
			case SYNC_BIT_OR:
				ReduceData(t,mask,0,UnpackSyncMarkerOR);
				for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
				{
					if( it->GetStatus() & (Element::Ghost | Element::Shared) )
					{
						if( !it->GetMarker(marker) && it->HaveData(t) ) it->SetMarker(marker);
					}
				}
				break;
			case SYNC_BIT_AND:
				ReduceData(t,mask,0,UnpackSyncMarkerAND);
				for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
				{
					if( it->GetStatus() & (Element::Ghost | Element::Shared))
					{
						if( it->HaveData(t) && it->Bulk(t) ) it->SetMarker(marker); else it->RemMarker(marker);
					}
				}
				break;
			case SYNC_BIT_XOR:
				ReduceData(t,mask,0,UnpackSyncMarkerXOR);
				for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
				{
					if( it->GetStatus() & (Element::Ghost | Element::Shared))
					{
						if( it->HaveData(t) && it->Bulk(t) ) it->SetMarker(marker); else it->RemMarker(marker);
					}
				}
				break;
			}
			

			DeleteTag(t,mask);
		}
#else//USE_MPI
		(void) marker;
		(void) mask;
		(void) op;
#endif//USE_MPI
	}

	ElementType Mesh::SynchronizeElementType(ElementType etype)
	{
		ElementType etypeout = etype;
#if defined(USE_MPI)
		MPI_Allreduce(&etype,&etypeout,1,INMOST_MPI_DATA_BULK_TYPE,MPI_BOR,GetCommunicator());
#endif//USE_MPI
		return etypeout;
	}
	Storage::integer Mesh::TotalNumberOf(ElementType mask)
	{
		Storage::integer number = 0, ret = 0;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			if( it->GetStatus() != Element::Ghost ) number++;
		ret = number;
#if defined(USE_MPI)
		MPI_Allreduce(&number,&ret,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
#endif//USE_MPI
		return ret;
	}
	Storage::integer Mesh::EnumerateSet(const ElementSet & set, const Tag & num_tag, Storage::integer start, bool define_sparse)
	{
		Storage::integer shift = 0, ret = 0;
		ElementType mask = CELL | FACE | EDGE | NODE;
#if defined(USE_MPI)
		Storage::integer number = 0;
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() != Element::Ghost && (define_sparse || it->HaveData(num_tag)) )
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
#endif//USE_MPI
		shift += start;
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() == Element::Owned && (define_sparse || it->HaveData(num_tag)) )
				it->Integer(num_tag) = shift++;
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() == Element::Shared && (define_sparse || it->HaveData(num_tag)) )
				it->Integer(num_tag) = shift++;
		ExchangeData(num_tag,mask,0);
		ret = shift;
#if defined(USE_MPI)
		MPI_Bcast(&ret,1,INMOST_MPI_DATA_INTEGER_TYPE,GetProcessorsNumber()-1,comm);
		//MPI_Allreduce(&shift,&ret,1,INMOST_DATA_INTEGER_TYPE,MPI_MAX,comm);
#endif//USE_MPI
		return ret;
	}
	Storage::integer Mesh::Enumerate(const HandleType * set, enumerator n, const Tag & num_tag, Storage::integer start, bool define_sparse)
	{
		Storage::integer shift = 0, ret = 0;
		ElementType mask = CELL | FACE | EDGE | NODE;
#if defined(USE_MPI)
		Storage::integer number = 0;
		for(const HandleType * it = set; it != set+n; ++it)
			if( GetStatus(*it) != Element::Ghost && (define_sparse || HaveData(*it,num_tag))) number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
#endif//USE_MPI
		shift += start;
		for(const HandleType * it = set; it != set+n; ++it)
			if( GetStatus(*it) == Element::Owned && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++;
		for(const HandleType * it = set; it != set+n; ++it) 
			if( GetStatus(*it) == Element::Shared && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++;
		ExchangeData(num_tag,mask,0);
		ret = shift;
#if defined(USE_MPI)
		MPI_Bcast(&ret,1,INMOST_MPI_DATA_INTEGER_TYPE,GetProcessorsNumber()-1,comm);
		//MPI_Allreduce(&shift,&ret,1,INMOST_DATA_INTEGER_TYPE,MPI_MAX,comm);
#endif//USE_MPI
		return ret;
	}
	Storage::integer Mesh::Enumerate(ElementType mask, Tag num_tag, Storage::integer start, bool define_sparse)
	{
		Storage::integer shift = 0, ret = 0;
#if defined(USE_MPI)
		Storage::integer number = 0;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			if( it->GetStatus() != Element::Ghost && (define_sparse || it->HaveData(num_tag)) )
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
#endif//USE_MPI
		shift += start;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			if( it->GetStatus() == Element::Owned && (define_sparse || it->HaveData(num_tag)) )
				it->Integer(num_tag) = shift++;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			if( it->GetStatus() == Element::Shared && (define_sparse || it->HaveData(num_tag)) )
				it->Integer(num_tag) = shift++;
		ExchangeData(num_tag,mask,0);
		ret = shift;
#if defined(USE_MPI)
		MPI_Bcast(&ret,1,INMOST_MPI_DATA_INTEGER_TYPE,GetProcessorsNumber()-1,comm);
		//MPI_Allreduce(&shift,&ret,1,INMOST_DATA_INTEGER_TYPE,MPI_MAX,comm);
#endif//USE_MPI
		return ret;
	}
	
	Storage::real Mesh::Integrate(Storage::real input)
	{
		Storage::real output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		return output;
	}


	
	
	Storage::integer Mesh::Integrate(Storage::integer input)
	{
		Storage::integer output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		return output;
	}

	void Mesh::Integrate(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::Integrate(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}
	
	Storage::integer Mesh::ExclusiveSum(Storage::integer input)
	{
		Storage::integer output = 0;
#if defined(USE_MPI)
		MPI_Scan(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		output -= input;
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		return output;
	}
	
	Storage::real Mesh::Integrate(const Tag & t, enumerator entry, ElementType mask)
	{
		Storage::real output = 0, input = 0;
		for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			if( GetStatus(*it) != Element::Ghost && HaveData(*it,t) ) 
			{
				real_array arr = RealArray(*it,t);
				if( arr.size() > entry ) input += arr[entry];
			}
		output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
#endif//USE_MPI
		return output;
	}
	
	Storage::real Mesh::AggregateMax(Storage::real input)
	{
		Storage::real output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm);
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		return output;
	}
	
	Storage::integer Mesh::AggregateMax(Storage::integer input)
	{
		Storage::integer output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm);
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		return output;
	}

	void Mesh::AggregateMax(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::AggregateMax(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	Storage::real Mesh::AggregateMin(Storage::real input)
	{
		Storage::real output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_MIN,comm);
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		return output;
	}
	
	Storage::integer Mesh::AggregateMin(Storage::integer input)
	{
		Storage::integer output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MIN,comm);
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		return output;
	}

	void Mesh::AggregateMin(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MIN,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::AggregateMin(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MIN,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}


	INMOST_MPI_Comm Mesh::GetCommunicator()
	{
#if defined(USE_MPI)
		return comm;
#else //USE_MPI
		return 0;
#endif //USE_MPI
	}
	
	int Mesh::GetProcessorRank()
	{
#if defined(USE_MPI)
		int rank;
		MPI_Comm_rank(comm,&rank);
		return rank;
#else //USE_MPI
		return 0;
#endif //USE_MPI
	}
	
	int Mesh::GetProcessorsNumber()
	{
#if defined(USE_MPI)
		int size;
		MPI_Comm_size(comm,&size);
		return size;
#else //USE_MPI
		return 1;
#endif //USE_MPI
	}	
	
	void Mesh::Initialize(int * argc, char *** argv)
	{
#if defined(USE_MPI)
		int test;
		MPI_Initialized(&test);
        if( test == 0 ) MPI_Init(argc,argv);
#else //USE_MPI
		(void) argc;
		(void) argv;
#endif //USE_MPI
#if defined(USE_PARALLEL_WRITE_TIME)
		atexit(Mesh::AtExit);
#endif
	}
	void Mesh::Finalize()
	{
#if defined(USE_MPI)
		int test = 0;
		MPI_Finalized(&test);
		if( !test )
			MPI_Finalize();
#endif //USE_MPI
	}
	
	
	
	void Mesh::SetCommunicator(INMOST_MPI_Comm _comm)
	{
		ENTER_FUNC();
		tag_shared = CreateTag("PROTECTED_STATUS",DATA_BULK,CELL | FACE | EDGE | NODE,NONE,1);
		tag_owner = CreateTag("OWNER_PROCESSOR",DATA_INTEGER, CELL | FACE | EDGE | NODE,NONE,1);
		tag_processors = CreateTag("PROCESSORS_LIST",DATA_INTEGER, MESH | NODE | EDGE | FACE | CELL,NONE);
		tag_layers = CreateTag("LAYERS",DATA_INTEGER,MESH,NONE,1);
		tag_bridge = CreateTag("BRIDGE",DATA_INTEGER,MESH,NONE,1);
		tag_sendto = CreateTag("PROTECTED_SENDTO",DATA_INTEGER, CELL | FACE | EDGE | NODE, CELL | FACE | EDGE | NODE);

#if defined(USE_MPI)
		
		
		parallel_strategy = 1;
		parallel_file_strategy = 1;

		

		Mesh::Initialize(NULL,NULL);
		//~ MPI_Comm_dup(_comm,&comm);
		comm = _comm;
		{
			INMOST_DATA_BIG_ENUM_TYPE t = pmid;
			REPORT_MPI(MPI_Allreduce(&t,&parallel_mesh_unique_id,1,INMOST_MPI_DATA_BIG_ENUM_TYPE,MPI_MAX,MPI_COMM_WORLD));
			pmid = parallel_mesh_unique_id+1;
		}
		m_state = Mesh::Parallel;


#if defined(USE_MPI_P2P)
		int err;
		REPORT_MPI(err = MPI_Alloc_mem(GetProcessorsNumber()*sizeof(unsigned)*2,MPI_INFO_NULL,&shared_space));
		if( err )
		{
			int errclass;
			MPI_Error_class(err,&errclass);
			std::cout << "Cannot allocate shared space of size " << GetProcessorsNumber()*sizeof(unsigned)*2 << " reason is "; 
			switch(err)
			{
				case MPI_SUCCESS: std::cout << "success"; break;
				case MPI_ERR_INFO: std::cout << "bad info"; break;
				case MPI_ERR_ARG: std::cout << "bad argument"; break;
				case MPI_ERR_NO_MEM: std::cout << "no memory"; break;
			}
			std::cout << std::endl;
			MPI_Abort(comm,err);
		}
		REPORT_MPI(MPI_Win_create(shared_space,sizeof(unsigned)*GetProcessorsNumber()*2,sizeof(unsigned),MPI_INFO_NULL,comm,&window));
#endif //USE_MPI_P2P
#else //USE_MPI
		(void) _comm;
#endif //USE_MPI
		EXIT_FUNC();
	}
	
#if defined(USE_PARALLEL_WRITE_TIME)
	void Mesh::Enter() { tab++; }
	void Mesh::Exit() {tab--; }
	std::fstream & Mesh::WriteTab(std::fstream & f)
	{
		for(int i = 0; i < tab; i++)
			f << "   ";
		return f;
	}
	std::fstream & Mesh::GetStream()
	{
		return out_time;
	}
	void Mesh::FinalizeFile()
	{
		if( tab > 1 )
		{
			REPORT_STR("Finalizing file with nonempty function stack - probably due to error.");
		}
		while(tab > 1)
		{
			EXIT_FUNC_DIE();
		}
		out_time << "</Debug>" << std::endl;
	}
#endif //USE_PARALLEL_WRITE_TIME
	
	
	
	void determine_my_procs_low(Mesh * m, HandleType h, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
	{
		Element::adj_type const & subelements = m->LowConn(h);
		Element::adj_type::const_iterator i = subelements.begin();
		Storage::integer_array p = m->IntegerArrayDV(*i,m->ProcessorsTag());
		result.clear();
		result.insert(result.end(),p.begin(),p.end());
		i++;
		while(i != subelements.end())
		{
			Storage::integer_array q = m->IntegerArrayDV(*i,m->ProcessorsTag());
			intersection.resize(std::max(static_cast<unsigned>(result.size()),static_cast<unsigned>(q.size())));
			dynarray<Storage::integer,64>::iterator qt = std::set_intersection(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
			intersection.resize(qt-intersection.begin());
			result.swap(intersection);
			i++;
		}
	}
	
	void determine_my_procs_high(Mesh * m, HandleType h, const Tag & procs, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
	{
		Element::adj_type const & overelements = m->HighConn(h);
		if( overelements.empty() ) return;
		Element::adj_type::const_iterator i = overelements.begin();
		result.clear();
		while(i != overelements.end())
		{
			Storage::integer_array q = m->IntegerArrayDV(*i,procs);
			intersection.resize(result.size()+q.size());
			dynarray<Storage::integer,64>::iterator qt = std::set_union(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
			intersection.resize(qt-intersection.begin());
			result.swap(intersection);
			i++;
		}
	}
	
	/*
	void determine_my_procs_low(adjacent<Element> & subelements,Storage::integer_array & my_procs,Tag & procs_tag)
	{
		if( subelements.empty() ) return;
		adjacent<Element>::iterator i = subelements.begin();
		std::vector<Storage::integer> result;
		Storage::integer_array p = i->IntegerArray(procs_tag);
		result.insert(result.begin(),p.begin(),p.end());
		i++;
		while(i != subelements.end())
		{
			Storage::integer_array q = i->IntegerArray(procs_tag);
			std::vector<Storage::integer> intersection(result.size());
			std::vector<Storage::integer>::iterator qt = std::set_intersection(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
			intersection.resize(qt-intersection.begin());
			result.swap(intersection);
			i++;
		}
		my_procs.clear();
		my_procs.insert(my_procs.begin(),result.begin(),result.end());
	}
	
	void determine_my_procs_high(adjacent<Element> & overelements,Storage::integer_array & my_procs,Tag & procs_tag)
	{
		if( overelements.empty() ) return;
		adjacent<Element>::iterator i = overelements.begin();
		std::vector<Storage::integer> result;
		Storage::integer_array p = i->IntegerArray(procs_tag);
		result.insert(result.begin(),p.begin(),p.end());
		i++;
		while(i != overelements.end())
		{
			Storage::integer_array q = i->IntegerArray(procs_tag);
			std::vector<Storage::integer> intersection(result.size()+q.size());
			std::vector<Storage::integer>::iterator qt = std::set_union(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
			intersection.resize(qt-intersection.begin());
			result.swap(intersection);
			i++;
		}
		my_procs.clear();
		my_procs.insert(my_procs.begin(),result.begin(),result.end());
	}
	*/
	/*
	void Mesh::MarkShared(ElementType mask)
	{
		for(ElementType etype = NODE; etype <= CELL; etype = etype << 1 ) if( etype & mask )
		for(Mesh::iteratorElement it = BeginElement(etype); it != EndElement(); it++)
		{
			int owner;
			Storage::integer_array v = it->IntegerArrayDV(tag_processors);
			std::sort(v.begin(),v.end());
			if( v.empty() )
			{
				owner = mpirank;
				v.push_back(mpirank);
			}
			else
				owner = std::min(mpirank,v[0]);
			
			it->IntegerDF(tag_owner) = owner;
			
			if( mpirank == owner )
			{
				if( v.size() == 1 )
					it->BulkDF(tag_shared) = Element::Owned;
				else
					it->BulkDF(tag_shared) = Element::Shared;
			}
			else
				it->BulkDF(tag_shared) = Element::Ghost;
		}
		
		ComputeSharedProcs();
		RecomputeParallelStorage(mask);
		AssignGlobalID(mask);
		
		//have to do it for all types
#if defined(USE_PARALLEL_STORAGE)
		for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
			qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
		for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)			
			qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
#endif
	}
	*/
	
	
	__INLINE bool point_in_bbox(Storage::real * p, Storage::real bbox[6], unsigned dim)
	{
		bool ret = true;
		for(unsigned k = 0; k < dim; k++) 
			ret &= (p[k] >= bbox[k] && p[k] <= bbox[dim+k]);
		return ret;
	}
	
	__INLINE bool compare_coord(Storage::real * a, Storage::real * b, INMOST_DATA_ENUM_TYPE dim, Storage::real eps)
	{
		for(INMOST_DATA_ENUM_TYPE i = 0; i <  dim; i++)
			if( fabs(a[i]-b[i]) > eps ) return a[i] <= b[i];
		return true;
	}

	class MappingComparator
	{
	public: bool operator () (const std::pair<int,int> & a, const std::pair<int,int> & b) {return a.first < b.first;}
	};
	
	void Mesh::ResolveShared()
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		if( m_state == Mesh::Serial ) SetCommunicator(INMOST_MPI_COMM_WORLD);
		integer dim = GetDimensions();
		int sendsize;
		int mpirank = GetProcessorRank(),mpisize = GetProcessorsNumber();
#if defined(USE_PARALLEL_STORAGE)
		shared_elements.clear();
		ghost_elements.clear();
#endif //USE_PARALLEL_STORAGE
		//determine which bboxes i intersect
		dynarray<int,64> procs;
		Storage::real bbox[6];
		dynarray<Storage::real,384> bboxs(mpisize*6);
		
		for(integer k = 0; k < dim; k++)
		{
			bbox[k] = 1e20;
			bbox[k+dim] = -1e20;
		}
		
		for(iteratorNode it = BeginNode(); it != EndNode(); it++)
		{
			Storage::real_array arr = it->Coords();
			for(integer k = 0; k < dim; k++)
			{
				if( arr[k] < bbox[k] ) bbox[k] = arr[k];
				if( arr[k] > bbox[k+dim] ) bbox[k+dim] = arr[k];
			}
		}
		for(integer k = 0; k < dim; k++)
		{
			REPORT_VAL("min",bbox[k]);
			REPORT_VAL("max",bbox[dim+k]);
		}
		REPORT_MPI(MPI_Allgather(&bbox[0],dim*2,INMOST_MPI_DATA_REAL_TYPE,&bboxs[0],dim*2,INMOST_MPI_DATA_REAL_TYPE,comm));
		for(int k = 0; k < mpisize; k++)
			if( k != mpirank )
			{
				bool flag = true;
				for(integer q = 0; q < dim; q++)
					flag &= !((bbox[q] > bboxs[k*dim*2+q+dim]) || (bbox[dim+q] < bboxs[k*dim*2+q]));
				if( flag ) procs.push_back(k);
			}
		REPORT_VAL("neighbour processors",procs.size());
		
		//~ if( procs.empty() )
		//~ {
			//~ REPORT_STR("no processors around - all elements are owned");
			//~ for(Mesh::iteratorElement it = BeginElement(CELL | EDGE | FACE | NODE); it != EndElement(); it++)
			//~ {
				//~ it->IntegerArrayDV(tag_processors).resize(1);
				//~ it->IntegerDF(tag_owner) = it->IntegerDV(tag_processors) = mpirank;
				//~ it->BulkDF(tag_shared) = Element::Owned;
			//~ }
			//~ ComputeSharedProcs();
			//~ RecomputeParallelStorage(CELL | EDGE | FACE | NODE);
			//~ AssignGlobalID(CELL | EDGE | FACE | NODE);
		//~ }
		//~ else
		{
			bool same_boxes = true, same_box;
			for(int k = 0; k < mpisize && same_boxes; k++)
			{
				same_box = true;
				for(integer j = 0; j < dim*2; j++)
					same_box &= fabs(bbox[j] - bboxs[k*dim*2+j]) < epsilon;
				same_boxes &= same_box;
			}
			
			if( same_boxes )
			{
				REPORT_STR("All bounding boxes are the same - assuming that mesh is replicated over all nodes");
				for(Mesh::iteratorElement it = BeginElement(CELL | EDGE | FACE | NODE); it != EndElement(); it++)
				{
					Storage::integer_array arr = it->IntegerArrayDV(tag_processors);
					arr.resize(mpisize);
					for(int k = 0; k < mpisize; k++) arr[k] = k;
					it->IntegerDF(tag_owner) = 0;
					if( mpirank == 0 ) SetStatus(*it,Element::Shared);
					else SetStatus(*it,Element::Ghost);
				}
				ComputeSharedProcs();
				RecomputeParallelStorage(CELL | EDGE | FACE | NODE);
				AssignGlobalID(CELL | EDGE | FACE | NODE);
#if defined(USE_PARALLEL_STORAGE)
				for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
					for(int i = 0; i < 4; i++)
					{
						if( !it->second[i].empty() )
							std::sort(it->second[i].begin(),it->second[i].end(),GlobalIDComparator(this));
							//qsort(&it->second[i][0],it->second[i].size(),sizeof(Element *),CompareElementsCGID);
					}
				for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)			
					for(int i = 0; i < 4; i++)
					{
						if( !it->second[i].empty() )
							std::sort(it->second[i].begin(),it->second[i].end(),GlobalIDComparator(this));
							//qsort(&it->second[i][0],it->second[i].size(),sizeof(Element *),CompareElementsCGID);
					}
#endif //USE_PARALLEL_STORAGE
			}
			else
			{
				double time = Timer();
				Storage::real epsilon = GetEpsilon();
			
				GeomParam table;
				for(ElementType etype = EDGE; etype <= CELL; etype = etype << 1)
					if( !HaveGeometricData(CENTROID,etype) )
						table[CENTROID] |= etype;
				PrepareGeometricData(table);
				
				time = Timer() - time;
			
				REPORT_STR("Prepare geometric data");
				REPORT_VAL("time",time);
			
				
			
				for(iteratorNode it = BeginNode(); it != EndNode(); it++)
				{
					Storage::integer_array arr = it->IntegerArrayDV(tag_processors);
					arr.resize(1);
					arr[0] = mpirank;
				}
			

				buffer_type exch_data;
				std::vector< INMOST_DATA_REAL_TYPE > unpack_real;
				std::vector< INMOST_DATA_REAL_TYPE > pack_real;
				element_set sorted_nodes;
				
				sorted_nodes.reserve(NumberOfNodes());
				
				
				time = Timer();
				
				
				for(iteratorNode n = BeginNode(); n != EndNode(); n++)
				{
					real_array c = n->Coords();
					for(real_array::size_type k = 0; k < procs.size(); k++)
						if( point_in_bbox(c.data(),bboxs.data()+procs[k]*dim*2,dim) )
						{
							sorted_nodes.push_back(*n);
							break;
						}
				}
				
				time = Timer() - time;
				REPORT_STR("Prepare array of nodes");
				REPORT_VAL("time",time);
				REPORT_VAL("share nodes", sorted_nodes.size());
				REPORT_VAL("total nodes", NumberOfNodes());
				
				
				
				time = Timer();
				if( !sorted_nodes.empty() )
					std::sort(sorted_nodes.begin(),sorted_nodes.end(),CentroidComparator(this));
					//qsort(&sorted_nodes[0],sorted_nodes.size(),sizeof(Element *),CompareElementsCCentroid);
				time = Timer() - time;
				REPORT_STR("Sort nodes");
				REPORT_VAL("time",time);
				
				
				pack_real.reserve(sorted_nodes.size()*dim);
				
				
				time = Timer();
				for(element_set::iterator it = sorted_nodes.begin(); it != sorted_nodes.end(); it++)
				{
					Storage::real_array arr = RealArrayDF(*it,CoordsTag());
					pack_real.insert(pack_real.end(),arr.begin(),arr.end());
				}
				time = Timer() - time;
				REPORT_STR("Gather coordinates");
				REPORT_VAL("time",time);
				
				
				time = Timer();
			
				MPI_Pack_size(static_cast<int>(pack_real.size()),INMOST_MPI_DATA_REAL_TYPE,comm,&sendsize);
				exch_data.resize(sendsize);
				int position = 0;
				if( sendsize > 0 ) MPI_Pack(&pack_real[0],static_cast<INMOST_MPI_SIZE>(pack_real.size()),INMOST_MPI_DATA_REAL_TYPE,&exch_data[0],static_cast<INMOST_MPI_SIZE>(exch_data.size()),&position,comm);
				
				
				time = Timer() - time;
				REPORT_STR("Pack coordinates");
				REPORT_VAL("time",time);

				
				{
					std::vector< MPI_Request > send_reqs, recv_reqs;
					exch_buffer_type send_buffs(procs.size()), recv_buffs(procs.size());
					std::vector<int> done;
//~ #if defined(USE_MPI_P2P)
					//~ unsigned * sendsizeall = shared_space;
					//~ unsigned usend[2] = {sendsize,pack_real.size()};
					//~ REPORT_MPI(MPI_Win_fence(MPI_MODE_NOPRECEDE,window)); //start exchange session
					//~ for(unsigned k = 0; k < procs.size(); k++)
						//~ REPORT_MPI(MPI_Put(usend,2,MPI_UNSIGNED,procs[k],mpirank*2,2,MPI_UNSIGNED,window));
					//~ REPORT_MPI(MPI_Win_fence(MPI_MODE_NOSUCCEED,window)); //end exchange session
//~ #else
					std::vector<unsigned> sendsizeall(mpisize*2);
					int pack_size2 = 0;
					unsigned usend[2] = {static_cast<unsigned>(sendsize),static_cast<unsigned>(pack_real.size())};
					MPI_Pack_size(2,MPI_UNSIGNED,comm,&pack_size2);
					for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
					{
						send_buffs[k].first = procs[k];
						send_buffs[k].second.resize(pack_size2);
						position = 0;
						MPI_Pack(usend,2,MPI_UNSIGNED,&send_buffs[k].second[0],static_cast<INMOST_MPI_SIZE>(send_buffs[k].second.size()),&position,comm);
						recv_buffs[k].first = procs[k];
						recv_buffs[k].second.resize(pack_size2);
					}
					ExchangeBuffersInner(send_buffs,recv_buffs,send_reqs,recv_reqs);
					while( !(done = FinishRequests(recv_reqs)).empty() )
					{
						for(std::vector<int>::iterator qt = done.begin(); qt != done.end(); qt++)
						{
							position = 0;
							MPI_Unpack(&recv_buffs[*qt].second[0],static_cast<INMOST_MPI_SIZE>(recv_buffs[*qt].second.size()),&position,&sendsizeall[procs[*qt]*2],2,MPI_UNSIGNED,comm);
						}
					}
					if( !send_reqs.empty() )
					{
						REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
					}
					//~ REPORT_MPI(MPI_Allgather(usend,2,MPI_UNSIGNED,&sendsizeall[0],2,MPI_UNSIGNED,comm));
//~ #endif
					double time2 = Timer();
					{
						
						
						for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
						{
							send_buffs[k].first = procs[k];
							send_buffs[k].second = exch_data;
							recv_buffs[k].first = procs[k];
							recv_buffs[k].second.resize(sendsizeall[procs[k]*2]);
						}
						
						//PrepareReceiveInner(send_buffs,recv_buffs);
						ExchangeBuffersInner(send_buffs,recv_buffs,send_reqs,recv_reqs);
						
						
						
						while( !(done = FinishRequests(recv_reqs)).empty() )
						{
							for(std::vector<int>::iterator qt = done.begin(); qt != done.end(); qt++)
							{
								time = Timer();
								REPORT_STR("receive node coordinates");
								REPORT_VAL("processor",recv_buffs[*qt].first);
								int count = 0;
								int position = 0;
								unpack_real.resize(sendsizeall[recv_buffs[*qt].first*2+1]);
								MPI_Unpack(&recv_buffs[*qt].second[0],static_cast<INMOST_MPI_SIZE>(recv_buffs[*qt].second.size()),&position,&unpack_real[0],static_cast<INMOST_MPI_SIZE>(unpack_real.size()),INMOST_MPI_DATA_REAL_TYPE,comm);
								std::vector<Storage::real>::iterator it1 = pack_real.begin() , it2 = unpack_real.begin();
								while(it1 != pack_real.end() && it2 != unpack_real.end() )
								{
									int res = 0;
									for(integer k = 0; k < dim; k++)
										if( fabs((*(it1+k))-(*(it2+k))) > epsilon )
										{
											if( (*(it1+k)) < (*(it2+k)) ) res = -1;
											else res = 1;
											break;
										}
									if( res < 0 ) 
										it1 += dim;
									else if( res > 0 ) 
										it2 += dim;
									else
									{
										count++;
										IntegerArrayDV(sorted_nodes[(it1-pack_real.begin())/dim],tag_processors).push_back(recv_buffs[*qt].first);
										it1 += dim;
										it2 += dim;
									}
								}
								REPORT_VAL("intersected coords",count);
								time = Timer() - time;
								REPORT_STR("Intersect coordinates");
								REPORT_VAL("time",time);
							}
						}
						if( !send_reqs.empty() )
						{
							REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
						}
					}			
					
					time2 = Timer() - time2;
					REPORT_STR("Intersect all coordinates");
					REPORT_VAL("time",time2);
					
					time = Timer();
					Element::Status estat;
					for(Mesh::iteratorElement it = BeginElement(NODE); it != EndElement(); it++)
					{
						int owner;
						Storage::integer_array v = it->IntegerArrayDV(tag_processors);
						std::sort(v.begin(),v.end());
						if( v.empty() )
						{
							owner = mpirank;
							v.push_back(mpirank);
						}
						else
							owner = std::min(mpirank,v[0]);
						
						it->IntegerDF(tag_owner) = owner;
						
						if( mpirank == owner )
						{
							if( v.size() == 1 ) 
								estat = Element::Owned;
							else 
								estat = Element::Shared;
						}
						else
							estat = Element::Ghost;
						SetStatus(*it,estat);
					}
					
					ComputeSharedProcs();
					RecomputeParallelStorage(NODE);
					AssignGlobalID(NODE);
					
#if defined(USE_PARALLEL_STORAGE)
					for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
					{
						if( !it->second[0].empty() )
							std::sort(it->second[0].begin(),it->second[0].end(),GlobalIDComparator(this));
							//qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
					}
					for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
					{
						if( !it->second[0].empty() )
							std::sort(it->second[0].begin(),it->second[0].end(),GlobalIDComparator(this));
							//qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
					}
#endif //USE_PARALLEL_STORAGE
					
					
					time = Timer() - time;
					REPORT_STR("Set parallel info for nodes");
					REPORT_VAL("time",time);
					
					
					dynarray<Storage::integer,64> result, intersection;
					
					for(ElementType current_mask = EDGE; current_mask <= CELL; current_mask = current_mask << 1 )
					{
						REPORT_STR("Set parallel info for");
						REPORT_VAL("type",ElementTypeName(current_mask));
						
						
            //int owned_elems = 0;
            //int shared_elems = 0;
						int owner;
						Element::Status estat;
						
						time = Timer();
						//Determine what processors potentially share the element
						for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
						{
							determine_my_procs_low(this,*it, result, intersection);
							Storage::integer_array p = it->IntegerArrayDV(tag_processors);
							if( result.empty() ) 
							{
								p.clear();
								p.push_back(mpirank);
                //++owned_elems;
							}
							else 
              {
                p.replace(p.begin(),p.end(),result.begin(),result.end());
                //if( result.size() == 1 && result[0] == mpirank ) 
                //  ++owned_elems;
                //else ++shared_elems;
              }
						}
						time = Timer() - time;
            //REPORT_VAL("predicted owned elements",owned_elems);
            //REPORT_VAL("predicted shared elements",shared_elems);
						REPORT_STR("Predict processors for elements");
						REPORT_VAL("time",time);
						
						
						time = Timer();
						//Initialize mapping that helps get local id by global id
						std::vector<std::pair<int,int> > mapping;
            REPORT_VAL("mapping type",ElementTypeName(current_mask >> 1));
						for(Mesh::iteratorElement it = BeginElement(current_mask >> 1); it != EndElement(); it++)
						{
							mapping.push_back(std::make_pair(it->GlobalID(),it->LocalID()));
						}
						if( !mapping.empty() ) 
							std::sort(mapping.begin(),mapping.end(),MappingComparator());
						time = Timer() - time;
            REPORT_VAL("mapping size",mapping.size())
						REPORT_STR("Compute global to local indexes mapping");
						REPORT_VAL("time",time);
						//Initialize arrays
						
						time = Timer();
						Storage::integer_array procs = IntegerArrayDV(GetHandle(),tag_processors);
						Storage::integer_array::iterator p = procs.begin();
						std::vector<int> message_send;
						std::vector< std::vector<int> > message_recv(procs.size());
						std::vector< element_set > elements(procs.size());
						
						std::vector< MPI_Request > send_reqs,recv_reqs;
						exch_buffer_type send_buffs(procs.size()), recv_buffs(procs.size());
						
						//Gather all possible shared elements and send global ids of their connectivity to the neighbouring proccessors
						for(p = procs.begin(); p != procs.end(); p++)
						{
							int m = static_cast<int>(p-procs.begin());
							{
								message_send.clear();
								message_send.push_back(0);
								message_send.push_back(0);
								for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
								{
									Storage::integer_array pr = it->IntegerArrayDV(tag_processors);
									if( std::binary_search(pr.begin(),pr.end(),*p) )
									{
										Element::adj_type & sub = LowConn(*it);
										if( sub.size() == 0 ) throw Impossible;
										message_send.push_back(static_cast<int>(sub.size()));
										for(Element::adj_type::iterator kt = sub.begin(); kt != sub.end(); kt++)
											message_send.push_back(GlobalID(*kt));
										message_send[1]++;
										elements[m].push_back(*it);
									}
								}
                REPORT_VAL("for processor",*p);
                REPORT_VAL("gathered elements",elements[m].size());
								message_send[0] = static_cast<int>(message_send.size());
								MPI_Pack_size(static_cast<INMOST_MPI_SIZE>(message_send.size()),MPI_INT,comm,&sendsize);
								send_buffs[m].first = *p;
								send_buffs[m].second.resize(sendsize);
								int position = 0;
								MPI_Pack(&message_send[0],static_cast<INMOST_MPI_SIZE>(message_send.size()),MPI_INT,&send_buffs[m].second[0],static_cast<INMOST_MPI_SIZE>(send_buffs[m].second.size()),&position,comm);
								send_buffs[m].second.resize(position);
								recv_buffs[m].first = *p;
							}
						}
						
						PrepareReceiveInner(UnknownSize,send_buffs,recv_buffs);
						ExchangeBuffersInner(send_buffs,recv_buffs,send_reqs,recv_reqs);
						
						time = Timer() - time;
						REPORT_STR("Pack messages for other processors");
						REPORT_VAL("time",time);
						
						
						time = Timer();
						std::vector<int> done;
						
						while( !(done = FinishRequests(recv_reqs)).empty() )
						{
							for(std::vector<int>::iterator qt = done.begin(); qt != done.end(); qt++)
							{
								int position = 0;
								int size;
								int pos = -1;
								for(p = procs.begin(); p != procs.end(); p++)
									if( *p == recv_buffs[*qt].first )
									{
										pos = static_cast<int>(p - procs.begin());
										break;
									}
								if( pos == -1 ) throw Impossible;
								MPI_Unpack(&recv_buffs[*qt].second[0],static_cast<INMOST_MPI_SIZE>(recv_buffs[*qt].second.size()),&position,&size,1,MPI_INT,comm);
                REPORT_VAL("unpacked message size",size-1);
								message_recv[pos].resize(size-1);
								MPI_Unpack(&recv_buffs[*qt].second[0],static_cast<INMOST_MPI_SIZE>(recv_buffs[*qt].second.size()),&position,&message_recv[pos][0],static_cast<INMOST_MPI_SIZE>(message_recv[pos].size()),MPI_INT,comm);
							}
						}
						
						time = Timer() - time;
						REPORT_STR("Exchange messages");
						REPORT_VAL("time",time);
						
						time = Timer();
						//Now find the difference of local elements with given processor number and remote elements
						for(p = procs.begin(); p != procs.end(); p++)
						{
							int m = static_cast<int>(p-procs.begin());
							element_set remote_elements;
							int pos = 0;
							if( message_recv[m].empty() ) continue;
							int num_remote_elements = message_recv[m][pos++];
							for(int i = 0; i < num_remote_elements; i++)
							{
								int conn_size = message_recv[m][pos++], flag = 1;
								
								dynarray<HandleType,64> sub_elements;
								for(int j = 0; j < conn_size; j++)
								{
									int global_id = message_recv[m][pos++];
									int find = -1;
									std::vector<std::pair<int,int> >::iterator it = std::lower_bound(mapping.begin(),mapping.end(),std::make_pair(global_id,0),MappingComparator());
									if( it != mapping.end() && it->first == global_id) 
										find = static_cast<int>(it-mapping.begin());
									if( find == -1 ) 
									{
										flag = 0;
										pos += conn_size-j-1;
										break;
									}
									int find_local_id = mapping[find].second;
									sub_elements.push_back(ComposeHandle(PrevElementType(current_mask), find_local_id));
									
								}
								if( flag )
								{
									HandleType e = FindSharedAdjacency(sub_elements.data(),static_cast<enumerator>(sub_elements.size()));
									if( e == InvalidHandle() ) continue;
									remote_elements.push_back(e);
								}
							}
              REPORT_VAL("number of unpacked remote elements",remote_elements.size());
              if( !remote_elements.empty() )
              {
                REPORT_VAL("first",remote_elements.front());
                REPORT_VAL("first type",ElementTypeName(GetHandleElementType(remote_elements.front())));
                REPORT_VAL("last",remote_elements.back());
                REPORT_VAL("last type",ElementTypeName(GetHandleElementType(remote_elements.back())));
              }
							std::sort(remote_elements.begin(),remote_elements.end());
              REPORT_VAL("original elements size",elements[m].size());
              if( !elements[m].empty() )
              {
                REPORT_VAL("first",elements[m].front());
                REPORT_VAL("first type",ElementTypeName(GetHandleElementType(elements[m].front())));
                REPORT_VAL("last",elements[m].back());
                REPORT_VAL("last type",ElementTypeName(GetHandleElementType(elements[m].back())));
              }
							std::sort(elements[m].begin(),elements[m].end());
							element_set result;
							element_set::iterator set_end;
							
							result.resize(elements[m].size());
							set_end = std::set_difference(elements[m].begin(),elements[m].end(),remote_elements.begin(),remote_elements.end(), result.begin());
							result.resize(set_end-result.begin());

              REPORT_VAL("set difference size",result.size());
							
							//elements in result are wrongly marked as ghost
							for(element_set::iterator qt = result.begin(); qt != result.end(); qt++)
							{
								integer_array pr = IntegerArrayDV(*qt,tag_processors);
								integer_array::iterator find = std::lower_bound(pr.begin(),pr.end(),*p);
								pr.erase(find);
							}
						}
						
						time = Timer() - time;
						REPORT_STR("Determine true shared based on remote info");
						REPORT_VAL("time",time);
						
						time = Timer();
						//Now mark all the processors status
						for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
						{
							Storage::integer_array pr = it->IntegerArrayDV(tag_processors);
							if( pr.empty() )
							{
								owner = mpirank;
								pr.push_back(mpirank);
							}
							else
								owner = std::min(mpirank,pr[0]);
							it->IntegerDF(tag_owner) = owner;
							if( mpirank == owner )
							{
								if( pr.size() == 1 )
									estat = Element::Owned;
								else
									estat = Element::Shared;
							}
							else
								estat = Element::Ghost;
							SetStatus(*it, estat);
						}
						RecomputeParallelStorage(current_mask);
						if( !send_reqs.empty() )
						{
							REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
						}
						AssignGlobalID(current_mask);					
						integer num = ElementNum(current_mask);
#if defined(USE_PARALLEL_STORAGE)
						for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
						{
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
						}
						for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
						{
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));