ParallelDebug.cc 23.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
#include "ParallelDebug.h"
#include "MeshDistributor.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
15 16 17 18 19
#include "ProblemVec.h"
#include "DOFVector.h"
#include "FixVec.h"
#include "StdMpi.h"
#include "Debug.h"
20
#include "MpiHelper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
21 22 23

namespace AMDiS {

24
  void ParallelDebug::testInteriorBoundary(MeshDistributor &pdb)
Thomas Witkowski's avatar
Thomas Witkowski committed
25
  {
26
    FUNCNAME("ParallelDebug::testInteriorBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
27

28
    typedef MeshDistributor::RankToBoundMap RankToBoundMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
29 30 31 32

    std::vector<int*> sendBuffers, recvBuffers;

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

37 38 39

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

Thomas Witkowski's avatar
Thomas Witkowski committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53
    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);
    }

54 55 56

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

Thomas Witkowski's avatar
Thomas Witkowski committed
57 58 59 60 61 62 63 64 65 66 67
    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);
    }


68 69
    // === To the last, do the same of periodic boundaries. ===

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

75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
      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
95 96 97
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];

98 99 100 101 102

    // === 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
103 104 105 106
    int bufCounter = 0;
    for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
	 rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {

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

111
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
112 113 114 115 116 117 118 119
	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++];
    }
120 121 122 123


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

127 128 129 130
      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;

131 132 133 134
	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);
135 136 137 138
      }

      delete [] recvBuffers[bufCounter++];
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
139 140 141
  }


142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");

    typedef MeshDistributor::PeriodicDofMap PeriodicDofMap;

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

    stdMpi.startCommunication<int>(MPI_INT);

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
      std::map<int, PeriodicDofMap> rankToMaps;
      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) {
	  for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
	       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) {
	for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
	     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");
  }


225
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
226
  {
227
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244

    clock_t first = clock();

    int testCommonDofs = 1;
    GET_PARAMETER(0, "dbg->test common dofs", "%d", &testCommonDofs);
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

    /// Defines a mapping type from rank numbers to sets of coordinates.
    typedef std::map<int, std::vector<WorldVector<double> > > RankToCoords;

    /// Defines a mapping type from rank numbers to sets of DOFs.
    typedef std::map<int, DofContainer> RankToDofContainer;

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

250 251
    // 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
252 253 254 255 256
    // this rank expectes to receive from.
    RankToCoords recvCoords;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin(); 
	   dofIt != it->second.end(); ++dofIt)
	sendCoords[it->first].push_back(coords[**dofIt]);

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

    std::vector<int> sendSize(pdb.mpiSize, 0);
    std::vector<int> recvSize(pdb.mpiSize, 0);
    std::vector<int> recvSizeBuffer(pdb.mpiSize, 0);
    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;

286 287
      request[requestCounter++] = 
	pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289 290 291 292 293
    }   

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

294 295
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
296 297 298 299 300
    }

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


301
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
302 303 304 305 306
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
307 308 309
	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
310 311
      }
    }
312 313
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330

    // === 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.                                ===

    typedef std::vector<WorldVector<double> > CoordsVec;
    StdMpi<CoordsVec> stdMpi(pdb.mpiComm, true);
    stdMpi.send(sendCoords);
    stdMpi.recv(recvCoords);   
    stdMpi.startCommunication<double>(MPI_DOUBLE);

    int dimOfWorld = Global::getGeo(WORLD);

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

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

335 336 337 338
	if (norm(tmp) > 1e-13) {
	  // === Print error message if the coordinates are not the same. ===
	  if (printCoords) {
	    MSG("[DBG] i = %d\n", i);	  
339 340 341
	    std::stringstream oss;
	    oss.precision(5);
	    oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
342
		<< " expect coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
343
	    for (int k = 0; k < dimOfWorld; k++) {
344
	      oss << recvCoords[it->first][i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
345
	      if (k + 1 < dimOfWorld)
346
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
347
	    }
348
	    oss << ")  received coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
349
	    for (int k = 0; k < dimOfWorld; k++) {
350
	      oss << (it->second)[i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	      if (k + 1 < dimOfWorld)
352
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	    }
354 355
	    oss << ")";
	    MSG("%s\n", oss.str().c_str());
356
	    
Thomas Witkowski's avatar
Thomas Witkowski committed
357 358
	    debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
	  }
359 360 361 362
	  ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
	  foundError = 1;
	}	 
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
363
    }
364 365
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
366 367 368 369 370

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


371 372 373 374 375 376 377 378 379 380 381 382
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

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

    typedef std::map<WorldVector<double>, int> CoordsIndexMap;
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it)
383
      coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()];    
384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423

    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      stdMpi.send(it->first, coordsToIndex);
    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it)
      stdMpi.recv(it->first);
   
    stdMpi.startCommunication<double>(MPI_DOUBLE); 

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

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
	  std::stringstream oss;
	  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");
  }

424

425
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
426
  {
427
    FUNCNAME("ParallelDebug::testAllElements()");
428

429 430 431 432
    std::set<int> macroElements;
    int minElementIndex = std::numeric_limits<int>::max();
    int maxElementIndex = std::numeric_limits<int>::min();   

433
    TraverseStack stack;
434
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
435
    while (elInfo) {
436 437 438 439 440 441 442 443 444 445
      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);
446

447 448 449
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
450 451 452 453
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
454 455 456
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
457

458 459 460
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
461 462 463 464 465
      }
    }
  }


466 467 468
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
						    RankToDofContainer &sendDofs,
						    RankToDofContainer &recvDofs)
469
  {
470
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
471 472 473 474 475 476 477 478 479 480

    std::map<int, int> sendNumber;
    for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it)
      sendNumber[it->first] = it->second.size();
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
    for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it)    
      stdMpi.recv(it->first);
    stdMpi.startCommunication<int>(MPI_INT);    
481 482
     
    int foundError = 0;
483
    for (std::map<int, int>::iterator it = stdMpi.getRecvData().begin();
484 485 486 487 488 489 490 491 492 493
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second != static_cast<int>(recvDofs[it->first].size())) {
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
	      recvDofs[it->first].size(), it->first, it->second);
	foundError = 1;
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
494 495 496
  }


497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

    std::map<WorldVector<double>, DegreeOfFreedom> cMap;
    int foundError = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int i = 0; i < 3; i++) {
	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)) {
	    MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f!\n",
		cMap[c], elInfo->getElement()->getDof(i, 0), c[0], c[1]);
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


528
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
529 530 531 532 533 534 535 536 537 538
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << std::endl;
      
      typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
      typedef std::map<int, DofContainer> RankToDofContainer;

      for (DofMapping::iterator it = pdb.mapLocalGlobalDofs.begin();
	   it != pdb.mapLocalGlobalDofs.end(); it++) {
	DegreeOfFreedom localdof = -1;
539 540
	if (pdb.mapLocalDofIndex.count(it->first) > 0)
	  localdof = pdb.mapLocalDofIndex[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
	
	std::cout << "DOF " << it->first << " " 
		  << it->second << " " 
		  << localdof << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, coords);
	coords.print();

	for (RankToDofContainer::iterator rankit = pdb.sendDofs.begin();
	     rankit != pdb.sendDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "SEND DOF TO " << rankit->first << std::endl;	  
	}

	for (RankToDofContainer::iterator rankit = pdb.recvDofs.begin();
	     rankit != pdb.recvDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "RECV DOF FROM " << rankit->first << std::endl;	  
	}

	std::cout << "------" << std::endl;
      }
    }
  }


571
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
572
  {
573
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
574

575 576 577
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
Thomas Witkowski's avatar
Thomas Witkowski committed
578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604
    typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef std::map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;

    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== DOF MAP PERIODIC ====== " << std::endl;

      for (PeriodicDofMap::iterator it = pdb.periodicDof.begin();
	   it != pdb.periodicDof.end(); ++it) {
	std::cout << "DOF MAP " << it->first << ": ";
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
	  std::cout << *dofit << " ";
	std::cout << std::endl;

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

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

	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(localdof, pdb.feSpace, coords);
	coords.print();
      }
    }
605
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
606 607 608
  }

  
609 610 611 612
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636
  {
    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== RANK DOF INFORMATION ====== " << std::endl;

      std::cout << "  RANK OWNED DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }

      std::cout << "  RANK ALL DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }      
    }
  }

637

638
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
639
  {
640
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656

    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);
    }
657 658

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
659 660
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
661 662 663 664 665
      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);
    }    
666 667
  }

668

669 670
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
				     std::string prefix, std::string postfix)
671
  {
672
    FUNCNAME("ParallelDebug::writeCoordsFile()");
673 674 675 676 677 678 679

    std::stringstream filename;
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

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

680 681 682 683 684 685 686 687 688 689 690 691 692 693 694
    typedef std::map<int, std::vector<DegreeOfFreedom> > ElDofMap;
    ElDofMap elDofMap;
    TraverseStack stack;
    const BasisFunction *basisFcts = pdb.feSpace->getBasisFcts();
    std::vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
				 pdb.feSpace->getAdmin(), localIndices);
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

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

695 696 697 698 699
    std::ofstream file;
    file.open(filename.str().c_str());
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
700 701 702
      file << it.getDOFIndex() << " " 
	   << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " "
	   << pdb.getIsRankDof(it.getDOFIndex());
703 704 705 706
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
707 708 709 710 711 712 713 714 715 716 717 718

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

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

719 720 721
    file.close();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
722
}