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


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

namespace AMDiS {

25
26
27
  using namespace std;


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

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

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

39
40
41

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

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

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

56
57
58

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

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

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


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

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

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

      sendBuffers.push_back(sBuffer);

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

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

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

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

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

100
101
102
103
104

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

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

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

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

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

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


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

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

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

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


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

148
149
    return;

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

152
153
    // === 1. check: All periodic DOFs should have at least a correct number ===
    // === of periodic associations.                                         ===
154
   
155
156
157
    for (map<int, std::set<BoundaryType> >::iterator it = pdb.periodicDofAssociations.begin();
	 it != pdb.periodicDofAssociations.end(); ++it) {
      WorldVector<double> c;
158
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
159
160
161
162
163
164
165
166
167
168
      int nAssoc = it->second.size();
      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7))
	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
	 it->first, c[0], c[1], (pdb.mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
    }    


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

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

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

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

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
185
      // Stores to each rank the periodic DOF mappings of this rank.
186
      map<int, PeriodicDofMap> rankToMaps;
187
188
189
190
191
192
193
194
195
      PeriodicDofMap dofMap = pdb.periodicDof;
      rankToMaps[0] = dofMap;

      for (int i = 1; i < pdb.mpiSize; i++) {
	PeriodicDofMap &otherMap = stdMpi.getRecvData(i);
	rankToMaps[i] = otherMap;
	
	for (PeriodicDofMap::iterator it = otherMap.begin(); 
	     it != otherMap.end(); ++it) {
196
	  for (DofMapping::iterator dofIt = it->second.begin();
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
	       dofIt != it->second.end(); ++dofIt) {
	    if (dofMap.count(it->first) == 1 &&
		dofMap[it->first].count(dofIt->first) == 1) {
	      TEST_EXIT_DBG(dofMap[it->first][dofIt->first] == dofIt->second)
		("Should not happen!\n");
	    } else {
	      dofMap[it->first][dofIt->first] = dofIt->second;
	    }
	  }
	}
      }


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

      for (PeriodicDofMap::iterator it = dofMap.begin(); 
	   it != dofMap.end(); ++it) {
215
	for (DofMapping::iterator dofIt = it->second.begin();
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
	     dofIt != it->second.end(); ++dofIt) {
	  if (it->second[dofIt->second] != dofIt->first) {
	    MSG("[DBG]  For boundary type %d: DOF %d -> %d, but %d -> %d!\n",
		it ->first, 
		dofIt->first, dofIt->second, 
		dofIt->second, it->second[dofIt->second]);

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

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

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
244
245
246
247
248
249
250


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

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

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

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

	DofContainer dofs;
264
	boundIt->rankObj.el->getAllDofs(pdb.feSpaces[0], boundIt->rankObj, dofs);
265
266
267
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
268
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpaces[0], c);
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
	  sendCoords[it->first].push_back(c);
	  rankToDofType[it->first].push_back(boundIt->type);
	}
      }
    }

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


    StdMpi<CoordsVec> stdMpiCoords(pdb.mpiComm, true);
    stdMpiCoords.send(sendCoords);
    stdMpiCoords.recv(recvCoords);   
286
    stdMpiCoords.startCommunication();
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316


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


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

	if ((rankToDofType[it->first][i] >= -3 && nEqual < 2) ||
	    (rankToDofType[it->first][i] < -3 && nEqual == 0)) {
	  MSG("[DBG]  %d-ith periodic DOF in boundary between ranks %d <-> %d is not correct!\n",
	      i, pdb.mpiRank, it->first);
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      pdb.mpiRank, c0[0], c0[1], (pdb.mesh->getDim() == 3 ? c0[2] : 0.0));
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      it->first, c1[0], c1[1], (pdb.mesh->getDim() == 3 ? c1[2] : 0.0));

	  foundError = 1;
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Wrond periodic coordinates found on at least on rank!\n");
317
318
319
  }


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

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

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

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

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

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

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

352
353
354
355
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
    for (it_type it = pdb.sendDofs.begin(); it != pdb.sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second[feSpace].begin(); 
	   dofIt != it->second[feSpace].end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
	sendCoords[it->first].push_back(coords[**dofIt]);

358
359
360
    for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
	recvCoords[it->first].push_back(coords[**dofIt]);

363
364
365
    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
366
367
368
369
370
371
372
373
374
375
376
377
378
    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;

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

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

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

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


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

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

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

    int dimOfWorld = Global::getGeo(WORLD);

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

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

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

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


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

467
468
469
470
471
    // 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);
472

473
    typedef map<WorldVector<double>, int> CoordsIndexMap;
474
475
476
477
    CoordsIndexMap coordsToIndex;

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

    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
482
483
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
    for (it_type it = pdb.sendDofs.begin(); it != pdb.sendDofs.end(); ++it)
484
      stdMpi.send(it->first, coordsToIndex);
485
    for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it)
486
487
      stdMpi.recv(it->first);
   
488
    stdMpi.startCommunication();
489
490

    int foundError = 0;
491
    for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it) {
492
493
494
495
496
497
      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) {
498
	  stringstream oss;
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
	  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");
  }

518

519
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
520
  {
521
    FUNCNAME("ParallelDebug::testAllElements()");
522

523
    std::set<int> macroElements;
524
525
    int minElementIndex = numeric_limits<int>::max();
    int maxElementIndex = numeric_limits<int>::min();   
526

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

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

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

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


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

566
567
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

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

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


599
600
601
602
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

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


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

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

651
652
653
654
	typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
	for (it_type rankit = pdb.sendDofs.begin(); rankit != pdb.sendDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second[feSpace].begin();
	       dofit != rankit->second[feSpace].end(); ++dofit)
Thomas Witkowski's avatar
Thomas Witkowski committed
655
	    if (**dofit == it->first)
656
	      cout << "SEND DOF TO " << rankit->first << endl;	  
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
	}

659
660
661
	for (it_type rankit = pdb.recvDofs.begin(); rankit != pdb.recvDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second[feSpace].begin();
	       dofit != rankit->second[feSpace].end(); ++dofit)
Thomas Witkowski's avatar
Thomas Witkowski committed
662
	    if (**dofit == it->first)
663
	      cout << "RECV DOF FROM " << rankit->first << endl;	  
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
	}

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


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

676
677
678
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
679
680
    const FiniteElemSpace* feSpace = pdb.feSpaces[0];

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

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

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

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

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

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

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

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

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

740

741
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
742
  {
743
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
744

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

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

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
767
768
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
769
770
771
772
773
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }    
774
775
  }

776

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

782
783
    const FiniteElemSpace *feSpace = pdb.feSpaces[0];

784
    stringstream filename;
785
786
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

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

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

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

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

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

834
835
836
    file.close();
  }

837
838
839
840
841
842
843
844
845
846
847
848

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
856

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

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

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !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
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
  }


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