mesh_parallel.cpp 157 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
#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;}
Kirill Terekhov's avatar
Kirill Terekhov committed
22 23 24
#define ENTER_FUNC() double all_time = Timer(); {WriteTab(out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << func_id++ << "\">" << std::endl; Enter();}
#define EXIT_FUNC() {WriteTab(out_time) << "<TIME>" << Timer() - all_time << "</TIME>" << std::endl; Exit(); WriteTab(out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC_DIE() {WriteTab(out_time) << "<TIME>" << -1 << "</TIME>" << std::endl; Exit(); WriteTab(out_time) << "</FUNCTION>" << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
25 26
#else
#define REPORT_MPI(x) x
Kirill Terekhov's avatar
Kirill Terekhov committed
27 28 29 30 31
#define REPORT_STR(x) {}
#define REPORT_VAL(str,x) {}
#define ENTER_FUNC() {}
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
Kirill Terekhov's avatar
Kirill Terekhov committed
32 33 34 35
#endif

namespace INMOST
{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
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
	//////////////////////////////
	/// 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);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
93
			for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it+=2)
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
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 146 147
				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
148

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
149
	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
150
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
151
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
152 153 154
		element->Bulk(tag) |= *data;
	}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
155
	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
156
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
157
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
158 159 160 161
		element->Bulk(tag) ^= *data;
	}

	
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
162
	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
163
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
164
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
165 166 167
		element->Bulk(tag) &= *data;
	}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
168 169 170 171 172
	//////////////////////////////
	/// REDUCTION FUNCTIONS END///
	//////////////////////////////


Kirill Terekhov's avatar
Kirill Terekhov committed
173

Kirill Terekhov's avatar
Kirill Terekhov committed
174
	void Mesh::SynchronizeMarker(MarkerType marker, ElementType mask, SyncBitOp op)
Kirill Terekhov's avatar
Kirill Terekhov committed
175 176 177 178 179 180
	{
#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
181 182 183
			//workaround for old gcc compiler
			const Element::Status SGhost = Element::Ghost;
			const Element::Status SAny = Element::Any;
Kirill Terekhov's avatar
Kirill Terekhov committed
184
			Element::Status Expr = (Element::Shared | ((op != SYNC_BIT_SET) ? SGhost : SAny));
Kirill Terekhov's avatar
Kirill Terekhov committed
185

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

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

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


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

#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
502
#if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520
		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
521 522 523 524
#endif //USE_MPI_P2P
#else //USE_MPI
		(void) _comm;
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
		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
Kirill Terekhov committed
541 542 543 544 545 546 547 548 549 550 551 552
	void Mesh::FinalizeFile()
	{
		if( tab > 1 )
		{
			REPORT_STR("Finalizing file with nonempty function stack - probably due to error.");
		}
		while(tab > 1)
		{
			EXIT_FUNC_DIE();
		}
		out_time << "</Debug>" << std::endl;
	}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
553
#endif //USE_PARALLEL_WRITE_TIME
Kirill Terekhov's avatar
Kirill Terekhov committed
554 555 556
	
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
557
	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
558
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
559 560 561
		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
562 563 564 565 566
		result.clear();
		result.insert(result.end(),p.begin(),p.end());
		i++;
		while(i != subelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
567 568
			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
569 570 571 572 573 574 575
			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
576
	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
577
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
578
		Element::adj_type const & overelements = m->HighConn(h);
Kirill Terekhov's avatar
Kirill Terekhov committed
579
		if( overelements.empty() ) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
580
		Element::adj_type::const_iterator i = overelements.begin();
Kirill Terekhov's avatar
Kirill Terekhov committed
581 582 583
		result.clear();
		while(i != overelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
584
			Storage::integer_array q = m->IntegerArrayDV(*i,procs);
Kirill Terekhov's avatar
Kirill Terekhov committed
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 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694
			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
695 696 697 698 699

	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
700 701 702 703 704 705
	
	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
706
		integer dim = GetDimensions();
Kirill Terekhov's avatar
Kirill Terekhov committed
707 708 709 710 711
		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
712
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
713 714 715 716 717
		//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
718
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
719 720 721 722 723 724 725 726
		{
			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
727
			for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
728 729 730 731 732
			{
				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
733
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
734 735 736 737 738 739 740 741 742
		{
			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
743
				for(integer q = 0; q < dim; q++)
Kirill Terekhov's avatar
Kirill Terekhov committed
744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767
					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
768
				for(integer j = 0; j < dim*2; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
769 770 771 772 773 774 775 776 777 778 779 780 781
					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
782 783
					if( mpirank == 0 ) SetStatus(*it,Element::Shared);
					else SetStatus(*it,Element::Ghost);
Kirill Terekhov's avatar
Kirill Terekhov committed
784 785 786 787 788 789 790 791 792
				}
				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
793 794
							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
795 796 797 798 799
					}
				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
800 801
							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
802
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
803
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
804 805 806 807 808 809
			}
			else
			{
				double time = Timer();
				Storage::real epsilon = GetEpsilon();
			
Kirill Terekhov's avatar
Kirill Terekhov committed
810
				GeomParam table;
Kirill Terekhov's avatar
Kirill Terekhov committed
811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830
				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
831
				buffer_type exch_data;
Kirill Terekhov's avatar
Kirill Terekhov committed
832 833
				std::vector< INMOST_DATA_REAL_TYPE > unpack_real;
				std::vector< INMOST_DATA_REAL_TYPE > pack_real;
Kirill Terekhov's avatar
Kirill Terekhov committed
834
				element_set sorted_nodes;
Kirill Terekhov's avatar
Kirill Terekhov committed
835 836 837 838 839 840 841 842 843
				
				sorted_nodes.reserve(NumberOfNodes());
				
				
				time = Timer();
				
				
				for(iteratorNode n = BeginNode(); n != EndNode(); n++)
				{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
844 845
					real_array c = n->Coords();
					for(real_array::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
846
						if( point_in_bbox(c.data(),bboxs.data()+procs[k]*dim*2,dim) )
Kirill Terekhov's avatar
Kirill Terekhov committed
847
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
848
							sorted_nodes.push_back(*n);
Kirill Terekhov's avatar
Kirill Terekhov committed
849 850 851 852 853 854 855 856 857 858 859 860 861 862
							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
863 864
					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
865 866 867 868 869 870 871 872 873
				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
874
				for(element_set::iterator it = sorted_nodes.begin(); it != sorted_nodes.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
875
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
876
					Storage::real_array arr = RealArrayDF(*it,CoordsTag());
Kirill Terekhov's avatar
Kirill Terekhov committed
877 878 879 880 881 882 883 884 885
					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
886
				MPI_Pack_size(static_cast<int>(pack_real.size()),INMOST_MPI_DATA_REAL_TYPE,comm,&sendsize);
Kirill Terekhov's avatar
Kirill Terekhov committed
887 888
				exch_data.resize(sendsize);
				int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
889
				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
890 891 892 893 894 895 896 897 898 899 900
				
				
				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
901
//~ #if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
902 903 904 905 906 907 908 909 910
					//~ 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
911
					unsigned long usend[2] = {static_cast<unsigned>(sendsize),static_cast<unsigned>(pack_real.size())};
Kirill Terekhov's avatar
Kirill Terekhov committed
912
					MPI_Pack_size(2,MPI_UNSIGNED,comm,&pack_size2);
Kirill Terekhov's avatar
Kirill Terekhov committed
913
					for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
914 915 916 917
					{
						send_buffs[k].first = procs[k];
						send_buffs[k].second.resize(pack_size2);
						position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
918
						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
919 920 921 922 923 924 925 926 927
						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
928
							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
929 930 931
						}
					}
					if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
932
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
933
						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
934
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
935 936 937 938 939 940
					//~ 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
941
						for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963
						{
							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
964
								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
965 966 967 968
								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
969
									for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
970 971 972 973 974 975 976 977 978 979 980 981 982
										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
983
										IntegerArrayDV(sorted_nodes[(it1-pack_real.begin())/dim],tag_processors).push_back(recv_buffs[*qt].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
984 985 986 987 988 989 990 991 992 993 994
										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
995
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
996
							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
997
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
998 999 1000 1001 1002 1003 1004
					}			
					
					time2 = Timer() - time2;
					REPORT_STR("Intersect all coordinates");
					REPORT_VAL("time",time2);
					
					time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
1005
					Element::Status estat;
Kirill Terekhov's avatar
Kirill Terekhov committed
1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022
					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
1023 1024 1025 1026
							if( v.size() == 1 ) 
								estat = Element::Owned;
							else 
								estat = Element::Shared;
Kirill Terekhov's avatar
Kirill Terekhov committed
1027 1028
						}
						else
Kirill Terekhov's avatar
Kirill Terekhov committed
1029 1030
							estat = Element::Ghost;
						SetStatus(*it,estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1031 1032 1033 1034 1035 1036
					}
					
					ComputeSharedProcs();
					RecomputeParallelStorage(NODE);
					AssignGlobalID(NODE);
					
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1037
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1038 1039 1040
					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
1041 1042
							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
1043 1044 1045 1046
					}
					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
1047 1048
							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
1049
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1050
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072
					
					
					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
1073
							determine_my_procs_low(this,*it, result, intersection);
Kirill Terekhov's avatar
Kirill Terekhov committed
1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
							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
1089
						std::vector<std::pair<int,int> > mapping;
Kirill Terekhov's avatar
Kirill Terekhov committed
1090 1091
						for(Mesh::iteratorElement it = BeginElement(current_mask >> 1); it != EndElement(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1092
							mapping.push_back(std::make_pair(it->GlobalID(),it->LocalID()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1093 1094
						}
						if( !mapping.empty() ) 
Kirill Terekhov's avatar
Kirill Terekhov committed
1095
							std::sort(mapping.begin(),mapping.end(),MappingComparator());
Kirill Terekhov's avatar
Kirill Terekhov committed
1096 1097 1098 1099 1100 1101
						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
1102
						Storage::integer_array procs = IntegerArrayDV(GetHandle(),tag_processors);
Kirill Terekhov's avatar
Kirill Terekhov committed
1103 1104 1105
						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
1106
						std::vector< element_set > elements(procs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1107 1108 1109 1110 1111 1112 1113
						
						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
1114
							int m = static_cast<int>(p-procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1115 1116 1117 1118 1119 1120 1121 1122 1123
							{
								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
1124
										Element::adj_type & sub = LowConn(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1125
										if( sub.size() == 0 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1126 1127 1128
										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
1129
										message_send[1]++;
Kirill Terekhov's avatar
Kirill Terekhov committed
1130
										elements[m].push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1131 1132
									}
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
1133 1134
								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
1135 1136 1137
								send_buffs[m].first = *p;
								send_buffs[m].second.resize(sendsize);
								int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1138
								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
1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164
								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
1165
										pos = static_cast<int>(p - procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1166 1167 1168
										break;
									}
								if( pos == -1 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1169
								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
1170
								message_recv[pos].resize(size-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1171
								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
1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182
							}
						}
						
						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
1183 1184
							int m = static_cast<int>(p-procs.begin());
							element_set remote_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1185 1186 1187 1188 1189 1190 1191
							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
1192
								dynarray<HandleType,64> sub_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1193 1194 1195
								for(int j = 0; j < conn_size; j++)
								{
									int global_id = message_recv[m][pos++];
Kirill Terekhov's avatar
Kirill Terekhov committed
1196 1197 1198 1199
									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
1200 1201 1202 1203 1204 1205
									if( find == -1 ) 
									{
										flag = 0;
										pos += conn_size-j-1;
										break;
									}
Kirill Terekhov's avatar
Kirill Terekhov committed
1206 1207
									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
1208 1209 1210 1211
									
								}
								if( flag )
								{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1212
									HandleType e = FindSharedAdjacency(sub_elements.data(),static_cast<enumerator>(sub_elements.size()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1213
									if( e == InvalidHandle() ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
1214 1215 1216 1217 1218
									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
1219 1220
							element_set result;
							element_set::iterator set_end;
Kirill Terekhov's avatar
Kirill Terekhov committed
1221 1222 1223 1224 1225 1226
							
							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
1227
							for(element_set::iterator qt = result.begin(); qt != result.end(); qt++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1228
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
1229 1230
								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
1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260
								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
1261
							SetStatus(*it, estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1262 1263 1264
						}
						RecomputeParallelStorage(current_mask);
						if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1265
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1266
							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
1267
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1268
						AssignGlobalID(current_mask);					
Kirill Terekhov's avatar
Kirill Terekhov committed
1269
						integer num = ElementNum(current_mask);
Kirill Terekhov's avatar
Kirill Terekhov committed
1270 1271 1272
#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
1273 1274
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1275 1276 1277
						}
						for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1278 1279
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1280
						}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1281
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292
	
						time = Timer() - time;
						REPORT_STR("Set parallel info");
						REPORT_VAL("time",time);
					}
				}
				RemoveGeometricData(table);	
			}
		}
		//RemoveGhost();
		
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1293
#else //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1294
		AssignGlobalID(CELL | FACE | EDGE | NODE);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1295
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
1296 1297 1298 1299 1300 1301 1302 1303 1304
		EXIT_FUNC();
	}
	
	
	void Mesh::RemoveGhost()
	{
		if( m_state == Mesh::Serial ) return;
		ENTER_FUNC()
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
1305
		element_set delete_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1306
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1307
		proc_elements del_shared, del_ghost;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1308
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1309 1310 1311
		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
1312
		for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1313
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1314
			Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1315 1316 1317
			if( estat == Element::Ghost ) 
			{
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1318
				del_ghost[it->IntegerDF(tag_owner)].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1319
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1320 1321 1322 1323
				delete_elements.push_back(*it);
			}
			else if( estat == Element::Shared )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1324 1325
				SetStatus(*it,Element::Owned);
				Storage::integer_array it_procs = it->IntegerArrayDV(tag_processors);
Kirill Terekhov's avatar
Kirill Terekhov committed
1326 1327 1328 1329 1330
				
#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
1331
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1332 1333 1334 1335 1336 1337 1338
				
				it_procs.resize(1);
				it_procs[0] = mpirank;
			}
		}
		
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1339
		for(proc_elements::iterator it = del_ghost.begin(); it != del_ghost.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1340
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354
			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
1355 1356 1357 1358
			result.resize(end-result.begin());
			ref.swap(result);
		}
		del_ghost.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1359
		for(proc_elements::iterator it = del_shared.begin(); it != del_shared.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1360
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374
			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
1375 1376 1377 1378
			result.resize(end-result.begin());
			ref.swap(result);
		}
		del_shared.clear();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1379
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1380

Kirill Terekhov's avatar
Kirill Terekhov committed
1381
		for(element_set::iterator it = delete_elements.begin(); it != delete_elements.end(); it++) Destroy(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1382
		
Kirill Terekhov's avatar
Kirill Terekhov committed
1383
		delete_elements.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1384 1385 1386 1387
		for(ElementType mask = FACE; mask >= NODE; mask = mask >> 1)
		{
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1388
				Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1389 1390 1391
				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
1392
			ReduceData(tag_delete,mask,0,DeleteUnpack);
Kirill Terekhov's avatar
Kirill Terekhov committed
1393
			ExchangeData(tag_delete,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
1394 1395
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1396
				Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1397 1398 1399 1400
				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
1401
					del_ghost[it->IntegerDF(tag_owner)].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1402
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1403
					delete_elements.push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1404 1405 1406 1407 1408 1409 1410 1411 1412
				}
				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
1413
						del_shared[*vit].push_back(*it);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1414
#endif //USE_PARALLEL_STORAGE				
Kirill Terekhov's avatar
Kirill Terekhov committed
1415 1416 1417 1418 1419 1420
					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
1421
						SetStatus(*it,Element::Owned);
Kirill Terekhov's avatar
Kirill Terekhov committed
1422 1423 1424 1425
				}
			}
			
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1426
			for(proc_elements::iterator it = del_ghost.begin(); it != del_ghost.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1427
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441
				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
1442 1443 1444 1445
				result.resize(end-result.begin());
				ref.swap(result);
			}
			del_ghost.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
1446
			for(proc_elements::iterator it = del_shared.begin(); it != del_shared.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1447
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461
				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
1462 1463 1464 1465
				result.resize(end-result.begin());
				ref.swap(result);
			}
			del_shared.clear();
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1466
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1467
			for(element_set::iterator it = delete_elements.begin(); it != delete_elements.end(); it++) Destroy(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1468 1469 1470 1471
			delete_elements.clear();
		}		
		DeleteTag(tag_delete);
		ComputeSharedProcs();
Kirill Terekhov's avatar
Kirill Terekhov committed
1472 1473
		Integer(GetHandle(),tag_layers) = 0;
		Integer(GetHandle(),tag_bridge) = NONE;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1474
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
1475 1476 1477
		EXIT_FUNC();
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1478
	void Mesh::RemoveGhostElements(const HandleType * ghost, enumerator num)
Kirill Terekhov's avatar
Kirill Terekhov committed
1479 1480 1481 1482 1483
	{
		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
1484
		element_set delete_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1485
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1486
		proc_elements del_shared, del_ghost;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1487
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1488 1489 1490 1491 1492
		
		for(ElementType mask = CELL; mask >= NODE; mask = mask >> 1)
		{
			if( mask & CELL )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1493 1494
				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
1495 1496 1497 1498 1499
			}
			else 
			{
				for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1500
					Element::Status estat = GetStatus(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1501 1502 1503 1504
					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
1505 1506
			ReduceData(tag_delete,mask,0,DeleteUnpack);
			ExchangeData(tag_delete,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
1507 1508
			for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1509
				Element::Status estat = GetStatus(*it);