parallel.cpp

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

#include "inmost.h"


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

//using namespace std;

#if defined(USE_PARALLEL_STORAGE)
bool allow_pack_by_gids = true;
#else // USE_PARALLEL_STORAGE
bool allow_pack_by_gids = false;
#endif // USE_PARALLEL_STORAGE


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

__INLINE std::string NameSlash(std::string input)
{
	for(size_t l = input.size(); l > 0; --l)
		if( input[l-1] == '/' || input[l-1] == '\\' )
			return std::string(input.c_str() + l);
	return input;
}

#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 ENTER_BLOCK() { double btime = Timer(); WriteTab(out_time) << "<FUNCTION name=\"" << __FUNCTION__ << ":" << NameSlash(__FILE__) << ":" << __LINE__ << "\" id=\"func" << GetFuncID()++ << "\">" << std::endl; Enter();
#define EXIT_BLOCK() WriteTab(out_time) << "<TIME>" << Timer() - btime << "</TIME>" << std::endl; Exit(); WriteTab(out_time) << "</FUNCTION>" << std::endl;}

#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 // USE_PARALLEL_WRITE_TIME
#define REPORT_MPI(x) x
#define REPORT_STR(x) {}
#define REPORT_VAL(str,x) {}
#define ENTER_FUNC() {}
#define ENTER_BLOCK()
#define EXIT_BLOCK()
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
#endif // USE_PARALLEL_WRITE_TIME


#if defined(__LINUX__) || defined(__linux__) || defined(__APPLE__)
#include <unistd.h>
#define PROCESSID getpid()
#else
#define PROCESSID -1
#endif


#if 1
#define MPI_Pack_call(data,size,type,buf,buf_size,pos,comm)   MPI_Pack(data,size,type,buf,buf_size,pos,comm)
#define MPI_Pack_size_call(size,type,comm,ret)                MPI_Pack_size(size,type,comm,ret)
#define MPI_Unpack_call(buf,buf_size,pos,data,size,type,comm) MPI_Unpack(buf,buf_size,pos,data,size,type,comm)
#define POS_TYPE int
#define SEND_AS MPI_PACKED
#else
#define MPI_Pack_call(data,size,type,buf,buf_size,pos,comm)   MPI_Pack_external("external32",data,size,type,buf,buf_size,pos)
#define MPI_Pack_size_call(size,type,comm,ret)                MPI_Pack_external_size("external32",size,type,ret)
#define MPI_Unpack_call(buf,buf_size,pos,data,size,type,comm) MPI_Unpack_external("external32",buf,buf_size,pos,data,size,type)
#define POS_TYPE MPI_Aint
#define SEND_AS MPI_BYTE
#endif


#if SIZE_MAX == UCHAR_MAX
#define MPI_SIZE_T MPI_UNSIGNED_CHAR
#elif SIZE_MAX == USHRT_MAX
#define MPI_SIZE_T MPI_UNSIGNED_SHORT
#elif SIZE_MAX == UINT_MAX
#define MPI_SIZE_T MPI_UNSIGNED
#elif SIZE_MAX == ULONG_MAX
#define MPI_SIZE_T MPI_UNSIGNED_LONG
#elif SIZE_MAX == ULLONG_MAX
#define MPI_SIZE_T MPI_UNSIGNED_LONG_LONG
#else
#error "Cannot detect size_t type"
#endif


namespace INMOST
{
#if defined(USE_MPI)
	template<typename T> struct MPIType;
	template<> struct MPIType<bool> {static MPI_Datatype Type() {return MPI_C_BOOL;}};
	template<> struct MPIType<char> {static MPI_Datatype Type() {return MPI_CHAR;}};
	template<> struct MPIType<const char> {static MPI_Datatype Type() {return MPI_CHAR;}};
	template<> struct MPIType<int> {static MPI_Datatype Type() {return MPI_INT;}};
	template<> struct MPIType<short> {static MPI_Datatype Type() {return MPI_SHORT;}};
	template<> struct MPIType<long> {static MPI_Datatype Type() {return MPI_LONG;}};
	template<> struct MPIType<long long> {static MPI_Datatype Type() {return MPI_LONG_LONG;}};
	template<> struct MPIType<unsigned char> {static MPI_Datatype Type() {return MPI_UNSIGNED_CHAR;}};
	template<> struct MPIType<unsigned short> {static MPI_Datatype Type() {return MPI_UNSIGNED_SHORT;}};
	template<> struct MPIType<unsigned> {static MPI_Datatype Type() {return MPI_UNSIGNED;}};
	template<> struct MPIType<unsigned long> {static MPI_Datatype Type() {return MPI_UNSIGNED_LONG;}};
	template<> struct MPIType<unsigned long long> {static MPI_Datatype Type() {return MPI_UNSIGNED_LONG_LONG;}};
	template<> struct MPIType<double> {static MPI_Datatype Type() {return MPI_DOUBLE;}};
	template<> struct MPIType<long double> {static MPI_Datatype Type() {return MPI_LONG_DOUBLE;}};
	template<> struct MPIType<float> {static MPI_Datatype Type() {return MPI_FLOAT;}};
	//template<> struct MPIType<size_t> {static MPI_Datatype Type() {return MPI_SIZE_T;}};
	template<typename T>
	void pack_data(Mesh::buffer_type & buf, const T & data, INMOST_MPI_Comm comm)
	{
		int ierr;
		POS_TYPE pack_size = 0, shift = 0;
		size_t write_pos = buf.size();
		ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size); 
		if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
		buf.resize(write_pos+pack_size);
		ierr = MPI_Pack_call((void*)&data,1,MPIType<T>::Type(),&buf[write_pos],pack_size,&shift,comm);
		if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
		buf.resize(write_pos+shift);
	}
	
	template<typename T>
	size_t pack_data_size(INMOST_MPI_Comm comm)
	{
		int ierr;
		POS_TYPE pack_size = 0;//, shift = 0;
		size_t ret = 0;
		ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size); 
		if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
		ret += pack_size;
		return ret;
	}
	
	
	template<typename T>
	void pack_data_array(Mesh::buffer_type & buf, const T * data, size_t N, INMOST_MPI_Comm comm)
	{
		int ierr;
		if( N )
		{
			POS_TYPE pack_size = 0;
			ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size);
			if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
			size_t rec_bytes = static_cast<size_t>(pack_size);
			size_t max_bytes = rec_bytes*N;
			size_t chunk_bytes = std::min(max_bytes,static_cast<size_t>(INT_MAX));
			size_t chunk_size = chunk_bytes / rec_bytes;
			size_t offset = 0;
			while( offset != N )
			{
				size_t write_pos = buf.size();
				size_t chunk = std::min(chunk_size, N-offset);
				ierr = MPI_Pack_size_call(static_cast<POS_TYPE>(chunk),MPIType<T>::Type(),comm,&pack_size); //pack_size is expected to be within INT_MAX
				if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
				buf.resize(write_pos+pack_size);
				POS_TYPE shift = 0;
				ierr = MPI_Pack_call((void *)&data[offset],static_cast<POS_TYPE>(chunk),MPIType<T>::Type(),&buf[write_pos],pack_size,&shift,comm);
				if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
				buf.resize(write_pos+shift);
				offset += chunk;
			}
		}
	}
	
	
	template<typename T>
	void pack_data_vector(Mesh::buffer_type & buf, const std::vector<T> & data, INMOST_MPI_Comm comm)
	{
		//int ierr;
		pack_data(buf,data.size(),comm);
		if( !data.empty() ) 
			pack_data_array(buf,&data[0],data.size(),comm);
	}
	
	template<typename T>
	size_t pack_data_array_size(size_t N, INMOST_MPI_Comm comm)
	{
		int ierr;
		size_t ret = 0;
		if( N )
		{
			POS_TYPE pack_size = 0;
			ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size);
			if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
			size_t rec_bytes = static_cast<size_t>(pack_size);
			size_t max_bytes = rec_bytes*N;
			size_t chunk_bytes = std::min(max_bytes,static_cast<size_t>(INT_MAX));
			size_t chunk_size = chunk_bytes / rec_bytes;
			size_t offset = 0;
			while( offset != N )
			{
				size_t chunk = std::min(chunk_size, N-offset);
				ierr = MPI_Pack_size_call(static_cast<POS_TYPE>(chunk),MPIType<T>::Type(),comm,&pack_size); //pack_size is expected to be within INT_MAX
				if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
				ret += pack_size;
				offset += chunk;
			}
		}
		return ret;
	}
	
	template<typename T>
	size_t pack_data_vector_size(size_t N, INMOST_MPI_Comm comm)
	{
		size_t ret = pack_data_size<size_t>(comm);
		if( N ) ret += pack_data_array_size<T>(N,comm);
		return ret;
	}
	
	
	template<typename T>
	void unpack_data(Mesh::buffer_type & buf, size_t & buf_pos, T & data, INMOST_MPI_Comm comm)
	{
		int ierr;
		POS_TYPE shift = 0, pack_size;
		ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size); 
		if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
		ierr = MPI_Unpack_call((void*)&buf[buf_pos],pack_size,&shift,&data,1,MPIType<T>::Type(),comm);
		if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
		buf_pos += shift;
	}
	
	
	template<typename T>
	void unpack_data_array(Mesh::buffer_type & buf, size_t & buf_pos, T * data, size_t N, INMOST_MPI_Comm comm)
	{
		int ierr;
		if( N )
		{
			POS_TYPE pack_size = 0;
			ierr = MPI_Pack_size_call(1,MPIType<T>::Type(),comm,&pack_size);
			if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
			size_t rec_bytes = static_cast<size_t>(pack_size);
			size_t max_bytes = rec_bytes*N;
			size_t chunk_bytes = std::min(max_bytes,static_cast<size_t>(INT_MAX));
			size_t chunk_size = chunk_bytes / rec_bytes;
			size_t offset = 0;
			while( offset != N )
			{
				size_t chunk = std::min(chunk_size, N-offset);
				ierr = MPI_Pack_size_call(static_cast<POS_TYPE>(chunk),MPIType<T>::Type(),comm,&pack_size); //pack_size is expected to be within INT_MAX
				if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
				POS_TYPE shift = 0;
				ierr = MPI_Unpack_call((void *)&buf[buf_pos],pack_size,&shift,&data[offset],static_cast<POS_TYPE>(chunk),MPIType<T>::Type(),comm);
				if( ierr != MPI_SUCCESS ) MPI_Abort(comm,__LINE__);
				buf_pos += shift;
				offset += chunk;
			}
		}
	}
	
	template<typename T>
	void unpack_data_vector(Mesh::buffer_type & buf, size_t & buf_pos, std::vector<T> & data, INMOST_MPI_Comm comm)
	{
		//int ierr;
		size_t unpack_size;
		unpack_data(buf,buf_pos,unpack_size,comm);
		data.resize(unpack_size);
		if( !data.empty() )
			unpack_data_array(buf,buf_pos,&data[0],data.size(),comm);
	}
#endif
    
	//////////////////////////////
	/// 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);
			//for(int k = 0; k < size; ++k)
			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 ) 
		{
			INMOST_DATA_ENUM_TYPE 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(static_cast<size_t>(p1.size())+size);
			result.resize(std::set_union(p1.begin(),p1.end(),p2,p2+size,result.begin())-result.begin());
			p1.replace(p1.begin(),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 )
		{
			INMOST_DATA_ENUM_TYPE 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 Mesh::CheckFaces()
    {
        //std::cout << "Check faces" << endl;

        for(Mesh::iteratorFace it = BeginFace(); it != EndFace(); ++it) 
        {
            std::set<int> set_nodes;
            ElementArray<Node> nodes = it->getNodes();
            bool suc = true;
            for (ElementArray<Node>::iterator node = nodes.begin(); node != nodes.end(); node++)
                if (set_nodes.find(node->LocalID()) != set_nodes.end())
                {
                    suc = false;
                    break;
                }
                else
                {
                    set_nodes.insert(node->LocalID());
                }
            //if (suc) std::cout << "-=== Good face: " << setw(2) << it->LocalID();
            //else std::cout << "-=== Error face: " << setw(2) << it->LocalID();

            //cout << ". Nodes = " << nodes.size() << ": ";
            //for (ElementArray<Node>::iterator node = nodes.begin(); node != nodes.end(); node++)
             //   std::cout << node->LocalID() << " ";
            //std::cout << endl;
        }
    }
*/

	static void OperationMinDistance(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void)size;
		assert(size == 2);
		INMOST_DATA_INTEGER_TYPE * idata = (INMOST_DATA_INTEGER_TYPE *)data;
		Storage::integer_array rdata = element->IntegerArray(tag);
		if (idata[1] < rdata[1])
		{
			rdata[0] = idata[0];
			rdata[1] = idata[1];
		}
		else if(idata[1] == rdata[1] && idata[0] < rdata[0])
		{
			rdata[0] = idata[0];
			rdata[1] = idata[1];
		}
	}
	
	static void OperationMinOwner(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void)size;
		assert(size == 1);
		INMOST_DATA_INTEGER_TYPE idata = *(INMOST_DATA_INTEGER_TYPE *)data;
		element->Integer(tag) = std::min(element->Integer(tag),idata);
	}
	
	void Mesh::EquilibrateGhost()//bool only_new)
	{
		static int counter = 0;
		if( GetProcessorsNumber() == 1 ) return;
		ENTER_FUNC();
#if defined(USE_MPI)
		//std::cout << "before_equiv"+std::to_string(counter)+".pvtk" << std::endl;
		//Save("before_equiv"+std::to_string(counter)+".pvtk");
		ElementType bridge = FACE;
		if( tag_bridge.isValid() ) bridge = ElementType(Integer(GetHandle(),tag_bridge));
		TagIntegerArray tag = CreateTag("TEMPORARY_OWNER_DISTANCE_TAG",DATA_INTEGER,CELL,CELL,2);
		std::queue<HandleType> cells_queue;
		ENTER_BLOCK();
		//TODO: compute graph distance only in new elements!!!!
		// Push nearest to owned cell into queue
		for(integer k = 0; k < CellLastLocalID(); k++) if( isValidCell(k) )
		{
			Cell c = CellByLocalID(k);
			if (!c.Hidden() && c.GetStatus() != Element::Owned)
			{
				assert( c.GetStatus() == Element::Shared || c.GetStatus() == Element::Ghost );
				tag[c][0] = -1;
				tag[c][1] = INT_MAX; 
				//ElementArray<Cell> cells = c.NeighbouringCells();
				ElementArray<Cell> cells = c.BridgeAdjacencies2Cell(bridge);
				for(INMOST_DATA_ENUM_TYPE l = 0; l < cells.size(); l++) if (cells[l].GetStatus() == Element::Owned)
				{
					cells_queue.push(c.GetHandle());
					tag[c][0] = cells[l].Integer(tag_owner);
					tag[c][1] = 1;
					break;
				}
			}
		}
		EXIT_BLOCK();
		int cur_dist = 1; // For assert
		ENTER_BLOCK();
		while (!cells_queue.empty())
		{
			Cell c(this,cells_queue.front());
			cells_queue.pop();
			if (cur_dist > tag[c][1]) assert(0); 
			if (cur_dist < tag[c][1]) cur_dist++;
			//ElementArray<Cell> cells = c.NeighbouringCells();
			ElementArray<Cell> cells = c.BridgeAdjacencies2Cell(bridge);
			for(INMOST_DATA_ENUM_TYPE l = 0; l < cells.size(); l++) if (cells[l].GetStatus() != Element::Owned)
			{
				assert( cells[l].GetStatus() == Element::Shared || cells[l].GetStatus() == Element::Ghost );
				if (tag[cells[l]][1] > tag[c][1] + 1)
				{
					tag[cells[l]][0] = tag[c][0];
					tag[cells[l]][1] = tag[c][1] + 1;
					cells_queue.push(cells[l].GetHandle());
				}
			}
		}
		EXIT_BLOCK();
		
		//static int equiv = 0;
		//std::string file;
		//file = "aequiv"+std::to_string(equiv)+".pvtk";
		//std::cout << "Save " << file << std::endl;
		//Save(file);
		
		//CheckCentroids(__FILE__,__LINE__);
		ReduceData(tag, CELL, 0, OperationMinDistance);
		ExchangeData(tag, CELL, 0);
		
		//std::cout << "comput_equiv"+std::to_string(counter)+".pvtk" << std::endl;
		//Save("comput_equiv"+std::to_string(counter)+".pvtk");
		//file = "bequiv"+std::to_string(equiv)+".pvtk";
		//std::cout << "Save " << file << std::endl;
		//Save(file);
		//equiv++;

		CheckGhostSharedCount(__FILE__,__LINE__);
		CheckCentroids(__FILE__,__LINE__);
		
		ENTER_BLOCK();
		//TagInteger tag_new_owner = CreateTag("TEMPORARY_NEW_OWNER",DATA_INTEGER,CELL,NONE,1);
		//for(int k = 0; k < CellLastLocalID(); k++) if( isValidCell(k) )
		//{
		//	Cell c = CellByLocalID(k);
		//	if( !c.Hidden() && c.GetStatus() != Element::Owned )
		//	{
		//		assert( c.GetStatus() == Element::Shared || c.GetStatus() == Element::Ghost );
		//		tag_new_owner[c] = tag[c][0];
		//	}
		//}
		//ExchangeData(tag_new_owner,CELL,0);

		for(integer k = 0; k < CellLastLocalID(); k++) if( isValidCell(k) )
		{
			Cell c = CellByLocalID(k);
			if( !c.Hidden() && c.GetStatus() != Element::Owned )
			{
				assert( c.GetStatus() == Element::Shared || c.GetStatus() == Element::Ghost );
				//if ( !only_new || it->nbAdjElements(NODE | EDGE | FACE | CELL, NewMarker()) != 0)
				{
					integer new_owner = tag[c][0];
					if( new_owner == -1 ) continue;
					//int new_owner = tag_new_owner[c];

					Storage::integer_array procs = c.IntegerArray(ProcessorsTag());
					if( !std::binary_search(procs.begin(),procs.end(),tag[c][0]) ) 
					{
						std::cout << "new owner " << tag[c][0] << " dist " << tag[c][1] << " is not among processors ";
						for(INMOST_DATA_ENUM_TYPE k = 0; k < procs.size();++k) std::cout << procs[k] << " ";
						std::cout << " for cell " << c.LocalID() << " " << c.GlobalID() << std::endl;
					}
					
					c.IntegerDF(tag_owner) = new_owner;
					if (GetProcessorRank() == new_owner)
						c.SetStatus(Element::Shared);
					else
						c.SetStatus(Element::Ghost);
				}
			}
		}
		//DeleteTag(tag_new_owner);
		//ExchangeData(tag_owner,CELL,0);

		RecomputeParallelStorage(CELL);
		CheckGhostSharedCount(__FILE__,__LINE__,CELL);
		EXIT_BLOCK();
		DeleteTag(tag);
		
		
		
		ENTER_BLOCK();
		
		TagInteger new_owner = CreateTag("TEMPORARY_NEW_OWNER",DATA_INTEGER,NODE|EDGE|FACE,NODE|EDGE|FACE,1);
		for(ElementType etype = FACE; etype >= NODE; etype = PrevElementType(etype))
		{
			REPORT_VAL("etype ", ElementTypeName(etype));
			ENTER_BLOCK();
			for(integer k = 0; k < LastLocalID(etype); k++) if( isValidElement(etype,k) )
			{
				Element e = ElementByLocalID(etype,k);
				if( !e.Hidden() && e.GetStatus() != Element::Owned )
				{
					assert( e.GetStatus() == Element::Shared || e.GetStatus() == Element::Ghost );
					//if( !only_new || e.nbAdjElements(NODE | EDGE | FACE | CELL,NewMarker()) )
					{
						Element::adj_type & hc = HighConn(e.GetHandle());
						new_owner[e] = INT_MAX;
						for(INMOST_DATA_ENUM_TYPE l = 0; l < hc.size(); ++l)
						{
							if( !GetMarker(hc[l],HideMarker()) )
								new_owner[e] = std::min(new_owner[e],Integer(hc[l],tag_owner));
						}
					}
				}
			}
			EXIT_BLOCK();
			ReduceData(new_owner,etype,0,OperationMinOwner);
			ExchangeData(new_owner,etype);
			ENTER_BLOCK();
			for(integer k = 0; k < LastLocalID(etype); k++) if( isValidElement(etype,k) )
			{
				Element e = ElementByLocalID(etype,k);
				if( !e.Hidden() && e.GetStatus() != Element::Owned )
				{
					assert( e.GetStatus() == Element::Shared || e.GetStatus() == Element::Ghost );
					//if( !only_new || e.nbAdjElements(NODE | EDGE | FACE | CELL,NewMarker()) )
					{
						if( new_owner[e] != INT_MAX &&
						   e.IntegerDF(tag_owner) != new_owner[e] )
						{
							e.IntegerDF(tag_owner) = new_owner[e];
							if (GetProcessorRank() == new_owner[e])
								e.SetStatus(Element::Shared);
							else
								e.SetStatus(Element::Ghost);
						}
					}
				}
			}
			EXIT_BLOCK();
			RecomputeParallelStorage(etype);
			CheckGhostSharedCount(__FILE__,__LINE__,etype);
		}
		DeleteTag(new_owner);
		//ComputeSharedProcs();
		//std::cout << "after_equiv"+std::to_string(counter)+".pvtk" << std::endl;
		//Save("after_equiv"+std::to_string(counter)+".pvtk");
		EXIT_BLOCK();
		//RecomputeParallelStorage(CELL | EDGE | FACE | NODE);
#endif
		EXIT_FUNC();
		counter++;
	}
	
	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///
	//////////////////////////////
	bool Mesh::HaveGlobalID(ElementType type) const
	{
		assert(OneType(type));
		Tag test;
		if( GlobalIDTag().isValid() ) test = GlobalIDTag();
		//else if( HaveTag("GLOBAL_ID") ) test = GetTag("GLOBAL_ID") );
		if( test.isValid() )
		{
			//bool ret = true;
			//for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype) ) if( etype & types )
			//	ret &= test.isDefined(etype);
			//return ret;
			return test.isDefined(type);
		}
		else return false;
	}


	void Mesh::SynchronizeMarker(MarkerType marker, ElementType mask, SyncBitOp op)
	{
		ENTER_FUNC();
#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
		EXIT_FUNC();
	}
	
	ElementType Mesh::SynchronizeElementType(ElementType etype)
	{
		ENTER_FUNC();
		ElementType etypeout = etype;
#if defined(USE_MPI)
		REPORT_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)
	{
		ENTER_FUNC();
		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)
		REPORT_MPI(MPI_Allreduce(&number,&ret,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm));
#endif//USE_MPI
		EXIT_FUNC();
		return ret;
	}
	Storage::integer Mesh::EnumerateSet(const ElementSet & set, const Tag & num_tag, Storage::integer start, bool define_sparse)
	{
		ENTER_FUNC();
		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++;
		REPORT_MPI(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 || it->GetStatus() == Element::Any) && (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)
		REPORT_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
		EXIT_FUNC();
		return ret;
	}
	Storage::integer Mesh::Enumerate(const HandleType * set, enumerator n, const Tag & num_tag, Storage::integer start, bool define_sparse)
	{
		ENTER_FUNC();
		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++;
		REPORT_MPI(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 || GetStatus(*it) == Element::Any) && (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)
		REPORT_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
		EXIT_FUNC();
		return ret;
	}
	Storage::integer Mesh::Enumerate(ElementType mask, Tag num_tag, Storage::integer start, bool define_sparse)
	{
		ENTER_FUNC();
		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++;
		REPORT_MPI(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 || it->GetStatus() == Element::Any) && (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)
		REPORT_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
		EXIT_FUNC();
		return ret;
	}
	
	Storage::real Mesh::Integrate(Storage::real input)
	{
		ENTER_FUNC();
		Storage::real output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::enumerator Mesh::Integrate(Storage::enumerator input)
	{
		ENTER_FUNC();
		Storage::enumerator output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	
	Storage::integer Mesh::Integrate(Storage::integer input)
	{
		ENTER_FUNC();
		Storage::integer output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	void Mesh::Integrate(Storage::real * input, Storage::integer size)
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		REPORT_MPI(MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	void Mesh::Integrate(Storage::integer * input, Storage::integer size)
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		REPORT_MPI(MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	void Mesh::Integrate(Storage::enumerator * input, Storage::integer size)
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		static dynarray<Storage::enumerator,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::enumerator)*size);
		REPORT_MPI(MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	Storage::integer Mesh::ExclusiveSum(Storage::integer input)
	{
		ENTER_FUNC();
		Storage::integer output = 0;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Scan(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm));
		output -= input;
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::enumerator Mesh::ExclusiveSum(Storage::enumerator input)
	{
		ENTER_FUNC();
		Storage::enumerator output = 0;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Scan(&input,&output,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,comm));
		output -= input;
#else//USE_MPI
		(void) input;
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::real Mesh::Integrate(const Tag & t, enumerator entry, ElementType mask)
	{
		ENTER_FUNC();
		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)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
#endif//USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::real Mesh::AggregateMax(Storage::real input)
	{
		ENTER_FUNC();
		Storage::real output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm));
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::integer Mesh::AggregateMax(Storage::integer input)
	{
		ENTER_FUNC();
		Storage::integer output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm));
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::enumerator Mesh::AggregateMax(Storage::enumerator input)
	{
		ENTER_FUNC();
		Storage::enumerator output = input;
#if defined(USE_MPI)
		REPORT_MPI(MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_MAX,comm));
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	void Mesh::AggregateMax(Storage::real * input, Storage::integer size)
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		REPORT_MPI(MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	void Mesh::AggregateMax(Storage::integer * input, Storage::integer size)
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		REPORT_MPI(MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm));
#else//USE_MPI
		(void) input;