Fixed several bugs in parallel code. Added more parallel check and debug functions.

AMDiS/src/Debug.cc

@@ -178,6 +178,21 @@ namespace AMDiS {
     }
 
 
+    Element* getLevel0ParentElement(Mesh *mesh, int elIndex)
+    {    
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+      while (elInfo) {
+	if (elInfo->getElement()->getIndex() == elIndex)
+	  return elInfo->getMacroElement()->getElement();      
+	
+	elInfo = stack.traverseNext(elInfo);
+      }
+      
+      return NULL;
+    }
+
+
     Element* getParentElement(Mesh *mesh, Element *el)
     {
       TraverseStack stack;
@@ -243,6 +258,22 @@ namespace AMDiS {
     }
 
     
+    void printElementCoords(const FiniteElemSpace *feSpace, Element *el)
+    {
+      FUNCNAME("debug:printElementCoords()");
+      
+      Mesh *mesh = feSpace->getMesh();
+
+      for (int i = 0; i <= mesh->getDim(); i++) {
+	DegreeOfFreedom dof = el->getDOF(i, 0);
+	WorldVector<double> coords;
+	mesh->getDofIndexCoords(dof, feSpace, coords);
+	MSG("%d-th DOF of element %d: %f %f %f\n", 
+	    i, el->getIndex(), coords[0], coords[1], coords[2]);
+      }
+    }
+
+
     void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
     {
       FUNCNAME("debug::printInfoByDof()");
@@ -557,6 +588,48 @@ namespace AMDiS {
     }
 
 
+    void importDofVectorByCoords(DOFVector<double>* vec, std::string filename)
+    { 
+      DOFVector<WorldVector<double> > coords(vec->getFeSpace(), "dofCoords");
+      vec->getFeSpace()->getMesh()->getDofIndexCoords(vec->getFeSpace(), coords);
+      int dim = vec->getFeSpace()->getMesh()->getDim();
+
+      std::ifstream file;
+      file.open(filename.c_str());
+      
+      int nEntries;
+      file >> nEntries;
+
+      for (int i = 0; i < nEntries; i++) {
+	WorldVector<double> dofCoords;
+	for (int j = 0; j < dim; j++)
+	  file >> dofCoords[j];
+
+	double value;
+	file >> value;
+
+      
+	DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
+	for (it.reset(); !it.end(); ++it) {
+	  bool found = true;
+	  for (int j = 0; j < dim; j++) {
+	    if (fabs((*it)[j] - dofCoords[j]) > 1e-8) {
+	      found = false;
+	      break;
+	    }
+	  }
+
+	  if (found) {
+	    (*vec)[it.getDOFIndex()] = value;
+	    break;
+	  }
+	}	
+      }
+
+      file.close();
+    }
+
+
     void exportDofVectorByCoords(const DOFVector<double>* vec, 
 				 std::string filename)
     {

AMDiS/src/Debug.h

@@ -86,6 +86,8 @@ namespace AMDiS {
     
     Element* getLevel0ParentElement(Mesh *mesh, Element *el);
 
+    Element* getLevel0ParentElement(Mesh *mesh, int elIndex);
+
     Element* getParentElement(Mesh *mesh, Element *el);
 
     Element* getParentElement(Mesh *mesh, int elIndex);
@@ -93,6 +95,8 @@ namespace AMDiS {
     Element* getElement(Mesh *mesh, int elIndex);
 
     void printElementInfo(Element *el);
+
+    void printElementCoords(const FiniteElemSpace *feSpace, Element *el);
     
     void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
 
@@ -141,6 +145,8 @@ namespace AMDiS {
 
     int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex);
 
+    void importDofVectorByCoords(DOFVector<double>* vec, std::string filename);
+
     void exportDofVectorByCoords(const DOFVector<double>* vec, 
 				 std::string filename);
 

AMDiS/src/MatrixVector.h

@@ -171,7 +171,7 @@ namespace AMDiS {
     /// Access to the i-th vector element for const vectors.
     inline const T& operator[] (int i) const 
     {
-      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
+      TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i);
       return valArray[i];
     }
 

AMDiS/src/Tetrahedron.cc

@@ -214,6 +214,10 @@ namespace AMDiS {
 					BoundaryObject bound,
 					DofContainer& dofs) const
   {
+    FUNCNAME("Tetrahedron::getVertexDofsAtFace()");
+
+    //    MSG("go at el %d with face %d\n", this->getIndex(), bound.ithObj);
+
     if (!child[0])
       return;
 
@@ -246,7 +250,7 @@ namespace AMDiS {
 	for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
 	     it != bound.excludedSubstructures.end(); ++it)
 	  if (it->first == EDGE && it->second == 0)
-	    addDof = false;
+	    addDof = false;	
 	
 	if (bound.reverseMode) {
 	  child[1]->getVertexDofs(feSpace, nextBound1, dofs);
@@ -327,8 +331,8 @@ namespace AMDiS {
     case FACE:      
       for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
    	   it != bound.excludedSubstructures.end(); ++it) {
-	if (it->first == EDGE)
-	  it->second = edgeOfChild[bound.elType][ithChild][it->second];
+	if (it->first == EDGE && it->second != -1)
+	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
       }
       
       bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];

AMDiS/src/ZeroOrderAssembler.cc

@@ -68,9 +68,9 @@ namespace AMDiS {
       newAssembler = new StandardZOA(op, assembler, quad);
     } else {
       if (pwConst)
-	newAssembler = new PrecalcZOA(op, assembler, quad);
+     	newAssembler = new PrecalcZOA(op, assembler, quad);
       else
-	newAssembler = new FastQuadZOA(op, assembler, quad);      
+     	newAssembler = new FastQuadZOA(op, assembler, quad);      
     }
 
     subAssemblers->push_back(newAssembler);

AMDiS/src/ZeroOrderTerm.cc

@@ -36,7 +36,7 @@ namespace AMDiS {
   { 
     if (f) {
       for (int iq = 0; iq < nPoints; iq++)
-	C[iq] += (*f)(vecAtQPs[iq]);      
+	C[iq] += (*f)(vecAtQPs[iq]);            
     } else {
       for (int iq = 0; iq < nPoints; iq++)
 	C[iq] += vecAtQPs[iq];      

AMDiS/src/parallel/MeshDistributor.cc

@@ -147,6 +147,7 @@ namespace AMDiS {
     debug::testSortedDofs(mesh, elMap);
     ParallelDebug::testInteriorBoundary(*this);
     ParallelDebug::testCommonDofs(*this, true);
+    ParallelDebug::testGlobalIndexByCoords(*this);
     MSG("Debug mode tests finished!\n");
 
     debug::writeMesh(feSpace, -1, "macro_mesh");   
@@ -595,19 +596,6 @@ namespace AMDiS {
 	if (elCode.getCode() != recvCodes[i].getCode()) {
 	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
 
-#if 0
-	  MSG("CHECK EL %d %d %d %d WITH El %d %d %d %d in RANK %d\n",
-	      boundIt->rankObj.elIndex,
-	      boundIt->rankObj.subObj,
-	      boundIt->rankObj.ithObj, 
-	      boundIt->rankObj.elType,	      
-	      boundIt->neighObj.elIndex,
-	      boundIt->neighObj.subObj,
-	      boundIt->neighObj.ithObj, 
-	      boundIt->neighObj.elType,	      
-	      it->first);
-#endif
-
 	  bool b = startFitElementToMeshCode(recvCodes[i], 
 					     boundIt->rankObj.el,
 					     boundIt->rankObj.subObj,
@@ -699,7 +687,6 @@ namespace AMDiS {
       refineManager->setMesh(el->getMesh());
       refineManager->setStack(&stack);
       refineManager->refineFunction(elInfo);
-      MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex());
       meshChanged = true;
     }
 
@@ -778,7 +765,6 @@ namespace AMDiS {
       refineManager->setMesh(el->getMesh());
       refineManager->setStack(&stack);
       refineManager->refineFunction(elInfo);
-      MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex());
       meshChanged = true;
     }
 
@@ -1221,9 +1207,14 @@ namespace AMDiS {
 	  bound.rankObj.subObj = geoIndex;
 	  bound.rankObj.ithObj = rankBoundEl.ithObject;
 
- 	  if (geoIndex == FACE)
- 	    for (int edgeNo = 0; edgeNo < 6; edgeNo++)
- 	      bound.rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
+ 	  if (geoIndex == FACE) {
+ 	    for (int edgeNo = 0; edgeNo < 3; edgeNo++) {
+	      int edgeOfFace = 
+		bound.rankObj.el->getEdgeOfFace(bound.rankObj.ithObj, edgeNo);
+
+	      bound.rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeOfFace));
+	    }
+	  }
 	  
 	  
 	  if (owner == mpiRank) {
@@ -1601,7 +1592,7 @@ namespace AMDiS {
       it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
 
       for (unsigned int i = 0; i < dofs.size(); i++)
-	sendDofs[it.getRank()].push_back(dofs[i]);          
+	sendDofs[it.getRank()].push_back(dofs[i]);
     }
 
 
@@ -1638,7 +1629,7 @@ namespace AMDiS {
     // the dofs in ranks partition that are owned by other rank are set to false.
     isRankDof.clear();
     for (int i = 0; i < nRankAllDofs; i++)
-      isRankDof[i] = true;
+      isRankDof[i] = true;    
 
     // Stores for all rank owned dofs a new global index.
     DofIndexMap rankDofsNewGlobalIndex;
@@ -1697,8 +1688,14 @@ namespace AMDiS {
     lastMeshChangeIndex = mesh->getChangeIndex();
     
 #if (DEBUG != 0)
+    std::stringstream oss;
+    oss << "elementIndex-" << mpiRank << ".vtu";
+    debug::writeElementIndexMesh(mesh, oss.str());
+
     debug::testSortedDofs(mesh, elMap);
     ParallelDebug::testCommonDofs(*this, true);
+    ParallelDebug::testGlobalIndexByCoords(*this);
+    ParallelDebug::writeDebugFile(*this, "mpi-dbg", "dat");
     
 #if 0
     MSG("------------- Debug information -------------\n");
@@ -1721,16 +1718,6 @@ namespace AMDiS {
 	MSG("%d recv DOF: %d\n", i, *(it->second[i]));
     }
     */
-
-#if 0
-    std::set<const DegreeOfFreedom*> testDofs;
-    debug::getAllDofs(feSpace, testDofs);
-    for (std::set<const DegreeOfFreedom*>::iterator it = testDofs.begin();
-	 it != testDofs.end(); ++it) {
-      MSG("DOF %d:   mapLocalGlobalDofs = %d   isRankDof = %d\n", **it, 
-	  mapLocalGlobalDofs[**it], isRankDof[**it]);
-    }
-#endif
 #endif
 #endif
   }

AMDiS/src/parallel/ParallelDebug.cc

@@ -265,6 +265,59 @@ namespace AMDiS {
   }
 
 
+  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)
+      coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()];
+
+    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");
+  }
+
   void ParallelDebug::testAllElements(MeshDistributor &pdb)
   {
     FUNCNAME("ParallelDebug::testAllElements()");
@@ -467,8 +520,8 @@ namespace AMDiS {
   }
 
 
-  void ParallelDebug::writeCoordsFile(MeshDistributor &pdb, 
-				      std::string prefix, std::string postfix)
+  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
+				     std::string prefix, std::string postfix)
   {
     FUNCNAME("ParallelDebug::writeCoordsFile()");
 
@@ -478,16 +531,45 @@ namespace AMDiS {
     DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp");
     pdb.mesh->getDofIndexCoords(pdb.feSpace, coords);
 
+    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. ===
+
     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) {
-      file << pdb.mapLocalGlobalDofs[it.getDOFIndex()];
+      file << it.getDOFIndex() << " " 
+	   << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " "
+	   << pdb.getIsRankDof(it.getDOFIndex());
       for (int i = 0; i < pdb.mesh->getDim(); i++)
 	file << " " << (*it)[i];
       file << "\n";
     }
+
+    // === 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";
+    }
+
     file.close();
   }
 

AMDiS/src/parallel/ParallelDebug.h

@@ -42,7 +42,7 @@ namespace AMDiS {
 
     /** \brief
      * This function is used for debugging only. It traverses all interior boundaries
-     * and compares the dof indices on them with the dof indices of the boundarys
+     * and compares the DOF indices on them with the dof indices of the boundarys
      * neighbours. The function fails, when dof indices on an interior boundary do
      * not fit together.
      *
@@ -52,6 +52,16 @@ namespace AMDiS {
      */
     static void testCommonDofs(MeshDistributor &pdb, bool printCoords = false);
 
+    /** \brief
+     * This function is used for debugging only. It checks if on all ranks DOFs with
+     * the same coordinates have the same global index. For this, the function genartes
+     * on all ranks a list of all DOFs with their coordinates and global indices and 
+     * sends the list to all neighbouring ranks.
+     *
+     * \param[in]  pdb           Parallel problem definition used for debugging.
+     */
+    static void testGlobalIndexByCoords(MeshDistributor &pdb);
+
     /** \brief
      * Tests if all elements in the macro mesh are memeber of exactly one rank.
      *
@@ -114,8 +124,8 @@ namespace AMDiS {
      */
     static void printBoundaryInfo(MeshDistributor &pdb);
 
-    static void writeCoordsFile(MeshDistributor &pdb, 
-				std::string prefix, std::string postfix);
+    static void writeDebugFile(MeshDistributor &pdb,
+			       std::string prefix, std::string postfix);
   };
   
 } // namespace AMDiS

AMDiS/src/parallel/PetscSolver.cc

@@ -30,7 +30,7 @@ namespace AMDiS {
     clock_t first = clock();
 
     fillPetscMatrix(systemMatrix, rhs);
-    solvePetscMatrix(*solution, adaptInfo);
+    solvePetscMatrix(*solution, adaptInfo);   
 
 #ifdef _OPENMP
     INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
@@ -94,7 +94,7 @@ namespace AMDiS {
 	// Calculate the exact position of the column index in the petsc matrix.
 	int colIndex = globalColDof * dispMult + dispAddCol;
 
-	// Set only non null values.
+	// Set only non zero values.
 	if (value(*icursor) != 0.0 || rowIndex == colIndex) {
 
 	  // If the current row is not periodic, but the current dof index is periodic,
@@ -116,7 +116,7 @@ namespace AMDiS {
  	      values.push_back(value(*icursor) * scalFactor);
 	    }
 
- 	  } else {
+ 	  } else {	    
 	    // The col dof index is not periodic, simple add entry.
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
@@ -409,7 +409,6 @@ namespace AMDiS {
     // === Create PETSc matrix with the computed nnz data structure. ===
 
     MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
-		    //		    0, PETSC_NULL, 0, PETSC_NULL, &petscMatrix);
 		    0, d_nnz, 0, o_nnz, &petscMatrix);
     
 #if (DEBUG != 0)
@@ -481,10 +480,6 @@ namespace AMDiS {
     // Do not delete the solution vector, use it for the initial guess.
     //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
 
-#if (DEBUG != 0)
-    //    ParallelDomainDbg::writeCoordsFile(*meshDistributor, "mpi-coords", "dat");
-#endif
-
 
     // === Run Petsc. ===
 

AMDiS/src/parallel/StdMpi.cc

@@ -12,6 +12,11 @@ namespace AMDiS {
     return data.size() * 2;
   }
 
+  int intSizeOf(std::map<WorldVector<double>, int> &data)
+  {
+    return data.size() * (Global::getGeo(WORLD) + 1);
+  }
+
   int intSizeOf(std::vector<MeshStructure> &data)
   {
     int s = 0;
@@ -91,6 +96,38 @@ namespace AMDiS {
       data.push_back(std::make_pair(buf[i * 2], buf[i * 2 + 1]));
   }
 
+  void makeBuf(std::map<WorldVector<double>, int> &data, double* buf)
+  {
+    int i = 0;
+    for (std::map<WorldVector<double>, int>::iterator it = data.begin();
+	 it != data.end(); ++it) {
+      for (int j = 0; j < Global::getGeo(WORLD); j++)
+	buf[i++] = it->first[j];
+      buf[i++] = static_cast<double>(it->second);      
+    }
+  }
+
+  void makeFromBuf(std::map<WorldVector<double>, int> &data, double* buf, int bufSize)
+  {
+    if (bufSize == 0)
+      return;
+
+    int oneEntrySize = Global::getGeo(WORLD) + 1;
+    int nEntry = bufSize / oneEntrySize;
+
+    TEST_EXIT(bufSize % oneEntrySize == 0)("This should not happen!\n");
+
+    data.clear();
+    int i = 0;
+    WorldVector<double> coords;
+
+    for (int j = 0; j < nEntry; j++) {
+      for (int k = 0; k < Global::getGeo(WORLD); k++)
+	coords[k] = buf[i++];
+      data[coords] = static_cast<int>(buf[i++]);
+    }
+  }
+
   void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf)
   {    
     int pos = 0;

AMDiS/src/parallel/StdMpi.h

@@ -33,6 +33,8 @@ namespace AMDiS {
 
   int intSizeOf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data);
 
+  int intSizeOf(std::map<WorldVector<double>, int> &data);
+
   int intSizeOf(std::vector<MeshStructure> &data);
 
   int intSizeOf(std::vector<WorldVector<double> > &data);
@@ -58,6 +60,10 @@ namespace AMDiS {
   void makeFromBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data, 
 		   int *buf, int bufSize);
 
+  void makeBuf(std::map<WorldVector<double>, int> &data, double* buf);
+
+  void makeFromBuf(std::map<WorldVector<double>, int> &data, double* buf, int bufSize);
+
   void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf);
 
   void makeFromBuf(std::vector<MeshStructure> &data,