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


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

namespace AMDiS {

25
26
27
  using namespace std;


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

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

    MPI::Request request[pdb.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
163
164

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

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

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

186
    stdMpi.startCommunication();
187
188
189
190
191
192

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
193
      // Stores to each rank the periodic DOF mappings of this rank.
194
      map<int, PeriodicDofMap> rankToMaps;
195
      PeriodicDofMap dofMap = perMap.periodicDofMap[feSpace];
196
197
198
199
200
201
202
203
      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) {
204
	  for (DofMapping::iterator dofIt = it->second.begin();
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
	       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) {
223
	for (DofMapping::iterator dofIt = it->second.begin();
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
250
251
	     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");
252
253


254
255
256
257
    // === 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.                                           ===
258
259

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

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


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

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


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

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

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

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

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

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

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

360
361
362
    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
363

364
365
366
    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
367

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

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

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

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

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


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

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

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

    int dimOfWorld = Global::getGeo(WORLD);

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

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

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

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


469
470
471
472
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

473
474
475
476
477
    // 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);
478

479
    typedef map<WorldVector<double>, int> CoordsIndexMap;
480
481
482
    CoordsIndexMap coordsToIndex;

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

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

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

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

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

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

527

528
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
529
  {
530
    FUNCNAME("ParallelDebug::testAllElements()");
531

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

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

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

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

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


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

575
576
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

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

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


608
609
610
611
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

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


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

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

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

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


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

680
681
682
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
683
684
    PeriodicMap &perMap = pdb.getPeriodicMap();

685
686
    const FiniteElemSpace* feSpace = pdb.feSpaces[0];

687
688
    typedef map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
689
690

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

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

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

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

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

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

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

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

746

747
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
748
  {
749
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
750

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

756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
    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);
    }
771
772

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

782

783
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
784
				     string prefix, string postfix)
785
  {
786
    FUNCNAME("ParallelDebug::writeCoordsFile()");
787

788
789
    const FiniteElemSpace *feSpace = pdb.feSpaces[0];

790
    stringstream filename;
791
792
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
830
//     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";
831
832
833
834
835
836
837
838
839
    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";
    }

840
841
842
    file.close();
  }

843
844
845
846
847
848
849
850
851
852
853
854

  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();
855
      vec[index] = pdb.partitionMap[index];
856
857
858
859
860
861
      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
862

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

  
  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, 
905
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
906
907
908
909
910

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


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