ParallelDebug.cc 30.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
#include "parallel/ParallelDebug.h"
#include "parallel/MeshDistributor.h"
#include "parallel/MpiHelper.h"
16
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
17
18
19
20
#include "DOFVector.h"
#include "FixVec.h"
#include "StdMpi.h"
#include "Debug.h"
21
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
22
23
24

namespace AMDiS {

25
26
27
  using namespace std;


28
  void ParallelDebug::testInteriorBoundary(MeshDistributor &pdb)
Thomas Witkowski's avatar
Thomas Witkowski committed
29
  {
30
    FUNCNAME("ParallelDebug::testInteriorBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
31

32
    vector<int*> sendBuffers, recvBuffers;
Thomas Witkowski's avatar
Thomas Witkowski committed
33

34
    MPI::Request request[pdb.rankIntBoundary.boundary.size() + 
35
36
			 pdb.otherIntBoundary.boundary.size() +
                         pdb.periodicBoundary.boundary.size() * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
37
38
    int requestCounter = 0;

39
40
41

    // === Send rank's boundary information. ===

42
43
    for (RankToBoundMap::iterator rankIt = pdb.rankIntBoundary.boundary.begin();
	 rankIt != pdb.rankIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
44
45
46
47
48
49
50
51
52
53
54
55

      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
	buffer[i] = (rankIt->second)[i].rankObj.elIndex;
      
      sendBuffers.push_back(buffer);
      
      request[requestCounter++] =
	pdb.mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
    }

56
57
58

    // === Receive information from other ranks about the interior boundaries. ====

Thomas Witkowski's avatar
Thomas Witkowski committed
59
60
61
62
63
64
65
66
67
68
69
    for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
	 rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
      recvBuffers.push_back(buffer);

      request[requestCounter++] = 
	pdb.mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
    }


70
71
    // === To the last, do the same of periodic boundaries. ===

72
    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();       
73
	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
74
75
76
      if (rankIt->first == pdb.mpiRank)
	continue;

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
      int nValues = rankIt->second.size();
      int* sBuffer = new int[nValues];
      for (int i = 0; i < nValues; i++)
	sBuffer[i] = (rankIt->second)[i].rankObj.elIndex;

      sendBuffers.push_back(sBuffer);

      request[requestCounter++] =
	pdb.mpiComm.Isend(sBuffer, nValues, MPI_INT, rankIt->first, 0);

      int *rBuffer = new int[nValues];
      recvBuffers.push_back(rBuffer);

      request[requestCounter++] = 
	pdb.mpiComm.Irecv(rBuffer, nValues, MPI_INT, rankIt->first, 0);      
    }

    // === Finish communication and delete all send buffers. ===

    MPI::Request::Waitall(requestCounter, request);
Thomas Witkowski's avatar
Thomas Witkowski committed
97
98
99
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];

100
101
102
103
104

    // === Finally, check the results, i.e., the indices of element at the     === 
    // === boundaries, if they fit together. First check the interior bounds,  ===
    // === and after this the periodic ones.                                   ===

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
107
108
    int bufCounter = 0;
    for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
	 rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {

109
110
      TEST_EXIT(rankIt->second.size() == 
		pdb.otherIntBoundary.boundary[rankIt->first].size())
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
	("Boundaries does not fit together!\n");      

113
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116
117
118
119
120
121
	int elIndex1 = recvBuffers[bufCounter][i];
	int elIndex2 = pdb.otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex;

	TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
      }

      delete [] recvBuffers[bufCounter++];
    }
122
123
124
125


    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
126
127
128
      if (rankIt->first == pdb.mpiRank)
	continue;

129
130
131
132
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
	int elIndex1 = recvBuffers[bufCounter][i];
	int elIndex2 = pdb.periodicBoundary.boundary[rankIt->first][i].neighObj.elIndex;

133
134
135
136
	TEST_EXIT(elIndex1 == elIndex2)
	  ("Wrong element index at periodic boundary el %d with rank %d: %d %d\n", 
	   pdb.periodicBoundary.boundary[rankIt->first][i].rankObj.elIndex,
	   rankIt->first, elIndex1, elIndex2);
137
138
139
140
      }

      delete [] recvBuffers[bufCounter++];
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
143
  }


144
145
146
147
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");

148
149
150
151
    for (unsigned int i = 0; i < pdb.feSpaces.size(); i++)
      testPeriodicBoundary(pdb, pdb.feSpaces[i]);
  }

152

153
154
155
156
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb,
					   const FiniteElemSpace *feSpace)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");
157

158
159
    // === 1. check: All periodic DOFs should have at least a correct number ===
    // === of periodic associations.                                         ===
160
161
162
163
164

    PeriodicMap &perMap = pdb.getPeriodicMap();
    for (map<int, std::set<BoundaryType> >::iterator it = 
	   perMap.periodicDofAssociations[feSpace].begin();
	 it != perMap.periodicDofAssociations[feSpace].end(); ++it) {
165
      WorldVector<double> c;
166
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
167
      int nAssoc = it->second.size();
168
#if 0
169
170
171
      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7))
	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
	 it->first, c[0], c[1], (pdb.mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
172
#endif
173
174
175
176
177
178
    }    


    // === 2. check: All periodic DOFs must be symmetric, i.e., if A is mapped ===
    // === to B, then B must be mapped to A.                                   ===

179
180
181
182
183
184
    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);

    if (pdb.mpiRank == 0) {
      for (int i = 1; i < pdb.mpiSize; i++)
	stdMpi.recv(i);
    } else {
185
      stdMpi.send(0, perMap.periodicDofMap[feSpace]);
186
187
    }

188
    stdMpi.startCommunication();
189
190
191
192
193
194

    int foundError = 0;

    // === The boundary DOFs are checked only on the zero rank. === 

    if (pdb.mpiRank == 0) {
195
      // Stores to each rank the periodic DOF mappings of this rank.
196
      map<int, PeriodicDofMap> rankToMaps;
197
      PeriodicDofMap dofMap = perMap.periodicDofMap[feSpace];
198
199
200
201
202
203
204
205
      rankToMaps[0] = dofMap;

      for (int i = 1; i < pdb.mpiSize; i++) {
	PeriodicDofMap &otherMap = stdMpi.getRecvData(i);
	rankToMaps[i] = otherMap;
	
	for (PeriodicDofMap::iterator it = otherMap.begin(); 
	     it != otherMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
206
	  for (DofMap::iterator dofIt = it->second.begin();
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
	       dofIt != it->second.end(); ++dofIt) {
	    if (dofMap.count(it->first) == 1 &&
		dofMap[it->first].count(dofIt->first) == 1) {
	      TEST_EXIT_DBG(dofMap[it->first][dofIt->first] == dofIt->second)
		("Should not happen!\n");
	    } else {
	      dofMap[it->first][dofIt->first] = dofIt->second;
	    }
	  }
	}
      }


      // === Now we test if global DOF A is mapped to B, then B must be mapped ===
      // === to A for the same boundary type.                                  ===

      for (PeriodicDofMap::iterator it = dofMap.begin(); 
	   it != dofMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
225
	for (DofMap::iterator dofIt = it->second.begin();
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
	     dofIt != it->second.end(); ++dofIt) {
	  if (it->second[dofIt->second] != dofIt->first) {
	    MSG("[DBG]  For boundary type %d: DOF %d -> %d, but %d -> %d!\n",
		it ->first, 
		dofIt->first, dofIt->second, 
		dofIt->second, it->second[dofIt->second]);

	    for (int i = 0; i < pdb.mpiSize; i++) {
	      if (rankToMaps[i][it->first].count(dofIt->first) == 1) {
		MSG("[DBG]    %d -> %d in rank %d\n", 
		    dofIt->first, rankToMaps[i][it->first][dofIt->first], i);
	      }

	      if (rankToMaps[i][it->first].count(dofIt->second) == 1) {
		MSG("[DBG]    %d -> %d in rank %d\n", 
		    dofIt->second, rankToMaps[i][it->first][dofIt->second], i);
	      }
	    }
	    
	    ERROR("Wrong periodic DOFs!\n");
	    foundError = 1;
	  }
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
254
255


256
257
258
259
    // === 3. check: On all edge and face periodic DOFs, at least one        ===
    // === componant of coordinates of each periodic DOF pair must be equal  ===
    // === (at least as long we consider periodic boundaries only on         ===
    // === rectangulars and boxes.                                           ===
260
261

    RankToCoords sendCoords;
262
    map<int, vector<BoundaryType> > rankToDofType;
263
264
265
266
267
268
269
270
271
272
273
274

    for (InteriorBoundary::RankToBoundMap::iterator it = pdb.periodicBoundary.boundary.begin();
	 it != pdb.periodicBoundary.boundary.end(); ++it) {
      if (it->first == pdb.mpiRank)
	continue;

      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {	
	if (boundIt->rankObj.subObj == VERTEX)
	  continue;

	DofContainer dofs;
275
	boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs);
276
277
278
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
279
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), feSpace, c);
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
	  sendCoords[it->first].push_back(c);
	  rankToDofType[it->first].push_back(boundIt->type);
	}
      }
    }

    // Each rank must receive exactly the same number of coordinates as it sends
    // to another rank.
    RankToCoords recvCoords;
    for (RankToCoords::iterator it = sendCoords.begin(); 
	 it != sendCoords.end(); ++it)
      recvCoords[it->first].resize(it->second.size());


    StdMpi<CoordsVec> stdMpiCoords(pdb.mpiComm, true);
    stdMpiCoords.send(sendCoords);
    stdMpiCoords.recv(recvCoords);   
297
    stdMpiCoords.startCommunication();
298
299
300
301
302
303
304
305
306
307
308
309
310
311


    for (RankToCoords::iterator it = sendCoords.begin(); 
	 it != sendCoords.end(); ++it) {
      for (unsigned int i = 0; i < it->second.size(); i++) {
	WorldVector<double> &c0 = it->second[i];
	WorldVector<double> &c1 = stdMpiCoords.getRecvData(it->first)[i];


	int nEqual = 0;
	for (int j = 0; j < pdb.mesh->getDim(); j++)
	  if (c0[j] == c1[j])
	    nEqual++;

312
	if (nEqual == 0) {
313
314
315
316
317
318
319
320
321
322
323
324
325
326
	  MSG("[DBG]  %d-ith periodic DOF in boundary between ranks %d <-> %d is not correct!\n",
	      i, pdb.mpiRank, it->first);
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      pdb.mpiRank, c0[0], c0[1], (pdb.mesh->getDim() == 3 ? c0[2] : 0.0));
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      it->first, c1[0], c1[1], (pdb.mesh->getDim() == 3 ? c1[2] : 0.0));

	  foundError = 1;
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Wrond periodic coordinates found on at least on rank!\n");
327
328
329
  }


330
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
331
  {
332
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334

    clock_t first = clock();
335
336
    // Get FE space with basis functions of the highest degree
    const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338

    int testCommonDofs = 1;
339
    Parameters::get("dbg->test common dofs", testCommonDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
342
343
344
345
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

    /// Defines a mapping type from rank numbers to sets of DOFs.
346
    typedef map<int, DofContainer> RankToDofContainer;
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348

    // Maps to each neighbour rank an array of WorldVectors. This array contains the 
349
    // coordinates of all DOFs this rank shares on the interior boundary with the 
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    // neighbour rank. A rank sends the coordinates to another rank, if it owns the
351
    // boundarys DOFs.
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
    RankToCoords sendCoords;

354
355
    // A rank receives all boundary DOFs that are at its interior boundaries but are
    // not owned by the rank. This map stores for each rank the coordinates of DOFs
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
    // this rank expectes to receive from.
    RankToCoords recvCoords;

359
360
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
    pdb.mesh->getDofIndexCoords(feSpace, coords);
361

362
363
364
    for (DofComm::Iterator it(pdb.sendDofs, feSpace); !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	sendCoords[it.getRank()].push_back(coords[it.getDofIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
365

366
367
368
    for (DofComm::Iterator it(pdb.recvDofs, feSpace); !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	recvCoords[it.getRank()].push_back(coords[it.getDofIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
369

370
371
372
    vector<int> sendSize(pdb.mpiSize, 0);
    vector<int> recvSize(pdb.mpiSize, 0);
    vector<int> recvSizeBuffer(pdb.mpiSize, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
375
376
377
378
379
380
381
382
383
384
385
    MPI::Request request[(pdb.mpiSize - 1) * 2];
    int requestCounter = 0;

    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it)
      sendSize[it->first] = it->second.size();

    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it)
      recvSize[it->first] = it->second.size();

    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

386
387
      request[requestCounter++] = 
	pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
390
391
392
393
    }   

    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

394
395
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
398
399
400
    }

    MPI::Request::Waitall(requestCounter, request);


401
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
404
405
406
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
407
408
409
	ERROR("MPI rank %d expectes to receive %d DOFs from rank %d. But this rank sends %d DOFs!\n", 
	      pdb.mpiRank, recvSize[i], i, recvSizeBuffer[i]);	
	foundError = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
      }
    }
412
413
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
418
419
420
421

    // === Now we know that the number of send and received DOFs fits together. ===
    // === So we can check if also the coordinates of the communicated DOFs are ===
    // === the same on both corresponding ranks.                                ===

    StdMpi<CoordsVec> stdMpi(pdb.mpiComm, true);
    stdMpi.send(sendCoords);
    stdMpi.recv(recvCoords);   
422
    stdMpi.startCommunication();
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
425
426
427
428
429

    int dimOfWorld = Global::getGeo(WORLD);

    // === Compare the received with the expected coordinates. ===

    for (RankToCoords::iterator it = stdMpi.getRecvData().begin(); 
	 it != stdMpi.getRecvData().end(); ++it) {
430
431
432
      for (unsigned int i = 0; i < it->second.size(); i++) {
	WorldVector<double> tmp = (it->second)[i];
	tmp -=  recvCoords[it->first][i];
433

434
	if (norm(tmp) > 1e-8) {
435
436
437
	  // === Print error message if the coordinates are not the same. ===
	  if (printCoords) {
	    MSG("[DBG] i = %d\n", i);	  
438
	    stringstream oss;
439
440
	    oss.precision(5);
	    oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
441
		<< " expect coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
442
	    for (int k = 0; k < dimOfWorld; k++) {
443
	      oss << recvCoords[it->first][i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
444
	      if (k + 1 < dimOfWorld)
445
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
446
	    }
447
	    oss << ")  received coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
448
	    for (int k = 0; k < dimOfWorld; k++) {
449
	      oss << (it->second)[i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
450
	      if (k + 1 < dimOfWorld)
451
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
452
	    }
453
454
	    oss << ")";
	    MSG("%s\n", oss.str().c_str());
455
	    
456
	    debug::printInfoByDof(feSpace, 
457
458
				  *(pdb.recvDofs.getDofContainer(it->first, 
								 feSpace)[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
459
	  }
460
461
462
463
	  ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
	  foundError = 1;
	}	 
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
464
    }
465
466
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
469
470
471

    INFO(pdb.info, 8)("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock()));
  }


472
473
474
475
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

476
477
478
479
480
    // Get FE space with basis functions of the highest degree
    const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];

    DOFVector<WorldVector<double> > coords(feSpace, "tmp");
    pdb.mesh->getDofIndexCoords(feSpace, coords);
481

482
    typedef map<WorldVector<double>, int> CoordsIndexMap;
483
484
485
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
486
    for (it.reset(); !it.end(); ++it) {
487
      coordsToIndex[(*it)] = 
488
	pdb.dofFeData[feSpace].mapDofToGlobal[it.getDOFIndex()];    
489

490
491
492
493
//       MSG("   CHECK FOR DOF %d AT COORDS %f %f %f\n",
// 	  coordsToIndex[(*it)], (*it)[0], (*it)[1], (*it)[2]);
    }

494
    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
495
496
497
498
    for (DofComm::Iterator it(pdb.sendDofs, feSpace); !it.end(); it.nextRank())
      stdMpi.send(it.getRank(), coordsToIndex);
    for (DofComm::Iterator it(pdb.recvDofs, feSpace); !it.end(); it.nextRank())
      stdMpi.recv(it.getRank());
499
   
500
    stdMpi.startCommunication();
501
502

    int foundError = 0;
503
504
    for (DofComm::Iterator it(pdb.recvDofs, feSpace); !it.end(); it.nextRank()) {
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it.getRank());
505
506
507
508
509

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
510
	  stringstream oss;
511
512
513
514
515
516
	  oss.precision(5);
	  oss << "DOF at coords ";
	  for (int i = 0; i < Global::getGeo(WORLD); i++)
	    oss << coordsIt->first[i] << " ";
	  oss << " do not fit together on rank " 
	      << pdb.getMpiRank() << " (global index: " 
517
	      << coordsToIndex[coordsIt->first] << ") and on rank "
518
	      << it.getRank() << " (global index: " << coordsIt->second << ")";
519
520
521
522
523
524
525
526
527
528
529

	  MSG("[DBG] %s\n", oss.str().c_str());
	  foundError = 1;
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
  }

530

531
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
532
  {
533
    FUNCNAME("ParallelDebug::testAllElements()");
534

535
    std::set<int> macroElements;
536
537
    int minElementIndex = numeric_limits<int>::max();
    int maxElementIndex = numeric_limits<int>::min();   
538

539
    TraverseStack stack;
540
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
541
    while (elInfo) {
542
543
544
545
546
547
548
549
550
551
      int elIndex = elInfo->getElement()->getIndex();
      minElementIndex = std::min(minElementIndex, elIndex);
      maxElementIndex = std::max(maxElementIndex, elIndex);
      macroElements.insert(elInfo->getElement()->getIndex());
      elInfo = stack.traverseNext(elInfo);      
    }
    
    int globalMinIndex, globalMaxIndex;
    pdb.mpiComm.Allreduce(&minElementIndex, &globalMinIndex, 1, MPI_INT, MPI_MIN);
    pdb.mpiComm.Allreduce(&maxElementIndex, &globalMaxIndex, 1, MPI_INT, MPI_MAX);
552

553
554
555
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
556
557
558
559
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
560
561
562
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
563

564
565
566
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
567
568
569
570
571
      }
    }
  }


572
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
573
574
						    map<int, map<const FiniteElemSpace*, DofContainer> > &sendDofs,
						    map<int, map<const FiniteElemSpace*, DofContainer> > &recvDofs)
575
  {
576
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
577

578
579
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

580
    map<int, int> sendNumber;
581
582
583
584
585
586
587
588
    for (it_type it = sendDofs.begin(); it != sendDofs.end(); ++it)
      for (map<const FiniteElemSpace*, DofContainer>::iterator dcIt = it->second.begin(); dcIt != it->second.end(); ++dcIt)
	sendNumber[it->first] += dcIt->second.size();
    
    map<int, int> recvNumber;
    for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it)
      for (map<const FiniteElemSpace*, DofContainer>::iterator dcIt = it->second.begin(); dcIt != it->second.end(); ++dcIt)
	recvNumber[it->first] += dcIt->second.size();
589
590
591
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
592
    for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it)
593
      stdMpi.recv(it->first);
594
    stdMpi.startCommunication();
595
596
     
    int foundError = 0;
597
    for (map<int, int>::iterator it = stdMpi.getRecvData().begin();
598
	 it != stdMpi.getRecvData().end(); ++it) {
599
      if (it->second != recvNumber[it->first]) {
600
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
601
	      recvNumber[it->first], it->first, it->second);
602
603
604
605
606
607
	foundError = 1;
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
608
609
610
  }


611
612
613
614
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

615
    map<WorldVector<double>, DegreeOfFreedom> cMap;
616
617
618
619
620
    int foundError = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
Thomas Witkowski's avatar
Merge    
Thomas Witkowski committed
621
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
622
623
624
625
626
	WorldVector<double> &c = elInfo->getCoord(i);
	if (cMap.count(c) == 0) {
	  cMap[c] = elInfo->getElement()->getDof(i, 0);
	} else {
	  if (cMap[c] != elInfo->getElement()->getDof(i, 0)) {
627
628
629
	    MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f %f!\n",
		cMap[c], elInfo->getElement()->getDof(i, 0), 
		c[0], c[1], mesh->getDim() == 3 ? c[2] : 0.0);
630
631
632
633
634
635
636
637
638
639
640
641
642
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


643
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
646
647
      const FiniteElemSpace *feSpace = pdb.feSpaces[0];

648
      cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
649
      
Thomas Witkowski's avatar
Thomas Witkowski committed
650
      for (DofMap::iterator it = pdb.dofFeData[feSpace].mapDofToGlobal.begin();
651
	   it != pdb.dofFeData[feSpace].mapDofToGlobal.end(); it++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
652
	DegreeOfFreedom localdof = -1;
653
654
	if (pdb.dofFeData[feSpace].mapLocalToDof.count(it->first) > 0)
	  localdof = pdb.dofFeData[feSpace].mapLocalToDof[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
655
	
656
	cout << "DOF " << it->first << " " 
Thomas Witkowski's avatar
Thomas Witkowski committed
657
		  << it->second << " " 
658
		  << localdof << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
659
	WorldVector<double> coords;
660
	pdb.mesh->getDofIndexCoords(it->first, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
	coords.print();

663
664
665
666
667
668
669
670
671
	for (DofComm::Iterator rit(pdb.sendDofs, feSpace); !rit.end(); rit.nextRank())
	  for (; !rit.endDofIter(); rit.nextDof())
	    if (it->first == rit.getDofIndex())
	      cout << "SEND DOF TO " << rit.getRank() << endl;
	
	for (DofComm::Iterator rit(pdb.recvDofs, feSpace); !rit.end(); rit.nextRank())
	  for (; !rit.endDofIter(); rit.nextDof())
	    if (it->first == rit.getDofIndex())
	      cout << "RECV DOF FROM " << rit.getRank() << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
672

673
	cout << "------" << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
676
677
678
      }
    }
  }


679
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
680
  {
681
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
682

683
684
685
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
686
687
    PeriodicMap &perMap = pdb.getPeriodicMap();

688
689
    const FiniteElemSpace* feSpace = pdb.feSpaces[0];

Thomas Witkowski's avatar
Thomas Witkowski committed
690
    typedef map<DegreeOfFreedom, DegreeOfFreedom> DofMap;
691
    typedef map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693

    if (rank == -1 || pdb.mpiRank == rank) {
694
      cout << "====== DOF MAP PERIODIC ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
695

696
697
      for (PeriodicDofMap::iterator it = perMap.periodicDofMap.begin();
	   it != perMap.periodicDofMap.end(); ++it) {
698
	cout << "DOF MAP " << it->first << ": ";
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
701
702
	  cout << *dofit << " ";
	cout << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
703
704

	DegreeOfFreedom localdof = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
705
	for (DofMap::iterator dofIt = pdb.dofFeData[feSpace].mapDofToGlobal.begin();
706
	     dofIt != pdb.dofFeData[feSpace].mapDofToGlobal.end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
710
711
712
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

	TEST_EXIT(localdof != -1)("There is something wrong!\n");

	WorldVector<double> coords;
713
	pdb.mesh->getDofIndexCoords(localdof, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
716
	coords.print();
      }
    }
717
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
  }

  
721
722
723
724
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
  {
    if (rank == -1 || pdb.mpiRank == rank) {
727
      cout << "====== RANK DOF INFORMATION ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
728

729
      cout << "  RANK OWNED DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
730
731
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
732
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
733
	WorldVector<double> coords;
734
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
737
	coords.print();
      }

738
      cout << "  RANK ALL DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
741
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
742
	WorldVector<double> coords;
743
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
746
747
748
	coords.print();
      }      
    }
  }

749

750
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
751
  {
752
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
753

Thomas Witkowski's avatar
Thomas Witkowski committed
754
    int tmp = 0;
755
    Parameters::get("parallel->debug->print boundary info", tmp);
756
    if (tmp <= 0)
757
758
      return;

759
    for (InteriorBoundary::iterator it(pdb.rankIntBoundary); !it.end(); ++it) {
760
761
762
763
764
765
766
767
768
769
770
771
772
773
      MSG("Rank owned boundary with rank %d: \n", it.getRank());
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it) {
      MSG("Other owned boundary with rank %d: \n", it.getRank());
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }
774
775

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
776
777
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
778
779
780
781
782
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }    
783
784
  }

785

786
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
787
				     string prefix, string postfix)
788
  {
789
    FUNCNAME("ParallelDebug::writeCoordsFile()");
790

791
    const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];
792

793
    stringstream filename;
794
795
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

796
797
    DOFVector<WorldVector<double> > coords(feSpace, "tmp");
    pdb.mesh->getDofIndexCoords(feSpace, coords);
798

799
    typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
800
801
    ElDofMap elDofMap;
    TraverseStack stack;
802
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
803
    vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
804
805
806
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
807
				 feSpace->getAdmin(), localIndices);
808
809
810
811
812
813
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

    // === Write informations about all DOFs. ===

814
    ofstream file;
815
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
816
817
//     file << "# First line contains number of DOFs, than each line has the format\n";
//     file << "# Local DOF index     Global DOF index     Is rank DOF     x-coord     y-coord     z-coord\n";
818
819
820
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
821
      file << it.getDOFIndex() << " " 
822
	   << pdb.dofFeData[feSpace].mapDofToGlobal[it.getDOFIndex()] << " "
823
	   << pdb.getIsRankDof(feSpace, it.getDOFIndex());
824
825
826
827
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
828
829
830

    // === Write to all elements in ranks mesh the included dofs. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
833
//     file << "\n\n";
//     file << "# First line containes number of elements in mesh, second line contain the number of DOFs per element.\n";
//     file << "# Than, each entry contains of two lines. The first is the element index, the second line is a list with the local DOF indices of this element.\n";
834
835
836
837
838
839
840
841
842
    file << elDofMap.size() << "\n";
    file << basisFcts->getNumber() << "\n";
    for (ElDofMap::iterator it = elDofMap.begin(); it != elDofMap.end(); ++it) {
      file << it->first << "\n";
      for (int i = 0; i < basisFcts->getNumber(); i++) 
	file << it->second[i] << " ";
      file << "\n";
    }

843
844
845
    file.close();
  }

846
847
848
849
850
851
852
853
854
855
856
857

  void ParallelDebug::writePartitioning(MeshDistributor &pdb, string filename)
  {
    FUNCNAME("ParallelDebug::writeParitioning()");

    map<int, double> vec;    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    
    while (elInfo) {		  
      int index = elInfo->getElement()->getIndex();
858
      vec[index] = pdb.partitionMap[index];
859
860
861
862
863
864
      elInfo = stack.traverseNext(elInfo);
    }

    ElementFileWriter::writeFile(vec, pdb.mesh, filename);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
865

866
  void ParallelDebug::writePartitioningFile(string filename, 
867
					    int counter,
868
					    const FiniteElemSpace *feSpace)
869
870
871
872
873
874
875
876
877
  {
    FUNCNAME("ParallelDebug::writePartitioningFile()");

    stringstream oss;
    oss << filename;
    if (counter >= 0)
      oss << "-" << counter;
    oss << ".vtu";

878
879
880
    DOFVector<double> tmp(feSpace, "tmp");
    tmp.set(MPI::COMM_WORLD.Get_rank());
    VtkWriter::writeFile(&tmp, oss.str());
881
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903

  
  bool ParallelDebug::followThisBound(int rankElIndex, int neighElIndex)
  {
    FUNCNAME("ParallelDebug::followThisBound()");

    int el0 = std::min(rankElIndex, neighElIndex);
    int el1 = std::max(rankElIndex, neighElIndex);

    vector<int> els;
    Parameters::get("parallel->debug->follow boundary", els);
    if (els.size() != 2)
      return false;

    return (el0 == els[0] && el1 == els[1]);
  }


  void ParallelDebug::followBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::followBoundary()");

904
    for (InteriorBoundary::iterator it(pdb.rankIntBoundary); !it.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
905
906
907
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank, 
				     it->rankObj.elIndex, 
908
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
909
910
911
912
913

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it)
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank,
				     it->rankObj.elIndex,
914
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
  }


  void ParallelDebug::followBoundary(Mesh *mesh, 
				     AtomicBoundary &bound, 
				     MeshStructure &code)
  {
    FUNCNAME("ParallelDebug::followBoundary()");

    if (mesh->getDim() != bound.rankObj.subObj)
      return;

    if (!followThisBound(bound.rankObj.elIndex, bound.neighObj.elIndex))
      return;

    MSG("Mesh structure code of bound %d/%d <-> %d/%d: %s\n",
	bound.rankObj.elIndex, mesh->getDim(),
	bound.neighObj.elIndex, mesh->getDim(), 
	code.toStr().c_str());
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
935
}