mesh_parallel.cpp 143 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1 2 3 4
#ifdef _MSC_VER //kill some warnings
#define _SCL_SECURE_NO_WARNINGS
#endif

Kirill Terekhov's avatar
Kirill Terekhov committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
#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() long 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;
#else
#define REPORT_MPI(x) x
#define REPORT_STR(x)
#define REPORT_VAL(str,x)
#define ENTER_FUNC()
#define EXIT_FUNC() 
#endif

namespace INMOST
{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
	//////////////////////////////
	/// 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 )
		{
			bool flag = true;
			const Storage::integer * recv = static_cast<const Storage::integer *>(static_cast<const void *>(data));
			Storage::integer_array arr = e->IntegerArray(tag);
			for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it++)
				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);
		}
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
146

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
147
	void UnpackSyncMarkerOR(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
Kirill Terekhov's avatar
Kirill Terekhov committed
148
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
149
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
150 151 152
		element->Bulk(tag) |= *data;
	}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
153
	void UnpackSyncMarkerXOR(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
Kirill Terekhov's avatar
Kirill Terekhov committed
154
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
155
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
156 157 158 159
		element->Bulk(tag) ^= *data;
	}

	
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
160
	void UnpackSyncMarkerAND(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
Kirill Terekhov's avatar
Kirill Terekhov committed
161
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
162
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
163 164 165
		element->Bulk(tag) &= *data;
	}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
166 167 168 169 170
	//////////////////////////////
	/// REDUCTION FUNCTIONS END///
	//////////////////////////////


Kirill Terekhov's avatar
Kirill Terekhov committed
171
	void Mesh::SynchronizeMarker(MarkerType marker, ElementType mask, SyncBitOp op)
Kirill Terekhov's avatar
Kirill Terekhov committed
172 173 174 175 176 177
	{
#if defined(USE_MPI)
		if( m_state == Mesh::Parallel )
		{
			Tag t = CreateTag("TEMP_SYNC_MARKER",DATA_BULK,mask,mask,1);

Kirill Terekhov's avatar
Kirill Terekhov committed
178 179 180
			//workaround for old gcc compiler
			const Element::Status SGhost = Element::Ghost;
			const Element::Status SAny = Element::Any;
Kirill Terekhov's avatar
Kirill Terekhov committed
181
			Element::Status Expr = (Element::Shared | ((op != SYNC_BIT_SET) ? SGhost : SAny));
Kirill Terekhov's avatar
Kirill Terekhov committed
182

Kirill Terekhov's avatar
Kirill Terekhov committed
183
			for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
184
				if( it->GetMarker(marker) && (it->GetStatus() & Expr) )
Kirill Terekhov's avatar
Kirill Terekhov committed
185 186 187 188 189 190
					it->Bulk(t) = 1;

			
			
			switch(op)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
191 192
			case SYNC_BIT_SET:
				ExchangeData(t,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
193 194 195 196 197 198 199 200 201
				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:
Kirill Terekhov's avatar
Kirill Terekhov committed
202
				ReduceData(t,mask,0,UnpackSyncMarkerOR);
Kirill Terekhov's avatar
Kirill Terekhov committed
203 204 205 206 207 208 209 210 211
				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:
Kirill Terekhov's avatar
Kirill Terekhov committed
212
				ReduceData(t,mask,0,UnpackSyncMarkerAND);
Kirill Terekhov's avatar
Kirill Terekhov committed
213 214 215 216 217 218 219 220 221
				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:
Kirill Terekhov's avatar
Kirill Terekhov committed
222
				ReduceData(t,mask,0,UnpackSyncMarkerXOR);
Kirill Terekhov's avatar
Kirill Terekhov committed
223 224 225 226 227 228 229 230 231 232 233 234 235
				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);
		}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
236 237 238 239 240
#else//USE_MPI
		(void) marker;
		(void) mask;
		(void) op;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
241 242 243 244 245 246 247
	}

	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());
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
248
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
249 250
		return etypeout;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
251 252 253 254 255 256 257 258
	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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
259
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
260 261
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
262
	Storage::integer Mesh::EnumerateSet(const ElementSet & set, const Tag & num_tag, Storage::integer start, bool define_sparse)
Kirill Terekhov's avatar
Kirill Terekhov committed
263 264 265 266 267
	{
		Storage::integer shift = 0, ret = 0;
		ElementType mask = CELL | FACE | EDGE | NODE;
#if defined(USE_MPI)
		Storage::integer number = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
268 269
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() != Element::Ghost && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
270 271 272
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
273
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
274
		shift += start;
Kirill Terekhov's avatar
Kirill Terekhov committed
275 276
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() == Element::Owned && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
277
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
278 279
		for(ElementSet::iterator it = set.Begin(); it != set.End(); it++)
			if( it->GetStatus() == Element::Shared && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
280
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
281
		ExchangeData(num_tag,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
282 283 284 285
		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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
286
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
287 288
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
289
	Storage::integer Mesh::Enumerate(const HandleType * set, enumerator n, const Tag & num_tag, Storage::integer start, bool define_sparse)
Kirill Terekhov's avatar
Kirill Terekhov committed
290 291 292 293 294
	{
		Storage::integer shift = 0, ret = 0;
		ElementType mask = CELL | FACE | EDGE | NODE;
#if defined(USE_MPI)
		Storage::integer number = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
295 296
		for(const HandleType * it = set; it != set+n; ++it)
			if( GetStatus(*it) != Element::Ghost && (define_sparse || HaveData(*it,num_tag))) number++;
Kirill Terekhov's avatar
Kirill Terekhov committed
297 298
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
299
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
300
		shift += start;
Kirill Terekhov's avatar
Kirill Terekhov committed
301 302 303 304 305
		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);
Kirill Terekhov's avatar
Kirill Terekhov committed
306 307 308 309
		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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
310
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
311 312
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
313
	Storage::integer Mesh::Enumerate(ElementType mask, Tag num_tag, Storage::integer start, bool define_sparse)
Kirill Terekhov's avatar
Kirill Terekhov committed
314 315 316 317 318
	{
		Storage::integer shift = 0, ret = 0;
#if defined(USE_MPI)
		Storage::integer number = 0;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
319
			if( it->GetStatus() != Element::Ghost && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
320 321 322
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
323
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
324 325
		shift += start;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
326
			if( it->GetStatus() != Element::Owned && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
327 328
				it->Integer(num_tag) = shift++;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
329
			if( it->GetStatus() != Element::Shared && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
330
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
331
		ExchangeData(num_tag,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
332 333 334 335
		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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
336
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
337 338 339 340 341 342 343 344
		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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
345 346 347
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
348 349 350 351 352 353 354 355
		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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
356 357 358
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
359 360 361 362 363 364 365 366 367
		return output;
	}
	
	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;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
368 369 370
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
371 372 373
		return output;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
374
	Storage::real Mesh::Integrate(const Tag & t, enumerator entry, ElementType mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
375 376 377
	{
		Storage::real output = 0, input = 0;
		for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
378 379 380 381 382
			if( GetStatus(*it) != Element::Ghost && HaveData(*it,t) ) 
			{
				real_array arr = RealArray(*it,t);
				if( arr.size() > entry ) input += arr[entry];
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
383 384 385
		output = input;
#if defined(USE_MPI)
		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
386
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
387 388 389
		return output;
	}
	
390 391 392 393 394
	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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
395 396 397
#else //USE_MPI
		(void) input;
#endif //USE_MPI
398 399
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
400
	
401 402 403 404 405
	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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
406 407 408
#else //USE_MPI
		(void) input;
#endif //USE_MPI
409 410
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
411 412 413 414 415 416


	INMOST_MPI_Comm Mesh::GetCommunicator()
	{
#if defined(USE_MPI)
		return comm;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
417
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
418
		return 0;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
419
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
420 421 422 423 424 425 426 427
	}
	
	int Mesh::GetProcessorRank()
	{
#if defined(USE_MPI)
		int rank;
		MPI_Comm_rank(comm,&rank);
		return rank;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
428
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
429
		return 0;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
430
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
431 432 433 434 435 436 437 438
	}
	
	int Mesh::GetProcessorsNumber()
	{
#if defined(USE_MPI)
		int size;
		MPI_Comm_size(comm,&size);
		return size;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
439
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
440
		return 1;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
441
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
442 443 444 445 446 447 448 449
	}	
	
	void Mesh::Initialize(int * argc, char *** argv)
	{
#if defined(USE_MPI)
		int test;
		MPI_Initialized(&test);
        if( test == 0 ) MPI_Init(argc,argv);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
450 451 452 453
#else //USE_MPI
		(void) argc;
		(void) argv;
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
454 455 456 457 458 459 460 461
	}
	void Mesh::Finalize()
	{
#if defined(USE_MPI)
		int test = 0;
		MPI_Finalized(&test);
		if( !test )
			MPI_Finalize();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
462
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
	}
	
	
	
	void Mesh::SetCommunicator(INMOST_MPI_Comm _comm)
	{
		ENTER_FUNC();
		tag_shared = CreateTag("PROTECTED_P_SHARED",DATA_BULK,CELL | FACE | EDGE | NODE,NONE,1);
		tag_owner = CreateTag("P_OWNER_PROCESSOR",DATA_INTEGER, CELL | FACE | EDGE | NODE,NONE,1);
		tag_processors = CreateTag("P_PROCESSORS_LIST",DATA_INTEGER, MESH | NODE | EDGE | FACE | CELL,NONE);
		tag_layers = CreateTag("P_LAYERS",DATA_INTEGER,MESH,NONE,1);
		tag_bridge = CreateTag("P_BRIDGE",DATA_INTEGER,MESH,NONE,1);
		tag_sendto = CreateTag("PROTECTED_P_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;


Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
496
#if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
		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));
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
515 516 517 518
#endif //USE_MPI_P2P
#else //USE_MPI
		(void) _comm;
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534
		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;
	}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
535
#endif //USE_PARALLEL_WRITE_TIME
Kirill Terekhov's avatar
Kirill Terekhov committed
536 537 538
	
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
539
	void determine_my_procs_low(Mesh * m, HandleType h, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
Kirill Terekhov's avatar
Kirill Terekhov committed
540
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
541 542 543
		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());
Kirill Terekhov's avatar
Kirill Terekhov committed
544 545 546 547 548
		result.clear();
		result.insert(result.end(),p.begin(),p.end());
		i++;
		while(i != subelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
549 550
			Storage::integer_array q = m->IntegerArrayDV(*i,m->ProcessorsTag());
			intersection.resize(std::max(static_cast<unsigned>(result.size()),static_cast<unsigned>(q.size())));
Kirill Terekhov's avatar
Kirill Terekhov committed
551 552 553 554 555 556 557
			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++;
		}
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
558
	void determine_my_procs_high(Mesh * m, HandleType h, const Tag & procs, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
Kirill Terekhov's avatar
Kirill Terekhov committed
559
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
560 561
		Element::adj_type const & overelements = m->HighConn(h);
		Element::adj_type::const_iterator i = overelements.begin();
Kirill Terekhov's avatar
Kirill Terekhov committed
562 563 564
		result.clear();
		while(i != overelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
565
			Storage::integer_array q = m->IntegerArrayDV(*i,procs);
Kirill Terekhov's avatar
Kirill Terekhov committed
566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
			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;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
676 677 678 679 680

	class MappingComparator
	{
	public: bool operator () (const std::pair<int,int> & a, const std::pair<int,int> & b) {return a.first < b.first;}
	};
Kirill Terekhov's avatar
Kirill Terekhov committed
681 682 683 684 685 686
	
	void Mesh::ResolveShared()
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		if( m_state == Mesh::Serial ) SetCommunicator(INMOST_MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
687
		integer dim = GetDimensions();
Kirill Terekhov's avatar
Kirill Terekhov committed
688 689 690 691 692
		int sendsize;
		int mpirank = GetProcessorRank(),mpisize = GetProcessorsNumber();
#if defined(USE_PARALLEL_STORAGE)
		shared_elements.clear();
		ghost_elements.clear();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
693
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
694 695 696 697 698
		//determine which bboxes i intersect
		dynarray<int,64> procs;
		Storage::real bbox[6];
		dynarray<Storage::real,384> bboxs(mpisize*6);
		
Kirill Terekhov's avatar
Kirill Terekhov committed
699
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
700 701 702 703 704 705 706 707
		{
			bbox[k] = 1e20;
			bbox[k+dim] = -1e20;
		}
		
		for(iteratorNode it = BeginNode(); it != EndNode(); it++)
		{
			Storage::real_array arr = it->Coords();
Kirill Terekhov's avatar
Kirill Terekhov committed
708
			for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
709 710 711 712 713
			{
				if( arr[k] < bbox[k] ) bbox[k] = arr[k];
				if( arr[k] > bbox[k+dim] ) bbox[k+dim] = arr[k];
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
714
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
715 716 717 718 719 720 721 722 723
		{
			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;
Kirill Terekhov's avatar
Kirill Terekhov committed
724
				for(integer q = 0; q < dim; q++)
Kirill Terekhov's avatar
Kirill Terekhov committed
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748
					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;
Kirill Terekhov's avatar
Kirill Terekhov committed
749
				for(integer j = 0; j < dim*2; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
750 751 752 753 754 755 756 757 758 759 760 761 762
					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;
Kirill Terekhov's avatar
Kirill Terekhov committed
763 764
					if( mpirank == 0 ) SetStatus(*it,Element::Shared);
					else SetStatus(*it,Element::Ghost);
Kirill Terekhov's avatar
Kirill Terekhov committed
765 766 767 768 769 770 771 772 773
				}
				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() )
Kirill Terekhov's avatar
Kirill Terekhov committed
774 775
							std::sort(it->second[i].begin(),it->second[i].end(),GlobalIDComparator(this));
							//qsort(&it->second[i][0],it->second[i].size(),sizeof(Element *),CompareElementsCGID);
Kirill Terekhov's avatar
Kirill Terekhov committed
776 777 778 779 780
					}
				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() )
Kirill Terekhov's avatar
Kirill Terekhov committed
781 782
							std::sort(it->second[i].begin(),it->second[i].end(),GlobalIDComparator(this));
							//qsort(&it->second[i][0],it->second[i].size(),sizeof(Element *),CompareElementsCGID);
Kirill Terekhov's avatar
Kirill Terekhov committed
783
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
784
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
785 786 787 788 789 790
			}
			else
			{
				double time = Timer();
				Storage::real epsilon = GetEpsilon();
			
Kirill Terekhov's avatar
Kirill Terekhov committed
791
				GeomParam table;
Kirill Terekhov's avatar
Kirill Terekhov committed
792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811
				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;
				}
			

Kirill Terekhov's avatar
Kirill Terekhov committed
812
				buffer_type exch_data;
Kirill Terekhov's avatar
Kirill Terekhov committed
813 814
				std::vector< INMOST_DATA_REAL_TYPE > unpack_real;
				std::vector< INMOST_DATA_REAL_TYPE > pack_real;
Kirill Terekhov's avatar
Kirill Terekhov committed
815
				element_set sorted_nodes;
Kirill Terekhov's avatar
Kirill Terekhov committed
816 817 818 819 820 821 822 823 824
				
				sorted_nodes.reserve(NumberOfNodes());
				
				
				time = Timer();
				
				
				for(iteratorNode n = BeginNode(); n != EndNode(); n++)
				{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
825 826
					real_array c = n->Coords();
					for(real_array::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
827
						if( point_in_bbox(c.data(),bboxs.data()+procs[k]*dim*2,dim) )
Kirill Terekhov's avatar
Kirill Terekhov committed
828
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
829
							sorted_nodes.push_back(*n);
Kirill Terekhov's avatar
Kirill Terekhov committed
830 831 832 833 834 835 836 837 838 839 840 841 842 843
							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() )
Kirill Terekhov's avatar
Kirill Terekhov committed
844 845
					std::sort(sorted_nodes.begin(),sorted_nodes.end(),CentroidComparator(this));
					//qsort(&sorted_nodes[0],sorted_nodes.size(),sizeof(Element *),CompareElementsCCentroid);
Kirill Terekhov's avatar
Kirill Terekhov committed
846 847 848 849 850 851 852 853 854
				time = Timer() - time;
				REPORT_STR("Sort nodes");
				REPORT_VAL("time",time);
				
				
				pack_real.reserve(sorted_nodes.size()*dim);
				
				
				time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
855
				for(element_set::iterator it = sorted_nodes.begin(); it != sorted_nodes.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
856
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
857
					Storage::real_array arr = RealArrayDF(*it,CoordsTag());
Kirill Terekhov's avatar
Kirill Terekhov committed
858 859 860 861 862 863 864 865 866
					pack_real.insert(pack_real.end(),arr.begin(),arr.end());
				}
				time = Timer() - time;
				REPORT_STR("Gather coordinates");
				REPORT_VAL("time",time);
				
				
				time = Timer();
			
Kirill Terekhov's avatar
Kirill Terekhov committed
867
				MPI_Pack_size(static_cast<int>(pack_real.size()),INMOST_MPI_DATA_REAL_TYPE,comm,&sendsize);
Kirill Terekhov's avatar
Kirill Terekhov committed
868 869
				exch_data.resize(sendsize);
				int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
870
				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);
Kirill Terekhov's avatar
Kirill Terekhov committed
871 872 873 874 875 876 877 878 879 880 881
				
				
				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;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
882
//~ #if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
883 884 885 886 887 888 889 890 891
					//~ 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;
Kirill Terekhov's avatar
Kirill Terekhov committed
892
					unsigned long usend[2] = {static_cast<unsigned>(sendsize),static_cast<unsigned>(pack_real.size())};
Kirill Terekhov's avatar
Kirill Terekhov committed
893
					MPI_Pack_size(2,MPI_UNSIGNED,comm,&pack_size2);
Kirill Terekhov's avatar
Kirill Terekhov committed
894
					for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
895 896 897 898
					{
						send_buffs[k].first = procs[k];
						send_buffs[k].second.resize(pack_size2);
						position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
899
						MPI_Pack(usend,2,MPI_UNSIGNED,&send_buffs[k].second[0],static_cast<INMOST_MPI_SIZE>(send_buffs[k].second.size()),&position,comm);
Kirill Terekhov's avatar
Kirill Terekhov committed
900 901 902 903 904 905 906 907 908
						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;
Kirill Terekhov's avatar
Kirill Terekhov committed
909
							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);
Kirill Terekhov's avatar
Kirill Terekhov committed
910 911 912
						}
					}
					if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
913
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
914
						REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
Kirill Terekhov's avatar
Kirill Terekhov committed
915
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
916 917 918 919 920 921
					//~ REPORT_MPI(MPI_Allgather(usend,2,MPI_UNSIGNED,&sendsizeall[0],2,MPI_UNSIGNED,comm));
//~ #endif
					double time2 = Timer();
					{
						
						
Kirill Terekhov's avatar
Kirill Terekhov committed
922
						for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944
						{
							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]);
Kirill Terekhov's avatar
Kirill Terekhov committed
945
								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);
Kirill Terekhov's avatar
Kirill Terekhov committed
946 947 948 949
								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;
Kirill Terekhov's avatar
Kirill Terekhov committed
950
									for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
951 952 953 954 955 956 957 958 959 960 961 962 963
										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++;
Kirill Terekhov's avatar
Kirill Terekhov committed
964
										IntegerArrayDV(sorted_nodes[(it1-pack_real.begin())/dim],tag_processors).push_back(recv_buffs[*qt].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
965 966 967 968 969 970 971 972 973 974 975
										it1 += dim;
										it2 += dim;
									}
								}
								REPORT_VAL("intersected coords",count);
								time = Timer() - time;
								REPORT_STR("Intersect coordinates");
								REPORT_VAL("time",time);
							}
						}
						if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
976
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
977
							REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
Kirill Terekhov's avatar
Kirill Terekhov committed
978
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
979 980 981 982 983 984 985
					}			
					
					time2 = Timer() - time2;
					REPORT_STR("Intersect all coordinates");
					REPORT_VAL("time",time2);
					
					time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
986
					Element::Status estat;
Kirill Terekhov's avatar
Kirill Terekhov committed
987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003
					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 )
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1004 1005 1006 1007
							if( v.size() == 1 ) 
								estat = Element::Owned;
							else 
								estat = Element::Shared;
Kirill Terekhov's avatar
Kirill Terekhov committed
1008 1009
						}
						else
Kirill Terekhov's avatar
Kirill Terekhov committed
1010 1011
							estat = Element::Ghost;
						SetStatus(*it,estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1012 1013 1014 1015 1016 1017
					}
					
					ComputeSharedProcs();
					RecomputeParallelStorage(NODE);
					AssignGlobalID(NODE);
					
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1018
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1019 1020 1021
					for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
					{
						if( !it->second[0].empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1022 1023
							std::sort(it->second[0].begin(),it->second[0].end(),GlobalIDComparator(this));
							//qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
Kirill Terekhov's avatar
Kirill Terekhov committed
1024 1025 1026 1027
					}
					for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
					{
						if( !it->second[0].empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1028 1029
							std::sort(it->second[0].begin(),it->second[0].end(),GlobalIDComparator(this));
							//qsort(&it->second[0][0],it->second[0].size(),sizeof(Element *),CompareElementsCGID);
Kirill Terekhov's avatar
Kirill Terekhov committed
1030
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1031
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053
					
					
					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 owner;
						Element::Status estat;
						
						time = Timer();
						//Determine what processors potentially share the element
						for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1054
							determine_my_procs_low(this,*it, result, intersection);
Kirill Terekhov's avatar
Kirill Terekhov committed
1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069
							Storage::integer_array p = it->IntegerArrayDV(tag_processors);
							if( result.empty() ) 
							{
								p.clear();
								p.push_back(mpirank);
							}
							else p.replace(p.begin(),p.end(),result.begin(),result.end());
						}
						time = Timer() - time;
						REPORT_STR("Predict processors for elements");
						REPORT_VAL("time",time);
						
						
						time = Timer();
						//Initialize mapping that helps get local id by global id
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1070
						std::vector<std::pair<int,int> > mapping;
Kirill Terekhov's avatar
Kirill Terekhov committed
1071 1072
						for(Mesh::iteratorElement it = BeginElement(current_mask >> 1); it != EndElement(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1073
							mapping.push_back(std::make_pair(it->GlobalID(),it->LocalID()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1074 1075
						}
						if( !mapping.empty() ) 
Kirill Terekhov's avatar
Kirill Terekhov committed
1076
							std::sort(mapping.begin(),mapping.end(),MappingComparator());
Kirill Terekhov's avatar
Kirill Terekhov committed
1077 1078 1079 1080 1081 1082
						time = Timer() - time;
						REPORT_STR("Compute global to local indexes mapping");
						REPORT_VAL("time",time);
						//Initialize arrays
						
						time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
1083
						Storage::integer_array procs = IntegerArrayDV(GetHandle(),tag_processors);
Kirill Terekhov's avatar
Kirill Terekhov committed
1084 1085 1086
						Storage::integer_array::iterator p = procs.begin();
						std::vector<int> message_send;
						std::vector< std::vector<int> > message_recv(procs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1087
						std::vector< element_set > elements(procs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1088 1089 1090 1091 1092 1093 1094
						
						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++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1095
							int m = static_cast<int>(p-procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1096 1097 1098 1099 1100 1101 1102 1103 1104
							{
								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) )
									{
Kirill Terekhov's avatar
Kirill Terekhov committed
1105
										Element::adj_type & sub = LowConn(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1106
										if( sub.size() == 0 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1107 1108 1109
										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));
Kirill Terekhov's avatar
Kirill Terekhov committed
1110
										message_send[1]++;
Kirill Terekhov's avatar
Kirill Terekhov committed
1111
										elements[m].push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1112 1113
									}
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
1114 1115
								message_send[0] = static_cast<int>(message_send.size());
								MPI_Pack_size(static_cast<INMOST_MPI_SIZE>(message_send.size()),MPI_INT,comm,&sendsize);
Kirill Terekhov's avatar
Kirill Terekhov committed
1116 1117 1118
								send_buffs[m].first = *p;
								send_buffs[m].second.resize(sendsize);
								int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1119
								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);
Kirill Terekhov's avatar
Kirill Terekhov committed
1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145
								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 )
									{
Kirill Terekhov's avatar
Kirill Terekhov committed
1146
										pos = static_cast<int>(p - procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1147 1148 1149
										break;
									}
								if( pos == -1 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1150
								MPI_Unpack(&recv_buffs[*qt].second[0],static_cast<INMOST_MPI_SIZE>(recv_buffs[*qt].second.size()),&position,&size,1,MPI_INT,comm);
Kirill Terekhov's avatar
Kirill Terekhov committed
1151
								message_recv[pos].resize(size-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1152
								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);
Kirill Terekhov's avatar
Kirill Terekhov committed
1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163
							}
						}
						
						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++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1164 1165
							int m = static_cast<int>(p-procs.begin());
							element_set remote_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1166 1167 1168 1169 1170 1171 1172
							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;
								
Kirill Terekhov's avatar
Kirill Terekhov committed
1173
								dynarray<HandleType,64> sub_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1174 1175 1176
								for(int j = 0; j < conn_size; j++)
								{
									int global_id = message_recv[m][pos++];
Kirill Terekhov's avatar
Kirill Terekhov committed
1177 1178 1179 1180
									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());
Kirill Terekhov's avatar
Kirill Terekhov committed
1181 1182 1183 1184 1185 1186
									if( find == -1 ) 
									{
										flag = 0;
										pos += conn_size-j-1;
										break;
									}
Kirill Terekhov's avatar
Kirill Terekhov committed
1187 1188
									int find_local_id = mapping[find].second;
									sub_elements.push_back(ComposeHandle(current_mask >> 1, find_local_id));
Kirill Terekhov's avatar
Kirill Terekhov committed
1189 1190 1191 1192
									
								}
								if( flag )
								{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1193
									HandleType e = FindSharedAdjacency(sub_elements.data(),static_cast<enumerator>(sub_elements.size()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1194
									if( e == InvalidHandle() ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
1195 1196 1197 1198 1199
									remote_elements.push_back(e);
								}
							}
							std::sort(remote_elements.begin(),remote_elements.end());
							std::sort(elements[m].begin(),elements[m].end());
Kirill Terekhov's avatar
Kirill Terekhov committed
1200 1201
							element_set result;
							element_set::iterator set_end;
Kirill Terekhov's avatar
Kirill Terekhov committed
1202 1203 1204 1205 1206 1207
							
							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());				
							
							//elements in result are wrongly marked as ghost
Kirill Terekhov's avatar
Kirill Terekhov committed
1208
							for(element_set::iterator qt = result.begin(); qt != result.end(); qt++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1209
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
1210 1211
								integer_array pr = IntegerArrayDV(*qt,tag_processors);
								integer_array::iterator find = std::lower_bound(pr.begin(),pr.end(),*p);
Kirill Terekhov's avatar
Kirill Terekhov committed
1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241
								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;
Kirill Terekhov's avatar
Kirill Terekhov committed
1242
							SetStatus(*it, estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1243 1244 1245
						}
						RecomputeParallelStorage(current_mask);
						if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1246
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1247
							REPORT_MPI(MPI_Waitall(static_cast<INMOST_MPI_SIZE>(send_reqs.size()),&send_reqs[0],MPI_STATUSES_IGNORE));
Kirill Terekhov's avatar
Kirill Terekhov committed
1248
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1249
						AssignGlobalID(current_mask);					
Kirill Terekhov's avatar
Kirill Terekhov committed
1250
						integer num = ElementNum(current_mask);
Kirill Terekhov's avatar
Kirill Terekhov committed
1251 1252 1253
#if defined(USE_PARALLEL_STORAGE)
						for(parallel_storage::iterator it = shared_elements.begin(); it != shared_elements.end(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1254 1255
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1256 1257 1258
						}
						for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1259 1260
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1261
						}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1262
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273
	
						time = Timer() - time;
						REPORT_STR("Set parallel info");
						REPORT_VAL("time",time);
					}
				}
				RemoveGeometricData(table);	
			}
		}
		//RemoveGhost();
		
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1274
#else //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1275
		AssignGlobalID(CELL | FACE | EDGE | NODE);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1276
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
1277 1278 1279 1280 1281 1282 1283 1284 1285
		EXIT_FUNC();
	}
	
	
	void Mesh::RemoveGhost()
	{
		if( m_state == Mesh::Serial ) return;
		ENTER_FUNC()
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
1286
		element_set delete_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1287
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1288
		proc_elements del_shared, del_ghost;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1289
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1290 1291 1292
		int mpirank = GetProcessorRank();
		Tag tag_delete = CreateTag("TEMPORARY_DELETE_GHOST_TAG",DATA_INTEGER,FACE | EDGE | NODE,FACE | EDGE | NODE);
		
Kirill Terekhov's avatar
Kirill Terekhov committed
1293
		for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1294
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1295
			Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1296 1297 1298
			if( estat == Element::Ghost ) 
			{
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1299
				del_ghost[it->IntegerDF(tag_owner)].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1300
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1301 1302 1303 1304
				delete_elements.push_back(*it);
			}
			else if( estat == Element::Shared )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1305 1306
				SetStatus(*it,Element::Owned);
				Storage::integer_array it_procs = it->IntegerArrayDV(tag_processors);
Kirill Terekhov's avatar
Kirill Terekhov committed
1307 1308 1309 1310 1311
				
#if defined(USE_PARALLEL_STORAGE)
				for(Storage::integer_array::iterator vit = it_procs.begin(); vit != it_procs.end(); vit++)
					if( *vit != mpirank )
						del_shared[*vit].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1312
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1313 1314 1315 1316 1317 1318 1319
				
				it_procs.resize(1);
				it_procs[0] = mpirank;
			}
		}
		
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1320
		for(proc_elements::iterator it = del_ghost.begin(); it != del_ghost.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1321
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335
			element_set & ref = ghost_elements[it->first][ElementNum(CELL)];
			if( !it->second.empty() ) 
			{
				if( have_global_id & CELL )
					std::sort(it->second.begin(),it->second.end(),GlobalIDComparator(this));
				else
					std::sort(it->second.begin(),it->second.end(),CentroidComparator(this));
			}
			element_set result(ref.size());
			element_set::iterator end;
			if( have_global_id & CELL )
				end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),GlobalIDComparator(this));
			else
				end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1336 1337 1338 1339
			result.resize(end-result.begin());
			ref.swap(result);
		}
		del_ghost.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1340
		for(proc_elements::iterator it = del_shared.begin(); it != del_shared.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1341
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355
			element_set & ref = shared_elements[it->first][ElementNum(CELL)];
			if( !it->second.empty() ) 
			{
				if( have_global_id & CELL )
					std::sort(it->second.begin(),it->second.end(),GlobalIDComparator(this));
				else
					std::sort(it->second.begin(),it->second.end(),CentroidComparator(this));
			}
			element_set result(ref.size());
			element_set::iterator end;
			if( have_global_id & CELL )
				end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),GlobalIDComparator(this));
			else
				end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1356 1357 1358 1359
			result.resize(end-result.begin());
			ref.swap(result);
		}
		del_shared.clear();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1360
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1361

Kirill Terekhov's avatar
Kirill Terekhov committed
1362
		for(element_set::iterator it = delete_elements.begin(); it != delete_elements.end(); it++) Destroy(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1363
		
Kirill Terekhov's avatar
Kirill Terekhov committed
1364
		delete_elements.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1365 1366 1367 1368
		for(ElementType mask = FACE; mask >= NODE; mask = mask >> 1)
		{
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1369
				Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1370 1371 1372
				if( estat == Element::Owned || estat == Element::Shared ) continue;
				if( it->nbAdjElements(mask<<1) == 0 ) it->Integer(tag_delete) = mpirank;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1373 1374
			ReduceData(tag_delete,mask,0,DeleteUnpack);
			ExchangeData(tag_delete,0,mask);
Kirill Terekhov's avatar
Kirill Terekhov committed
1375 1376
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1377
				Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1378 1379 1380 1381
				if( estat == Element::Owned ) continue;
				if( estat == Element::Ghost && it->nbAdjElements(mask<<1) == 0 ) 
				{
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1382
					del_ghost[it->IntegerDF(tag_owner)].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1383
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1384
					delete_elements.push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1385 1386 1387 1388 1389 1390 1391 1392 1393
				}
				else if( it->HaveData(tag_delete) )
				{
					Storage::integer_array procs = it->IntegerArrayDV(tag_processors);
					Storage::integer_array del_procs = it->IntegerArray(tag_delete);
					std::vector<Storage::integer> result(procs.size());
					std::sort(del_procs.begin(),del_procs.end());								
#if defined(USE_PARALLEL_STORAGE)				
					for(Storage::integer_array::iterator vit = del_procs.begin(); vit != del_procs.end(); vit++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1394
						del_shared[*vit].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1395
#endif //USE_PARALLEL_STORAGE				
Kirill Terekhov's avatar
Kirill Terekhov committed
1396 1397 1398 1399 1400 1401
					std::vector<Storage::integer>::iterator end = std::set_difference(procs.begin(),procs.end(),del_procs.begin(),del_procs.end(),result.begin());
					result.resize(end-result.begin());
					procs.clear();
					procs.insert(procs.begin(),result.begin(),result.end());
					
					if( procs.size() == 1 && procs[0] == mpirank )
Kirill Terekhov's avatar
Kirill Terekhov committed
1402
						SetStatus(*it,Element::Owned);
Kirill Terekhov's avatar
Kirill Terekhov committed
1403 1404 1405 1406
				}
			}
			
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1407
			for(proc_elements::iterator it = del_ghost.begin(); it != del_ghost.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1408
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422
				element_set & ref = ghost_elements[it->first][ElementNum(mask)];
				if( !it->second.empty() ) 
				{
					if( have_global_id & mask )
						std::sort(it->second.begin(),it->second.end(),GlobalIDComparator(this));
					else
						std::sort(it->second.begin(),it->second.end(),CentroidComparator(this));
				}
				element_set result(ref.size());
				element_set::iterator end;
				if( have_global_id & CELL )
					end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),GlobalIDComparator(this));
				else
					end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1423 1424 1425 1426
				result.resize(end-result.begin());
				ref.swap(result);
			}
			del_ghost.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1427
			for(proc_elements::iterator it = del_shared.begin(); it != del_shared.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1428
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442
				element_set & ref = shared_elements[it->first][ElementNum(mask)];
				if( !it->second.empty() ) 
				{
					if( have_global_id & mask )
						std::sort(it->second.begin(),it->second.end(),GlobalIDComparator(this));
					else
						std::sort(it->second.begin(),it->second.end(),CentroidComparator(this));
				}
				element_set result(ref.size());
				element_set::iterator end;
				if( have_global_id & CELL )
					end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),GlobalIDComparator(this));
				else
					end = std::set_difference(ref.begin(),ref.end(),it->second.begin(),it->second.end(),result.begin(),CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1443 1444 1445 1446
				result.resize(end-result.begin());
				ref.swap(result);
			}
			del_shared.clear();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1447
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1448
			for(element_set::iterator it = delete_elements.begin(); it != delete_elements.end(); it++) Destroy(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1449 1450 1451 1452
			delete_elements.clear();
		}		
		DeleteTag(tag_delete);
		ComputeSharedProcs();
Kirill Terekhov's avatar
Kirill Terekhov committed
1453 1454
		Integer(GetHandle(),tag_layers) = 0;
		Integer(GetHandle(),tag_bridge) = NONE;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1455
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
1456 1457 1458
		EXIT_FUNC();
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1459
	void Mesh::RemoveGhostElements(const HandleType * ghost, enumerator num)
Kirill Terekhov's avatar
Kirill Terekhov committed
1460 1461 1462 1463 1464
	{
		ENTER_FUNC();
#if defined(USE_MPI)
		int mpirank = GetProcessorRank();
		Tag tag_delete = CreateTag("TEMPORARY_DELETE_GHOST_ELEMENTS_TAG",DATA_INTEGER,CELL | FACE | EDGE | NODE,CELL | FACE | EDGE | NODE);
Kirill Terekhov's avatar
Kirill Terekhov committed
1465
		element_set delete_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1466
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1467
		proc_elements del_shared, del_ghost;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1468
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1469 1470 1471 1472 1473
		
		for(ElementType mask = CELL; mask >= NODE; mask = mask >> 1)
		{
			if( mask & CELL )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1474 1475
				for(const HandleType * it = ghost; it != ghost + num; it++)
					if( GetStatus(*it) == Element::Ghost ) Integer(*it,tag_delete) = mpirank;
Kirill Terekhov's avatar
Kirill Terekhov committed
1476 1477 1478 1479 1480
			}
			else 
			{
				for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1481
					Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1482 1483 1484 1485
					if( estat == Element::Owned || estat == Element::Shared ) continue;
					if( it->nbAdjElements(mask<<1) == 0 ) it->Integer(tag_delete) = mpirank;
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1486 1487
			ReduceData(tag_delete,mask,0,DeleteUnpack);
			ExchangeData(tag_delete,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
1488 1489
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1490
				Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1491 1492 1493 1494 1495 1496 1497 1498 1499
				if( estat == Element::Owned ) continue;
				if( it->HaveData(tag_delete) )
				{
						Storage::integer_array del_procs = it->IntegerArray(tag_delete);
						std::sort(del_procs.begin(),del_procs.end());
						
						if( estat == Element::Ghost && std::binary_search(del_procs.begin(),del_procs.end(),mpirank) )
						{
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1500
							del_ghost[it->IntegerDF(tag_owner)].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1501
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1502
							delete_elements.push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1503 1504 1505 1506 1507 1508 1509 1510 1511
						}
						else
						{
							Storage::integer_array procs = it->IntegerArrayDV(tag_processors);
							std::vector<Storage::integer> result(procs.size());
#if defined(USE_PARALLEL_STORAGE)
							if( estat == Element::Shared )
							{
								for(Storage::integer_array::iterator vit = del_procs.begin(); vit != del_procs.end(); vit++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1512
									del_shared[*vit].push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1513
							}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1514
#endif	//USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1515 1516 1517 1518 1519 1520
							std::vector<Storage::integer>::iterator end = std::set_difference(procs.begin(),procs.end(),del_procs.begin(),del_procs.end(),result.begin());
							result.resize(end-result.begin());
							procs.clear();
							procs.insert(procs.begin(),result.begin(),result.end());
							
							if( procs.size() == 1 && procs[0] == mpirank )