ParallelDebug.cc 30.4 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
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
    for (map<int, std::set<BoundaryType> >::iterator it = pdb.periodicDofAssociations[feSpace].begin();
	 it != pdb.periodicDofAssociations[feSpace].end(); ++it) {
163
      WorldVector<double> c;
164
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
165
166
167
168
169
170
171
172
173
174
      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.                                   ===

175
176
177
178
179
180
    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);

    if (pdb.mpiRank == 0) {
      for (int i = 1; i < pdb.mpiSize; i++)
	stdMpi.recv(i);
    } else {
181
      stdMpi.send(0, pdb.periodicDofMap[feSpace]);
182
183
    }

184
    stdMpi.startCommunication();
185
186
187
188
189
190

    int foundError = 0;

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

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


252
253
254
255
    // === 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.                                           ===
256
257

    RankToCoords sendCoords;
258
    map<int, vector<BoundaryType> > rankToDofType;
259
260
261
262
263
264
265
266
267
268
269
270

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


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

308
	if (nEqual == 0) {
309
310
311
312
313
314
315
316
317
318
319
320
321
322
	  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");
323
324
325
  }


326
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
327
  {
328
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330

    clock_t first = clock();
331
332
    // 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
333
334

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

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

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

350
351
    // 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
352
353
354
    // this rank expectes to receive from.
    RankToCoords recvCoords;

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

358
359
360
    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
361

362
363
364
    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
365

366
367
368
    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
369
370
371
372
373
374
375
376
377
378
379
380
381
    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;

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

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

390
391
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
    }

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


397
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

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

    // === 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);   
418
    stdMpi.startCommunication();
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
421
422
423
424
425

    int dimOfWorld = Global::getGeo(WORLD);

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

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

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

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


467
468
469
470
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

471
472
473
474
475
    // 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);
476

477
    typedef map<WorldVector<double>, int> CoordsIndexMap;
478
479
480
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
481
    for (it.reset(); !it.end(); ++it) {
482
483
      coordsToIndex[(*it)] = 
	pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()];    
484

485
486
487
488
//       MSG("   CHECK FOR DOF %d AT COORDS %f %f %f\n",
// 	  coordsToIndex[(*it)], (*it)[0], (*it)[1], (*it)[2]);
    }

489
    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
490
491
492
493
    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());
494
   
495
    stdMpi.startCommunication();
496
497

    int foundError = 0;
498
499
    for (DofComm::Iterator it(pdb.recvDofs, feSpace); !it.end(); it.nextRank()) {
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it.getRank());
500
501
502
503
504

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
505
	  stringstream oss;
506
507
508
509
510
511
	  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: " 
512
	      << coordsToIndex[coordsIt->first] << ") and on rank "
513
	      << it.getRank() << " (global index: " << coordsIt->second << ")";
514
515
516
517
518
519
520
521
522
523
524

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

525

526
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
527
  {
528
    FUNCNAME("ParallelDebug::testAllElements()");
529

530
    std::set<int> macroElements;
531
532
    int minElementIndex = numeric_limits<int>::max();
    int maxElementIndex = numeric_limits<int>::min();   
533

534
    TraverseStack stack;
535
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
536
    while (elInfo) {
537
538
539
540
541
542
543
544
545
546
      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);
547

548
549
550
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
551
552
553
554
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
555
556
557
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
558

559
560
561
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
562
563
564
565
566
      }
    }
  }


567
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
568
569
						    map<int, map<const FiniteElemSpace*, DofContainer> > &sendDofs,
						    map<int, map<const FiniteElemSpace*, DofContainer> > &recvDofs)
570
  {
571
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
572

573
574
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

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

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


606
607
608
609
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

610
    map<WorldVector<double>, DegreeOfFreedom> cMap;
611
612
613
614
615
    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
616
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
617
618
619
620
621
	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)) {
622
623
624
	    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);
625
626
627
628
629
630
631
632
633
634
635
636
637
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


638
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
641
642
      const FiniteElemSpace *feSpace = pdb.feSpaces[0];

643
      cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
644
      
645
646
      for (DofMapping::iterator it = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin();
	   it != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); it++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
647
	DegreeOfFreedom localdof = -1;
648
649
	if (pdb.dofFeData[feSpace].mapLocalDofIndex.count(it->first) > 0)
	  localdof = pdb.dofFeData[feSpace].mapLocalDofIndex[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
650
	
651
	cout << "DOF " << it->first << " " 
Thomas Witkowski's avatar
Thomas Witkowski committed
652
		  << it->second << " " 
653
		  << localdof << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
654
	WorldVector<double> coords;
655
	pdb.mesh->getDofIndexCoords(it->first, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
	coords.print();

658
659
660
661
662
663
664
665
666
	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
667

668
	cout << "------" << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
671
672
673
      }
    }
  }


674
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
675
  {
676
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
677

678
679
680
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
681
682
    const FiniteElemSpace* feSpace = pdb.feSpaces[0];

683
684
    typedef map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
685
686

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

689
690
      for (PeriodicDofMap::iterator it = pdb.periodicDofMap.begin();
	   it != pdb.periodicDofMap.end(); ++it) {
691
	cout << "DOF MAP " << it->first << ": ";
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
694
695
	  cout << *dofit << " ";
	cout << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
696
697

	DegreeOfFreedom localdof = -1;
698
699
	for (DofMapping::iterator dofIt = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin();
	     dofIt != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
702
703
704
705
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

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

	WorldVector<double> coords;
706
	pdb.mesh->getDofIndexCoords(localdof, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
	coords.print();
      }
    }
710
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
713
  }

  
714
715
716
717
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
  {
    if (rank == -1 || pdb.mpiRank == rank) {
720
      cout << "====== RANK DOF INFORMATION ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
721

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

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

742

743
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
744
  {
745
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
746

Thomas Witkowski's avatar
Thomas Witkowski committed
747
    int tmp = 0;
748
    Parameters::get("parallel->debug->print boundary info", tmp);
749
    if (tmp <= 0)
750
751
      return;

752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
    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);
    }
767
768

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
769
770
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
771
772
773
774
775
      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);
    }    
776
777
  }

778

779
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
780
				     string prefix, string postfix)
781
  {
782
    FUNCNAME("ParallelDebug::writeCoordsFile()");
783

784
785
    const FiniteElemSpace *feSpace = pdb.feSpaces[0];

786
    stringstream filename;
787
788
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

792
    typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
793
794
    ElDofMap elDofMap;
    TraverseStack stack;
795
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
796
    vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
797
798
799
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
800
				 feSpace->getAdmin(), localIndices);
801
802
803
804
805
806
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

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

807
    ofstream file;
808
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
//     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";
811
812
813
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
814
      file << it.getDOFIndex() << " " 
815
816
	   << pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()] << " "
	   << pdb.getIsRankDof(feSpace, it.getDOFIndex());
817
818
819
820
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
821
822
823

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

Thomas Witkowski's avatar
Thomas Witkowski committed
824
825
826
//     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";
827
828
829
830
831
832
833
834
835
    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";
    }

836
837
838
    file.close();
  }

839
840
841
842
843
844
845
846
847
848
849
850

  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();
851
      vec[index] = pdb.partitionMap[index];
852
853
854
855
856
857
      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
858

859
  void ParallelDebug::writePartitioningFile(string filename, 
860
					    int counter,
861
					    const FiniteElemSpace *feSpace)
862
863
864
865
866
867
868
869
870
871
872
873
874
  {
    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
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900

  
  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, 
901
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
905
906

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it)
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank,
				     it->rankObj.elIndex,
907
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
  }


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