ParallelDebug.cc 30.1 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.myIntBoundary.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. ===

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

      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
    return;

150
151
    TEST_EXIT(pdb.feSpaces.size() == 1)("Shoudl not happen!\n");

152
153
    // === 1. check: All periodic DOFs should have at least a correct number ===
    // === of periodic associations.                                         ===
154
   
155
156
157
    for (map<int, std::set<BoundaryType> >::iterator it = pdb.periodicDofAssociations.begin();
	 it != pdb.periodicDofAssociations.end(); ++it) {
      WorldVector<double> c;
158
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
159
160
161
162
163
164
165
166
167
168
      int nAssoc = it->second.size();
      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);
    }    


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

169
170
171
172
173
174
175
176
177
    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);

    if (pdb.mpiRank == 0) {
      for (int i = 1; i < pdb.mpiSize; i++)
	stdMpi.recv(i);
    } else {
      stdMpi.send(0, pdb.periodicDof);
    }

178
    stdMpi.startCommunication();
179
180
181
182
183
184

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
185
      // Stores to each rank the periodic DOF mappings of this rank.
186
      map<int, PeriodicDofMap> rankToMaps;
187
188
189
190
191
192
193
194
195
      PeriodicDofMap dofMap = pdb.periodicDof;
      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) {
196
	  for (DofMapping::iterator dofIt = it->second.begin();
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
	       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) {
215
	for (DofMapping::iterator dofIt = it->second.begin();
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
	     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");
244
245
246
247
248
249
250


    // === 3. check: On all edge and face periodic DOFs, at least on coordinate of ===
    // === each periodic DOF pair must be equal (at least as long we consider      ===
    // === periodic boundaries only on rectangulars and boxes.                     ===

    RankToCoords sendCoords;
251
    map<int, vector<BoundaryType> > rankToDofType;
252
253
254
255
256
257
258
259
260
261
262
263

    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;
264
	boundIt->rankObj.el->getAllDofs(pdb.feSpaces[0], boundIt->rankObj, dofs);
265
266
267
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
268
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpaces[0], c);
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
	  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);   
286
    stdMpiCoords.startCommunication();
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316


    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++;

	if ((rankToDofType[it->first][i] >= -3 && nEqual < 2) ||
	    (rankToDofType[it->first][i] < -3 && nEqual == 0)) {
	  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");
317
318
319
  }


320
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
321
  {
322
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324

    clock_t first = clock();
325
326
    // 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
327
328

    int testCommonDofs = 1;
329
    Parameters::get("dbg->test common dofs", testCommonDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
332
333
334
335
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

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

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

344
345
    // 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
346
347
348
    // this rank expectes to receive from.
    RankToCoords recvCoords;

349
350
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
    pdb.mesh->getDofIndexCoords(feSpace, coords);
351

352
353
354
    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
355

356
357
358
    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
359

360
361
362
    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
363
364
365
366
367
368
369
370
371
372
373
374
375
    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;

376
377
      request[requestCounter++] = 
	pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
381
382
383
    }   

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

384
385
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389
390
    }

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


391
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
397
398
399
	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
400
401
      }
    }
402
403
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
406
407
408
409
410
411

    // === 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);   
412
    stdMpi.startCommunication();
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
415
416
417
418
419

    int dimOfWorld = Global::getGeo(WORLD);

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

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

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

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


461
462
463
464
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

465
466
467
468
469
    // 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);
470

471
    typedef map<WorldVector<double>, int> CoordsIndexMap;
472
473
474
475
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it)
476
477
      coordsToIndex[(*it)] = 
	pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()];    
478
479

    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
480
481
482
483
    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());
484
   
485
    stdMpi.startCommunication();
486
487

    int foundError = 0;
488
489
    for (DofComm::Iterator it(pdb.recvDofs, feSpace); !it.end(); it.nextRank()) {
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it.getRank());
490
491
492
493
494

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
495
	  stringstream oss;
496
497
498
499
500
501
502
	  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: " 
	      << coordsToIndex[coordsIt->first] << " and on rank "
503
	      << it.getRank() << " (global index: " << coordsIt->second << ")";
504
505
506
507
508
509
510
511
512
513
514

	  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");
  }

515

516
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
517
  {
518
    FUNCNAME("ParallelDebug::testAllElements()");
519

520
    std::set<int> macroElements;
521
522
    int minElementIndex = numeric_limits<int>::max();
    int maxElementIndex = numeric_limits<int>::min();   
523

524
    TraverseStack stack;
525
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
526
    while (elInfo) {
527
528
529
530
531
532
533
534
535
536
      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);
537

538
539
540
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
541
542
543
544
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
545
546
547
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
548

549
550
551
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
552
553
554
555
556
      }
    }
  }


557
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
558
559
						    map<int, map<const FiniteElemSpace*, DofContainer> > &sendDofs,
						    map<int, map<const FiniteElemSpace*, DofContainer> > &recvDofs)
560
  {
561
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
562

563
564
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

565
    map<int, int> sendNumber;
566
567
568
569
570
571
572
573
    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();
574
575
576
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
577
    for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it)
578
      stdMpi.recv(it->first);
579
    stdMpi.startCommunication();
580
581
     
    int foundError = 0;
582
    for (map<int, int>::iterator it = stdMpi.getRecvData().begin();
583
	 it != stdMpi.getRecvData().end(); ++it) {
584
      if (it->second != recvNumber[it->first]) {
585
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
586
	      recvNumber[it->first], it->first, it->second);
587
588
589
590
591
592
	foundError = 1;
      }
    }

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


596
597
598
599
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

600
    map<WorldVector<double>, DegreeOfFreedom> cMap;
601
602
603
604
605
    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
606
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
607
608
609
610
611
	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)) {
612
613
614
	    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);
615
616
617
618
619
620
621
622
623
624
625
626
627
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


628
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
631
632
      const FiniteElemSpace *feSpace = pdb.feSpaces[0];

633
      cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
634
      
635
636
      for (DofMapping::iterator it = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin();
	   it != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); it++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
637
	DegreeOfFreedom localdof = -1;
638
639
	if (pdb.dofFeData[feSpace].mapLocalDofIndex.count(it->first) > 0)
	  localdof = pdb.dofFeData[feSpace].mapLocalDofIndex[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
640
	
641
	cout << "DOF " << it->first << " " 
Thomas Witkowski's avatar
Thomas Witkowski committed
642
		  << it->second << " " 
643
		  << localdof << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
644
	WorldVector<double> coords;
645
	pdb.mesh->getDofIndexCoords(it->first, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
646
647
	coords.print();

648
649
650
651
652
653
654
655
656
	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
657

658
	cout << "------" << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660
661
662
663
      }
    }
  }


664
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
665
  {
666
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
667

668
669
670
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
671
672
    const FiniteElemSpace* feSpace = pdb.feSpaces[0];

673
674
    typedef map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676

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

      for (PeriodicDofMap::iterator it = pdb.periodicDof.begin();
	   it != pdb.periodicDof.end(); ++it) {
681
	cout << "DOF MAP " << it->first << ": ";
Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
684
685
	  cout << *dofit << " ";
	cout << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687

	DegreeOfFreedom localdof = -1;
688
689
	for (DofMapping::iterator dofIt = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin();
	     dofIt != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
692
693
694
695
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

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

	WorldVector<double> coords;
696
	pdb.mesh->getDofIndexCoords(localdof, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
699
	coords.print();
      }
    }
700
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
703
  }

  
704
705
706
707
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
708
709
  {
    if (rank == -1 || pdb.mpiRank == rank) {
710
      cout << "====== RANK DOF INFORMATION ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
711

712
      cout << "  RANK OWNED DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
713
714
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
715
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
716
	WorldVector<double> coords;
717
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
	coords.print();
      }

721
      cout << "  RANK ALL DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
724
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
725
	WorldVector<double> coords;
726
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
729
730
731
	coords.print();
      }      
    }
  }

732

733
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
734
  {
735
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
736

Thomas Witkowski's avatar
Thomas Witkowski committed
737
    int tmp = 0;
738
    Parameters::get("parallel->debug->print boundary info", tmp);
739
    if (tmp <= 0)
740
741
      return;

742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
    for (InteriorBoundary::iterator it(pdb.myIntBoundary); !it.end(); ++it) {
      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);
    }
757
758

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
759
760
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
761
762
763
764
765
      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);
    }    
766
767
  }

768

769
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
770
				     string prefix, string postfix)
771
  {
772
    FUNCNAME("ParallelDebug::writeCoordsFile()");
773

774
775
    const FiniteElemSpace *feSpace = pdb.feSpaces[0];

776
    stringstream filename;
777
778
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

782
    typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
783
784
    ElDofMap elDofMap;
    TraverseStack stack;
785
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
786
    vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
787
788
789
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
790
				 feSpace->getAdmin(), localIndices);
791
792
793
794
795
796
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

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

797
    ofstream file;
798
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
//     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";
801
802
803
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
804
      file << it.getDOFIndex() << " " 
805
806
	   << pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()] << " "
	   << pdb.getIsRankDof(feSpace, it.getDOFIndex());
807
808
809
810
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
811
812
813

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

Thomas Witkowski's avatar
Thomas Witkowski committed
814
815
816
//     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";
817
818
819
820
821
822
823
824
825
    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";
    }

826
827
828
    file.close();
  }

829
830
831
832
833
834
835
836
837
838
839
840

  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();
841
      vec[index] = pdb.partitionMap[index];
842
843
844
845
846
847
      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
848

849
  void ParallelDebug::writePartitioningFile(string filename, 
850
					    int counter,
851
					    const FiniteElemSpace *feSpace)
852
853
854
855
856
857
858
859
860
861
862
863
864
  {
    FUNCNAME("ParallelDebug::writePartitioningFile()");

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

    DOFVector<double> tmpa(feSpace, "tmp");
    tmpa.set(MPI::COMM_WORLD.Get_rank());
    VtkWriter::writeFile(&tmpa, oss.str());
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890

  
  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()");

    for (InteriorBoundary::iterator it(pdb.myIntBoundary); !it.end(); ++it)
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank, 
				     it->rankObj.elIndex, 
891
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
892
893
894
895
896

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it)
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank,
				     it->rankObj.elIndex,
897
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
  }


  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
918
}