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

167
168
169
170
171
172
173
174
175
    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);
    }

176
    stdMpi.startCommunication();
177
178
179
180
181
182

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
183
      // Stores to each rank the periodic DOF mappings of this rank.
184
      map<int, PeriodicDofMap> rankToMaps;
185
186
187
188
189
190
191
192
193
      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) {
194
	  for (DofMapping::iterator dofIt = it->second.begin();
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
	       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) {
213
	for (DofMapping::iterator dofIt = it->second.begin();
214
215
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
	     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");
242
243
244
245
246
247
248


    // === 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;
249
    map<int, vector<BoundaryType> > rankToDofType;
250
251
252
253
254
255
256
257
258
259
260
261

    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;
262
	boundIt->rankObj.el->getAllDofs(pdb.feSpace, boundIt->rankObj, dofs);
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpace, c);
	  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);   
284
    stdMpiCoords.startCommunication();
285
286
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


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


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

    clock_t first = clock();

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

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

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

340
341
    // 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
342
343
344
345
346
    // this rank expectes to receive from.
    RankToCoords recvCoords;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
350
351
352
353
354
355
356
357
358
359
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin(); 
	   dofIt != it->second.end(); ++dofIt)
	sendCoords[it->first].push_back(coords[**dofIt]);

    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	recvCoords[it->first].push_back(coords[**dofIt]);

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
425
426
427
	if (norm(tmp) > 1e-13) {
	  // === 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
	    
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
	    debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
	  }
448
449
450
451
	  ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
	  foundError = 1;
	}	 
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
452
    }
453
454
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
457
458
459

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


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

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

467
    typedef map<WorldVector<double>, int> CoordsIndexMap;
468
469
470
471
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it)
472
      coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()];    
473
474
475
476
477
478
479
480
481

    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      stdMpi.send(it->first, coordsToIndex);
    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it)
      stdMpi.recv(it->first);
   
482
    stdMpi.startCommunication();
483
484
485
486
487
488
489
490
491
492

    int foundError = 0;
    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it) {
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it->first);

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

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

513

514
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
515
  {
516
    FUNCNAME("ParallelDebug::testAllElements()");
517

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

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

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

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

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


555
556
557
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
						    RankToDofContainer &sendDofs,
						    RankToDofContainer &recvDofs)
558
  {
559
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
560

561
    map<int, int> sendNumber;
562
563
564
565
566
567
568
    for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it)
      sendNumber[it->first] = it->second.size();
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
    for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it)    
      stdMpi.recv(it->first);
569
    stdMpi.startCommunication();
570
571
     
    int foundError = 0;
572
    for (map<int, int>::iterator it = stdMpi.getRecvData().begin();
573
574
575
576
577
578
579
580
581
582
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second != static_cast<int>(recvDofs[it->first].size())) {
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
	      recvDofs[it->first].size(), it->first, it->second);
	foundError = 1;
      }
    }

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


586
587
588
589
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

590
    map<WorldVector<double>, DegreeOfFreedom> cMap;
591
592
593
594
595
    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
596
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
597
598
599
600
601
	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)) {
602
603
604
	    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);
605
606
607
608
609
610
611
612
613
614
615
616
617
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


618
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
619
620
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
621
      cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
622
623
624
625
      
      for (DofMapping::iterator it = pdb.mapLocalGlobalDofs.begin();
	   it != pdb.mapLocalGlobalDofs.end(); it++) {
	DegreeOfFreedom localdof = -1;
626
627
	if (pdb.mapLocalDofIndex.count(it->first) > 0)
	  localdof = pdb.mapLocalDofIndex[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
628
	
629
	cout << "DOF " << it->first << " " 
Thomas Witkowski's avatar
Thomas Witkowski committed
630
		  << it->second << " " 
631
		  << localdof << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
635
636
637
638
639
640
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, coords);
	coords.print();

	for (RankToDofContainer::iterator rankit = pdb.sendDofs.begin();
	     rankit != pdb.sendDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
641
	      cout << "SEND DOF TO " << rankit->first << endl;	  
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
644
645
646
647
648
	}

	for (RankToDofContainer::iterator rankit = pdb.recvDofs.begin();
	     rankit != pdb.recvDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
649
	      cout << "RECV DOF FROM " << rankit->first << endl;	  
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
	}

652
	cout << "------" << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
655
656
657
      }
    }
  }


658
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
659
  {
660
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
661

662
663
664
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
665
666
    typedef map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668

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

      for (PeriodicDofMap::iterator it = pdb.periodicDof.begin();
	   it != pdb.periodicDof.end(); ++it) {
673
	cout << "DOF MAP " << it->first << ": ";
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
676
677
	  cout << *dofit << " ";
	cout << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
680
681
682
683
684
685
686
687
688
689
690
691

	DegreeOfFreedom localdof = -1;
	for (DofMapping::iterator dofIt = pdb.mapLocalGlobalDofs.begin();
	     dofIt != pdb.mapLocalGlobalDofs.end(); ++dofIt)
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

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

	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(localdof, pdb.feSpace, coords);
	coords.print();
      }
    }
692
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
693
694
695
  }

  
696
697
698
699
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
  {
    if (rank == -1 || pdb.mpiRank == rank) {
702
      cout << "====== RANK DOF INFORMATION ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
703

704
      cout << "  RANK OWNED DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
707
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
708
709
710
711
712
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }

713
      cout << "  RANK ALL DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
716
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
719
720
721
722
723
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }      
    }
  }

724

725
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
726
  {
727
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
728

Thomas Witkowski's avatar
Thomas Witkowski committed
729
    int tmp = 0;
730
    Parameters::get("parallel->debug->print boundary info", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
731
    if (tmp <= 1)
732
733
      return;

734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
    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);
    }
749
750

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
751
752
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
753
754
755
756
757
      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);
    }    
758
759
  }

760

761
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
762
				     string prefix, string postfix)
763
  {
764
    FUNCNAME("ParallelDebug::writeCoordsFile()");
765

766
    stringstream filename;
767
768
769
770
771
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

772
    typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
773
774
775
    ElDofMap elDofMap;
    TraverseStack stack;
    const BasisFunction *basisFcts = pdb.feSpace->getBasisFcts();
776
    vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
777
778
779
780
781
782
783
784
785
786
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
				 pdb.feSpace->getAdmin(), localIndices);
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

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

787
    ofstream file;
788
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
    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";
791
792
793
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
794
795
796
      file << it.getDOFIndex() << " " 
	   << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " "
	   << pdb.getIsRankDof(it.getDOFIndex());
797
798
799
800
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
801
802
803

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

804
805
806
    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";
807
808
809
810
811
812
813
814
815
    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";
    }

816
817
818
    file.close();
  }

819
820
821
822
823
824
825
826
827
828
829
830

  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();
831
      vec[index] = pdb.partitionMap[index];
832
833
834
835
836
837
      elInfo = stack.traverseNext(elInfo);
    }

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

838
  void ParallelDebug::writePartitioningFile(string filename, 
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
					    int counter,
					    FiniteElemSpace *feSpace)
  {
    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
854
}