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

static std::string NameSlash(std::string input)
{
	for(unsigned 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

namespace INMOST
{
    //static int block_recursion = 0;
	/*
    std::string ro()
    {
        int rank = 0;
#ifdef USE_MPI
        MPI_Comm_rank(MPI_COMM_WORLD,&rank);
#endif
        std::stringstream ss;
        for (int i = 0; i < rank; i++)
            ss << "   ";
        return ss.str();
    }
    */
	//////////////////////////////
	/// 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 ) 
		{
			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);
			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 )
		{
			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 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(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[c][0] = -1;
				tag[c][1] = INT_MAX; 
				//ElementArray<Cell> cells = c.NeighbouringCells();
				ElementArray<Cell> cells = c.BridgeAdjacencies2Cell(bridge);
				for(unsigned 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(unsigned 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(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 );
				//if ( !only_new || it->nbAdjElements(NODE | EDGE | FACE | CELL, NewMarker()) != 0)
				{
					int 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(int 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(int 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(int 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(int 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::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();
	}
	
	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::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;
	}
	
	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;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	Storage::real Mesh::AggregateMin(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_MIN,comm));
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	Storage::integer Mesh::AggregateMin(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_MIN,comm));
#else //USE_MPI
		(void) input;
#endif //USE_MPI
		EXIT_FUNC();
		return output;
	}
	
	void Mesh::AggregateMin(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_MIN,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	void Mesh::AggregateMin(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_MIN,comm));
#else//USE_MPI
		(void) input;
		(void) size;
#endif//USE_MPI
		EXIT_FUNC();
	}
	
	
	INMOST_MPI_Comm Mesh::GetCommunicator() const
	{
#if defined(USE_MPI)
		return comm;
#else //USE_MPI
		return 0;
#endif //USE_MPI
	}
	
	int Mesh::GetProcessorRank() const
	{
#if defined(USE_MPI)
		int rank;
		MPI_Comm_rank(comm,&rank);
		return rank;
#else //USE_MPI
		return 0;
#endif //USE_MPI
	}
	
	int Mesh::GetProcessorsNumber() const
	{
#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
	}
	
    INMOST_MPI_Group Mesh::GetGroup() const
    {
        INMOST_MPI_Group ret = INMOST_MPI_GROUP_EMPTY;
#if defined(USE_MPI)
        MPI_Comm_group(GetCommunicator(), &ret);
#endif
        return ret;
    }
	
	void Mesh::SetCommunicator(INMOST_MPI_Comm _comm)
	{
		ENTER_FUNC();
		tag_shared = CreateTag("PROTECTED_STATUS",DATA_BULK,  ESET |CELL | FACE | EDGE | NODE,NONE,1);
		tag_owner = CreateTag("OWNER_PROCESSOR",DATA_INTEGER, ESET | CELL | FACE | EDGE | NODE,NONE,1);
		tag_processors = CreateTag("PROCESSORS_LIST",DATA_INTEGER, ESET | 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, ESET | CELL | FACE | EDGE | NODE, ESET | CELL | FACE | EDGE | NODE);
		tag_an_id = CreateTag("PROTECTED_ID",DATA_INTEGER, CELL | FACE | EDGE | NODE, CELL | FACE | EDGE | NODE);

#if defined(USE_MPI)
		randomizer = Random();
		
		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,comm));
			pmid = parallel_mesh_unique_id+1;
		}
		m_state = Mesh::Parallel;


#if defined(USE_MPI_P2P)