parallel.cpp 168 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
	//////////////////////////////
	/// 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 )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
79
      const Storage::integer * recv = static_cast<const Storage::integer *>(static_cast<const void *>(data));
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
80 81 82 83 84 85 86 87 88 89
			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 )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
      for(INMOST_DATA_ENUM_TYPE k = 0; k < size; k+=2 )
      {
			  bool flag = true;
			  const Storage::integer * recv = static_cast<const Storage::integer *>(static_cast<const void *>(data))+k;
			  Storage::integer_array arr = e->IntegerArray(tag);
      
        for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it+=2)
        {
          if( *it == recv[0] )
				  {
					  flag = false;
					  break;
				  }
        }
			  if( flag ) 
			  {
				  arr.push_back(recv[0]);
				  arr.push_back(recv[1]);
			  }
      }
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
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 148 149 150 151 152 153
		}
	}
	


	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
154

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

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

	
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
168
	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
169
	{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
170
		(void) size;
Kirill Terekhov's avatar
Kirill Terekhov committed
171 172 173
		element->Bulk(tag) &= *data;
	}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
174 175 176 177 178
	//////////////////////////////
	/// REDUCTION FUNCTIONS END///
	//////////////////////////////


Kirill Terekhov's avatar
Kirill Terekhov committed
179

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

Kirill Terekhov's avatar
Kirill Terekhov committed
192
			for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
193
				if( it->GetMarker(marker) && (it->GetStatus() & Expr) )
Kirill Terekhov's avatar
Kirill Terekhov committed
194 195 196 197 198 199
					it->Bulk(t) = 1;

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

	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
257
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
258 259
		return etypeout;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
260 261 262 263 264 265 266 267
	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
268
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
269 270
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
271
	Storage::integer Mesh::EnumerateSet(const ElementSet & set, const Tag & num_tag, Storage::integer start, bool define_sparse)
Kirill Terekhov's avatar
Kirill Terekhov committed
272 273 274 275 276
	{
		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
277 278
		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
279 280 281
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
282
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
283
		shift += start;
Kirill Terekhov's avatar
Kirill Terekhov committed
284 285
		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
286
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
287 288
		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
289
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
290
		ExchangeData(num_tag,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
291 292 293 294
		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
295
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
296 297
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
298
	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
299 300 301 302 303
	{
		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
304 305
		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
306 307
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
308
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
309
		shift += start;
Kirill Terekhov's avatar
Kirill Terekhov committed
310
		for(const HandleType * it = set; it != set+n; ++it)
311
			if( GetStatus(*it) == Element::Owned && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
312
		for(const HandleType * it = set; it != set+n; ++it) 
313
			if( GetStatus(*it) == Element::Shared && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
314
		ExchangeData(num_tag,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
315 316 317 318
		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
319
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
320 321
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
322
	Storage::integer Mesh::Enumerate(ElementType mask, Tag num_tag, Storage::integer start, bool define_sparse)
Kirill Terekhov's avatar
Kirill Terekhov committed
323 324 325 326 327
	{
		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
328
			if( it->GetStatus() != Element::Ghost && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
329 330 331
				number++;
		MPI_Scan(&number,&shift,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
		shift -= number;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
332
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
333 334
		shift += start;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
335
			if( it->GetStatus() == Element::Owned && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
336 337
				it->Integer(num_tag) = shift++;
		for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++)
338
			if( it->GetStatus() == Element::Shared && (define_sparse || it->HaveData(num_tag)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
339
				it->Integer(num_tag) = shift++;
Kirill Terekhov's avatar
Kirill Terekhov committed
340
		ExchangeData(num_tag,mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
341 342 343 344
		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
345
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
346 347 348 349 350 351 352 353
		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
354 355 356
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
357 358
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
359 360 361


	
Kirill Terekhov's avatar
Kirill Terekhov committed
362 363 364 365 366 367
	
	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
368 369 370
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
371 372
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
373 374 375 376

	void Mesh::Integrate(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
377 378
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
379 380 381 382 383 384 385 386 387 388
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::Integrate(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
389 390
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
391 392 393 394 395 396
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
397 398 399 400 401 402 403
	
	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
404 405 406
#else//USE_MPI
		(void) input;
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
407 408 409
		return output;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
410
	Storage::real Mesh::Integrate(const Tag & t, enumerator entry, ElementType mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
411 412 413
	{
		Storage::real output = 0, input = 0;
		for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
414 415 416 417 418
			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
419 420 421
		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
422
#endif//USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
423 424 425
		return output;
	}
	
426 427 428 429 430
	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
431 432 433
#else //USE_MPI
		(void) input;
#endif //USE_MPI
434 435
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
436
	
437 438 439 440 441
	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
442 443 444
#else //USE_MPI
		(void) input;
#endif //USE_MPI
445 446
		return output;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
447

Kirill Terekhov's avatar
Kirill Terekhov committed
448 449 450
	void Mesh::AggregateMax(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
451 452
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
453 454 455 456 457 458 459 460 461 462
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::AggregateMax(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
463 464
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
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 496
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

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

	void Mesh::AggregateMin(Storage::real * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
497 498
		static dynarray<Storage::real,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
499 500 501 502 503 504 505 506 507 508
		memcpy(temp.data(),input,sizeof(Storage::real)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MIN,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

	void Mesh::AggregateMin(Storage::integer * input, Storage::integer size)
	{
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
509 510
		static dynarray<Storage::integer,64> temp;
		temp.resize(size);
Kirill Terekhov's avatar
Kirill Terekhov committed
511 512 513 514 515 516 517
		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MIN,comm);
#else//USE_MPI
		(void) input;
#endif//USE_MPI
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
518

Kirill Terekhov's avatar
Kirill Terekhov committed
519
	INMOST_MPI_Comm Mesh::GetCommunicator() const
Kirill Terekhov's avatar
Kirill Terekhov committed
520 521 522
	{
#if defined(USE_MPI)
		return comm;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
523
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
524
		return 0;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
525
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
526 527
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
528
	int Mesh::GetProcessorRank() const
Kirill Terekhov's avatar
Kirill Terekhov committed
529 530 531 532 533
	{
#if defined(USE_MPI)
		int rank;
		MPI_Comm_rank(comm,&rank);
		return rank;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
534
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
535
		return 0;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
536
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
537 538
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
539
	int Mesh::GetProcessorsNumber() const
Kirill Terekhov's avatar
Kirill Terekhov committed
540 541 542 543 544
	{
#if defined(USE_MPI)
		int size;
		MPI_Comm_size(comm,&size);
		return size;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
545
#else //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
546
		return 1;
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
547
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
548 549 550 551 552 553 554 555
	}	
	
	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
556 557 558 559
#else //USE_MPI
		(void) argc;
		(void) argv;
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
560 561 562
#if defined(USE_PARALLEL_WRITE_TIME)
		atexit(Mesh::AtExit);
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
563 564 565 566 567 568 569 570
	}
	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
571
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
572 573
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
574
    INMOST_MPI_Group Mesh::GetGroup() const
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
575 576 577 578 579 580 581
    {
        INMOST_MPI_Group ret = INMOST_MPI_GROUP_EMPTY;
#if defined(USE_MPI)
        MPI_Comm_group(GetCommunicator(), &ret);
#endif
        return ret;
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
582 583 584 585
	
	void Mesh::SetCommunicator(INMOST_MPI_Comm _comm)
	{
		ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
586 587 588 589 590 591
		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
592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611

#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
612
#if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630
		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
631 632 633 634
#endif //USE_MPI_P2P
#else //USE_MPI
		(void) _comm;
#endif //USE_MPI
Kirill Terekhov's avatar
Kirill Terekhov committed
635 636 637 638 639 640
		EXIT_FUNC();
	}
	
#if defined(USE_PARALLEL_WRITE_TIME)
	void Mesh::Enter() { tab++; }
	void Mesh::Exit() {tab--; }
Kirill Terekhov's avatar
Kirill Terekhov committed
641
	std::ostream & Mesh::WriteTab(std::ostream & f)
Kirill Terekhov's avatar
Kirill Terekhov committed
642 643 644 645 646 647 648 649 650
	{
		for(int i = 0; i < tab; i++)
			f << "   ";
		return f;
	}
	std::fstream & Mesh::GetStream()
	{
		return out_time;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
651 652
	void Mesh::FinalizeFile()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
653 654
    //std::stringstream str;
    if( tab > 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
655
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
656
      out_time << "<TEXT><![CDATA[Died!]]></TEXT>\n";// << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
657
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
658
    while(tab > 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
659
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
660 661
      out_time << "<TIME>-1</TIME>\n</FUNCTION>\n";// << std::endl; 
      Exit(); 
Kirill Terekhov's avatar
Kirill Terekhov committed
662 663
		}
		out_time << "</Debug>" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
664 665 666
    //out_time << str;
    //out_time.flush();
    out_time.close();
Kirill Terekhov's avatar
Kirill Terekhov committed
667
	}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
668
#endif //USE_PARALLEL_WRITE_TIME
Kirill Terekhov's avatar
Kirill Terekhov committed
669 670 671
	
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
672
	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
673
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
674 675 676
		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
677 678 679 680 681
		result.clear();
		result.insert(result.end(),p.begin(),p.end());
		i++;
		while(i != subelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
682 683
			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
684 685 686 687 688 689 690
			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
691
	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
692
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
693
		Element::adj_type const & overelements = m->HighConn(h);
Kirill Terekhov's avatar
Kirill Terekhov committed
694
		if( overelements.empty() ) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
695
		Element::adj_type::const_iterator i = overelements.begin();
Kirill Terekhov's avatar
Kirill Terekhov committed
696 697 698
		result.clear();
		while(i != overelements.end())
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
699
			Storage::integer_array q = m->IntegerArrayDV(*i,procs);
Kirill Terekhov's avatar
Kirill Terekhov committed
700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806
			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++)
Kirill Terekhov's avatar
Kirill Terekhov committed
807
			if( ::fabs(a[i]-b[i]) > eps ) return a[i] <= b[i];
Kirill Terekhov's avatar
Kirill Terekhov committed
808 809
		return true;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
810 811 812 813 814

	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
815 816 817 818 819 820
	
	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
821
		integer dim = GetDimensions();
Kirill Terekhov's avatar
Kirill Terekhov committed
822 823 824 825 826
		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
827
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
828 829
		//determine which bboxes i intersect
		dynarray<int,64> procs;
Kirill Terekhov's avatar
Kirill Terekhov committed
830 831 832 833 834
		Storage::real bbox[6]; //local bounding box
		std::vector<Storage::real> bboxs(mpisize*6);
		//Compute local bounding box containing nodes.
		//Will be more convinient to compute (or store)
		//and communicate local octree over all the nodes.
Kirill Terekhov's avatar
Kirill Terekhov committed
835
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
836 837 838 839
		{
			bbox[k] = 1e20;
			bbox[k+dim] = -1e20;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
840 841 842
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
843
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
844
			real bbox0[6];
Kirill Terekhov's avatar
Kirill Terekhov committed
845
			for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
846
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871
				bbox0[k] = 1e20;
				bbox0[k+dim] = -1e20;
			}
#if defined(USE_OMP)
#pragma omp for
#endif
			for(integer nit = 0; nit < NodeLastLocalID(); ++nit) if( isValidNode(nit) )
			{
				Node it = NodeByLocalID(nit);
				Storage::real_array arr = it->Coords();
				for(integer k = 0; k < dim; k++)
				{
					if( arr[k] < bbox0[k] ) bbox0[k] = arr[k];
					if( arr[k] > bbox0[k+dim] ) bbox0[k+dim] = arr[k];
				}
			}
#if defined(USE_OMP)
#pragma omp critical
#endif
			{
				for(integer k = 0; k < dim; k++)
				{
					if( bbox0[k] < bbox[k] ) bbox[k] = bbox0[k];
					if( bbox0[k] > bbox[k+dim] ) bbox[k+dim] = bbox0[k];
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
872 873
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
874
		// write down bounding boxes
Kirill Terekhov's avatar
Kirill Terekhov committed
875
		for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
876 877 878 879
		{
			REPORT_VAL("min",bbox[k]);
			REPORT_VAL("max",bbox[dim+k]);
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
880
		// communicate bounding boxes
Kirill Terekhov's avatar
Kirill Terekhov committed
881
		REPORT_MPI(MPI_Allgather(&bbox[0],dim*2,INMOST_MPI_DATA_REAL_TYPE,&bboxs[0],dim*2,INMOST_MPI_DATA_REAL_TYPE,comm));
Kirill Terekhov's avatar
Kirill Terekhov committed
882
		// find all processors that i communicate with
Kirill Terekhov's avatar
Kirill Terekhov committed
883 884 885 886
		for(int k = 0; k < mpisize; k++)
			if( k != mpirank )
			{
				bool flag = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
887
				for(integer q = 0; q < dim; q++)
Kirill Terekhov's avatar
Kirill Terekhov committed
888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911
					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
912
				for(integer j = 0; j < dim*2; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
913
					same_box &= ::fabs(bbox[j] - bboxs[k*dim*2+j]) < epsilon;
Kirill Terekhov's avatar
Kirill Terekhov committed
914 915 916 917 918 919
				same_boxes &= same_box;
			}
			
			if( same_boxes )
			{
				REPORT_STR("All bounding boxes are the same - assuming that mesh is replicated over all nodes");
Kirill Terekhov's avatar
Kirill Terekhov committed
920 921
				//for(Mesh::iteratorElement it = BeginElement(CELL | EDGE | FACE | NODE); it != EndElement(); it++)
				for(ElementType etype = NODE; etype <= CELL; etype = NextElementType(etype) )
Kirill Terekhov's avatar
Kirill Terekhov committed
922
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938
#if defined(USE_OMP)
#pragma omp parallel for
#endif
					for(integer eit = 0; eit < LastLocalID(etype); ++eit) if( isValidElement(etype,eit) )
					{
						Element it = ElementByLocalID(etype,eit);
						integer_array arr = it->IntegerArrayDV(tag_processors);
						arr.resize(mpisize);
						for(int k = 0; k < mpisize; k++)
							arr[k] = k;
						it->IntegerDF(tag_owner) = 0;
						if( mpirank == 0 )
							SetStatus(it->GetHandle(),Element::Shared);
						else
							SetStatus(it->GetHandle(),Element::Ghost);
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
939 940 941 942 943 944 945 946 947
				}
				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
948 949
							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
950
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
951
				for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
952 953 954
					for(int i = 0; i < 4; i++)
					{
						if( !it->second[i].empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
955 956
							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
957
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
958
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
959 960 961 962 963 964
			}
			else
			{
				double time = Timer();
				Storage::real epsilon = GetEpsilon();
			
Kirill Terekhov's avatar
Kirill Terekhov committed
965
				GeomParam table;
Kirill Terekhov's avatar
Kirill Terekhov committed
966 967 968 969 970 971 972 973 974 975 976
				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);
			
				
Kirill Terekhov's avatar
Kirill Terekhov committed
977 978 979 980 981 982
				
				//for(iteratorNode it = BeginNode(); it != EndNode(); it++)
#if defined(USE_OMP)
#pragma omp parallel for
#endif
				for(integer nit = 0; nit < NodeLastLocalID(); ++nit) if( isValidNode(nit) )
Kirill Terekhov's avatar
Kirill Terekhov committed
983
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
984
					Node it = NodeByLocalID(nit);
Kirill Terekhov's avatar
Kirill Terekhov committed
985 986 987 988 989 990
					Storage::integer_array arr = it->IntegerArrayDV(tag_processors);
					arr.resize(1);
					arr[0] = mpirank;
				}
			

Kirill Terekhov's avatar
Kirill Terekhov committed
991
				buffer_type exch_data;
Kirill Terekhov's avatar
Kirill Terekhov committed
992 993
				std::vector< INMOST_DATA_REAL_TYPE > unpack_real;
				std::vector< INMOST_DATA_REAL_TYPE > pack_real;
Kirill Terekhov's avatar
Kirill Terekhov committed
994
				element_set sorted_nodes;
Kirill Terekhov's avatar
Kirill Terekhov committed
995 996 997 998 999 1000 1001 1002 1003
				
				sorted_nodes.reserve(NumberOfNodes());
				
				
				time = Timer();
				
				
				for(iteratorNode n = BeginNode(); n != EndNode(); n++)
				{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1004 1005
					real_array c = n->Coords();
					for(real_array::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1006
						if( point_in_bbox(c.data(),bboxs.data()+procs[k]*dim*2,dim) )
Kirill Terekhov's avatar
Kirill Terekhov committed
1007
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1008
							sorted_nodes.push_back(*n);
Kirill Terekhov's avatar
Kirill Terekhov committed
1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022
							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
1023 1024
					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
1025 1026 1027 1028 1029 1030 1031 1032 1033
				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
1034
				for(element_set::iterator it = sorted_nodes.begin(); it != sorted_nodes.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1035
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1036
					Storage::real_array arr = RealArrayDF(*it,CoordsTag());
Kirill Terekhov's avatar
Kirill Terekhov committed
1037 1038 1039 1040 1041 1042 1043 1044 1045
					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
1046
				MPI_Pack_size(static_cast<int>(pack_real.size()),INMOST_MPI_DATA_REAL_TYPE,comm,&sendsize);
Kirill Terekhov's avatar
Kirill Terekhov committed
1047 1048
				exch_data.resize(sendsize);
				int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1049
				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
1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060
				
				
				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
1061
//~ #if defined(USE_MPI_P2P)
Kirill Terekhov's avatar
Kirill Terekhov committed
1062 1063
					//~ unsigned * sendsizeall = shared_space;
					//~ unsigned usend[2] = {sendsize,pack_real.size()};
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1064
					//~ REPORT_MPI(MPI_Win_fence(0,window)); //start exchange session
Kirill Terekhov's avatar
Kirill Terekhov committed
1065 1066
					//~ for(unsigned k = 0; k < procs.size(); k++)
						//~ REPORT_MPI(MPI_Put(usend,2,MPI_UNSIGNED,procs[k],mpirank*2,2,MPI_UNSIGNED,window));
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1067
					//~ REPORT_MPI(MPI_Win_fence(MPI_MODE_NOSTORE | MPI_MODE_NOSUCCEED,window)); //end exchange session
Kirill Terekhov's avatar
Kirill Terekhov committed
1068 1069 1070
//~ #else
					std::vector<unsigned> sendsizeall(mpisize*2);
					int pack_size2 = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1071
					unsigned usend[2] = {static_cast<unsigned>(sendsize),static_cast<unsigned>(pack_real.size())};
Kirill Terekhov's avatar
Kirill Terekhov committed
1072
					MPI_Pack_size(2,MPI_UNSIGNED,comm,&pack_size2);
Kirill Terekhov's avatar
Kirill Terekhov committed
1073
					for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1074 1075 1076 1077
					{
						send_buffs[k].first = procs[k];
						send_buffs[k].second.resize(pack_size2);
						position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1078
						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
1079 1080 1081 1082 1083 1084 1085 1086 1087
						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
1088
							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
1089 1090 1091
						}
					}
					if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1092
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1093
						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
1094
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1095 1096 1097 1098 1099 1100
					//~ 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
1101
						for(dynarray<integer,64>::size_type k = 0; k < procs.size(); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123
						{
							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
1124
								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
1125 1126 1127 1128
								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
1129
									for(integer k = 0; k < dim; k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1130
										if( ::fabs((*(it1+k))-(*(it2+k))) > epsilon )
Kirill Terekhov's avatar
Kirill Terekhov committed
1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142
										{
											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
1143
										IntegerArrayDV(sorted_nodes[(it1-pack_real.begin())/dim],tag_processors).push_back(recv_buffs[*qt].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154
										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
1155
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1156
							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
1157
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1158 1159 1160 1161 1162 1163 1164
					}			
					
					time2 = Timer() - time2;
					REPORT_STR("Intersect all coordinates");
					REPORT_VAL("time",time2);
					
					time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
1165
					Element::Status estat;
Kirill Terekhov's avatar
Kirill Terekhov committed
1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182
					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
1183 1184 1185 1186
							if( v.size() == 1 ) 
								estat = Element::Owned;
							else 
								estat = Element::Shared;
Kirill Terekhov's avatar
Kirill Terekhov committed
1187 1188
						}
						else
Kirill Terekhov's avatar
Kirill Terekhov committed
1189 1190
							estat = Element::Ghost;
						SetStatus(*it,estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1191 1192 1193 1194 1195 1196
					}
					
					ComputeSharedProcs();
					RecomputeParallelStorage(NODE);
					AssignGlobalID(NODE);
					
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1197
#if defined(USE_PARALLEL_STORAGE)
Kirill Terekhov's avatar
Kirill Terekhov committed
1198 1199 1200
					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
1201 1202
							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
1203 1204 1205 1206
					}
					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
1207 1208
							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
1209
					}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1210
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar
Kirill Terekhov committed
1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225
					
					
					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));
						
						
Kirill Terekhov's avatar
Kirill Terekhov committed
1226 1227
						//int owned_elems = 0;
						//int shared_elems = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1228 1229 1230 1231 1232
						int owner;
						Element::Status estat;
						
						time = Timer();
						//Determine what processors potentially share the element
Kirill Terekhov's avatar
Kirill Terekhov committed
1233 1234 1235 1236
#if defined(USE_OMP)
#pragma omp parallel for
#endif
						for(integer eit = 0; eit < LastLocalID(current_mask); ++eit)
Kirill Terekhov's avatar
Kirill Terekhov committed
1237
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1238
							if( isValidElement(current_mask,eit) )
Kirill Terekhov's avatar
Kirill Terekhov committed
1239
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255
								Element it = ElementByLocalID(current_mask,eit);
								determine_my_procs_low(this,it->GetHandle(), result, intersection);
								Storage::integer_array p = it->IntegerArrayDV(tag_processors);
								if( result.empty() )
								{
									p.clear();
									p.push_back(mpirank);
									//++owned_elems;
								}
								else
								{
									p.replace(p.begin(),p.end(),result.begin(),result.end());
									//if( result.size() == 1 && result[0] == mpirank )
									//  ++owned_elems;
									//else ++shared_elems;
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
1256 1257 1258
							}
						}
						time = Timer() - time;
Kirill Terekhov's avatar
Kirill Terekhov committed
1259 1260
						//REPORT_VAL("predicted owned elements",owned_elems);
						//REPORT_VAL("predicted shared elements",shared_elems);
Kirill Terekhov's avatar
Kirill Terekhov committed
1261 1262 1263 1264 1265 1266
						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
1267
						std::vector<std::pair<int,int> > mapping;
Kirill Terekhov's avatar
Kirill Terekhov committed
1268
						REPORT_VAL("mapping type",ElementTypeName(current_mask >> 1));
Kirill Terekhov's avatar
Kirill Terekhov committed
1269 1270
						for(Mesh::iteratorElement it = BeginElement(current_mask >> 1); it != EndElement(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1271
							mapping.push_back(std::make_pair(it->GlobalID(),it->LocalID()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1272 1273
						}
						if( !mapping.empty() ) 
Kirill Terekhov's avatar
Kirill Terekhov committed
1274
							std::sort(mapping.begin(),mapping.end(),MappingComparator());
Kirill Terekhov's avatar
Kirill Terekhov committed
1275
						time = Timer() - time;
Kirill Terekhov's avatar
Kirill Terekhov committed
1276
            REPORT_VAL("mapping size",mapping.size())
Kirill Terekhov's avatar
Kirill Terekhov committed
1277 1278 1279 1280 1281
						REPORT_STR("Compute global to local indexes mapping");
						REPORT_VAL("time",time);
						//Initialize arrays
						
						time = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
1282
						Storage::integer_array procs = IntegerArrayDV(GetHandle(),tag_processors);
Kirill Terekhov's avatar
Kirill Terekhov committed
1283 1284 1285
						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
1286
						std::vector< element_set > elements(procs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1287 1288 1289 1290 1291 1292 1293
						
						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
1294
							int m = static_cast<int>(p-procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1295 1296 1297 1298 1299 1300 1301 1302 1303
							{
								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
1304
										Element::adj_type & sub = LowConn(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1305
										if( sub.size() == 0 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1306 1307 1308
										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
1309
										message_send[1]++;
Kirill Terekhov's avatar
Kirill Terekhov committed
1310
										elements[m].push_back(*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
1311 1312
									}
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
1313 1314
                REPORT_VAL("for processor",*p);
                REPORT_VAL("gathered elements",elements[m].size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1315 1316
								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
1317 1318 1319
								send_buffs[m].first = *p;
								send_buffs[m].second.resize(sendsize);
								int position = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1320
								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
1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346
								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
1347
										pos = static_cast<int>(p - procs.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1348 1349 1350
										break;
									}
								if( pos == -1 ) throw Impossible;
Kirill Terekhov's avatar
Kirill Terekhov committed
1351
								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
1352
                REPORT_VAL("unpacked message size",size-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1353
								message_recv[pos].resize(size-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1354
								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
1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365
							}
						}
						
						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
1366 1367
							int m = static_cast<int>(p-procs.begin());
							element_set remote_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1368 1369 1370 1371 1372 1373 1374
							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
1375
								dynarray<HandleType,64> sub_elements;
Kirill Terekhov's avatar
Kirill Terekhov committed
1376 1377 1378
								for(int j = 0; j < conn_size; j++)
								{
									int global_id = message_recv[m][pos++];
Kirill Terekhov's avatar
Kirill Terekhov committed
1379 1380 1381 1382
									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
1383 1384 1385 1386 1387 1388
									if( find == -1 ) 
									{
										flag = 0;
										pos += conn_size-j-1;
										break;
									}
Kirill Terekhov's avatar
Kirill Terekhov committed
1389
									int find_local_id = mapping[find].second;
Kirill Terekhov's avatar
Kirill Terekhov committed
1390
									sub_elements.push_back(ComposeHandle(PrevElementType(current_mask), find_local_id));
Kirill Terekhov's avatar
Kirill Terekhov committed
1391 1392 1393 1394
									
								}
								if( flag )
								{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1395
									HandleType e = FindSharedAdjacency(sub_elements.data(),static_cast<enumerator>(sub_elements.size()));
Kirill Terekhov's avatar
Kirill Terekhov committed
1396
									if( e == InvalidHandle() ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
1397 1398 1399
									remote_elements.push_back(e);
								}
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
1400 1401 1402 1403 1404 1405 1406 1407
              REPORT_VAL("number of unpacked remote elements",remote_elements.size());
              if( !remote_elements.empty() )
              {
                REPORT_VAL("first",remote_elements.front());
                REPORT_VAL("first type",ElementTypeName(GetHandleElementType(remote_elements.front())));
                REPORT_VAL("last",remote_elements.back());
                REPORT_VAL("last type",ElementTypeName(GetHandleElementType(remote_elements.back())));
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
1408
							std::sort(remote_elements.begin(),remote_elements.end());
Kirill Terekhov's avatar
Kirill Terekhov committed
1409 1410 1411 1412 1413 1414 1415 1416
              REPORT_VAL("original elements size",elements[m].size());
              if( !elements[m].empty() )
              {
                REPORT_VAL("first",elements[m].front());
                REPORT_VAL("first type",ElementTypeName(GetHandleElementType(elements[m].front())));
                REPORT_VAL("last",elements[m].back());
                REPORT_VAL("last type",ElementTypeName(GetHandleElementType(elements[m].back())));
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
1417
							std::sort(elements[m].begin(),elements[m].end());
Kirill Terekhov's avatar
Kirill Terekhov committed
1418 1419
							element_set result;
							element_set::iterator set_end;
Kirill Terekhov's avatar
Kirill Terekhov committed
1420 1421 1422
							
							result.resize(elements[m].size());
							set_end = std::set_difference(elements[m].begin(),elements[m].end(),remote_elements.begin(),remote_elements.end(), result.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1423 1424 1425
							result.resize(set_end-result.begin());

              REPORT_VAL("set difference size",result.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
1426 1427
							
							//elements in result are wrongly marked as ghost
Kirill Terekhov's avatar
Kirill Terekhov committed
1428
							for(element_set::iterator qt = result.begin(); qt != result.end(); qt++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1429
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
1430 1431
								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
1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461
								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
1462
							SetStatus(*it, estat);
Kirill Terekhov's avatar
Kirill Terekhov committed
1463 1464 1465
						}
						RecomputeParallelStorage(current_mask);
						if( !send_reqs.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1466
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1467
							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
1468
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1469
						AssignGlobalID(current_mask);					
Kirill Terekhov's avatar
Kirill Terekhov committed
1470
						integer num = ElementNum(current_mask);
Kirill Terekhov's avatar
Kirill Terekhov committed
1471 1472 1473
#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
1474 1475
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1476 1477 1478
						}
						for(parallel_storage::iterator it = ghost_elements.begin(); it != ghost_elements.end(); it++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1479 1480
							if( !it->second[num].empty() )
								std::sort(it->second[num].begin(),it->second[num].end(),GlobalIDComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1481
						}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1482
#endif //USE_PARALLEL_STORAGE
Kirill Terekhov's avatar