Work on debugging functions wihtin ProblemVec

AMDiS/src/OEMSolver.h

@@ -134,7 +134,6 @@ namespace AMDiS {
 	counter++;
       }
 	
-      // std::cout << "Iterative Solver: A = \n" << A.getMatrix();
       int r = solveSystem(A.getMatrix(), xx, bb);
 
       // Copy solution vector to DOFVector
@@ -167,13 +166,12 @@ namespace AMDiS {
       }
 
       // Solver on DOFVector for single system
-      std::cout << "Iterative Solver (System): A = \n" << A.getMatrix();
-      int r= solveSystem(A.getMatrix(), xx, bb);
+      int r = solveSystem(A.getMatrix(), xx, bb);
 
       for (int i = 0, counter = 0; i < ns; i++) {
-	  DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);	  
-	  for (it.reset(); !it.end(); ++it)	
-	      *it = xx[counter++];	
+	DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);	  
+	for (it.reset(); !it.end(); ++it)	
+	  *it = xx[counter++];	
       }
 
       return r;

AMDiS/src/ProblemVec.cc

@@ -765,7 +765,6 @@ namespace AMDiS {
 	if (assembleMatrix)
 	  matrix->finishInsertion();
 
-	// TODO: ExitMatrix should be called after finishInsertion!
  	if (assembleMatrix && matrix->getBoundaryManager())
  	  matrix->getBoundaryManager()->exitMatrix(matrix);	  
 	
@@ -777,7 +776,7 @@ namespace AMDiS {
       assembleBoundaryConditions(rhs->getDOFVector(i),
 				 solution->getDOFVector(i),
 				 componentMeshes[i],
-				 assembleFlag);      
+				 assembleFlag);     
     }
 
     solverMatrix.setMatrix(*systemMatrix);
@@ -1416,7 +1415,7 @@ namespace AMDiS {
   {
     const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
     WorldVector<double> coords;
-    ElementDofIterator elDofIter(componentSpaces[0]);
+    ElementDofIterator elDofIter(componentSpaces[0], true);
 
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(componentMeshes[0], -1, 
@@ -1544,21 +1543,17 @@ namespace AMDiS {
     CoordToDof coordToDof;
     createCoordToDofMap(coordToDof);
 
-    DofToCoord dofToCoord;
-    WorldVector<double> coords;
-    DegreeOfFreedom dof;
     int dimOfWorld = Global::getGeo(WORLD);
-
-    std::vector<std::vector<std::map<std::pair<int, int>, double> > > nnzValues;
-    std::vector<std::map<int, double> > rhsValues;
-    std::vector<std::map<int, double> > solValues;
-
-    nnzValues.resize(nComponents);
-    rhsValues.resize(nComponents);
-    solValues.resize(nComponents);
+    std::vector<std::vector<std::map<std::pair<int, int>, double> > > nnzValues(nComponents);
+    std::vector<std::map<int, double> > rhsValues(nComponents);
+    std::vector<std::map<int, double> > solValues(nComponents);
     for (int i = 0; i < nComponents; i++)
       nnzValues[i].resize(nComponents);
 
+    // Stores to each dof index of this problem a map from rank indices (of each rank
+    // that also has this dof) to the corresponding local dof index.
+    std::map<DegreeOfFreedom, std::vector<DegreeOfFreedom> > dofMapHereToFiles;
+
     for (std::vector<std::string>::iterator fileIt = filenames.begin();
 	 fileIt != filenames.end(); ++fileIt) {
 
@@ -1566,16 +1561,38 @@ namespace AMDiS {
       std::ifstream in(fileIt->c_str());
       int nReadDofs;
       in >> nReadDofs;
-      
+
+      // Is used to map the dof indices in the files to the global coordinates.
+      DofToCoord dofToCoord;
+     
       // Read all DOF indices and their world coordinates.
-      dofToCoord.clear();
       for (int i = 0; i < nReadDofs; i++) {
+	DegreeOfFreedom dof;
+	WorldVector<double> coords;
+
 	in >> dof;
 	for (int j = 0; j < dimOfWorld; j++)
 	  in >> coords[j];
 	dofToCoord[dof] = coords;
       }
 
+      std::map<DegreeOfFreedom, DegreeOfFreedom> dofMapFileToHere;
+
+      for (DofToCoord::iterator it = dofToCoord.begin(); it != dofToCoord.end(); ++it) {
+	DegreeOfFreedom dofIndexInFile = it->first;
+	WorldVector<double> &dofCoords = it->second;
+
+	if (coordToDof.find(dofCoords) == coordToDof.end()) {
+	  std::cout << "Cannot find dof index for coords: " << std::endl;
+	  dofCoords.print();
+	  exit(0);
+	}	
+
+	DegreeOfFreedom dofIndexHere = coordToDof[dofCoords];
+	dofMapHereToFiles[dofIndexHere].push_back(dofIndexInFile);
+	dofMapFileToHere[dofIndexInFile] = dofIndexHere;
+      }
+
       // Now we traverse all matrices and check them.
       for (int i = 0; i < nComponents; i++) {
 	for (int j = 0; j < nComponents; j++) {
@@ -1597,25 +1614,21 @@ namespace AMDiS {
 	    in >> col;
 	    in >> value;
 	    
-	    WorldVector<double> rowCoords = dofToCoord[row];
-	    WorldVector<double> colCoords = dofToCoord[col];
-    
-	    if (coordToDof.find(rowCoords) == coordToDof.end()) {
-	      std::cout << "CANNOT ROW FIND" << std::endl;
-	      rowCoords.print();
+	    if (dofMapFileToHere.count(row) == 0) {
+	      std::cout << "Cannot find row index for: " << row << std::endl;
 	      exit(0);
 	    }
 	    
-	    if (coordToDof.find(colCoords) == coordToDof.end()) {
-	      std::cout << "CANNOT COL FIND" << std::endl;
-	      rowCoords.print();
+	    if (dofMapFileToHere.count(col) == 0) {
+	      std::cout << "Cannot find col index for: " << col << std::endl;
 	      exit(0);
 	    }
 
-	    DegreeOfFreedom rowHere = coordToDof[rowCoords];
-	    DegreeOfFreedom colHere = coordToDof[colCoords];
-
+	    // Get dof indices for row and col of this problem matrix.
+	    DegreeOfFreedom rowHere = dofMapFileToHere[row];
+	    DegreeOfFreedom colHere = dofMapFileToHere[col];
 	    std::pair<int, int> rowcol = make_pair(rowHere, colHere);
+
 	    if (nnzValues[i][j].count(rowcol) == 0)
 	      nnzValues[i][j][rowcol] = value;
 	    else
@@ -1637,7 +1650,7 @@ namespace AMDiS {
 	  WorldVector<double> rowCoords = dofToCoord[row];
 
 	  if (coordToDof.find(rowCoords) == coordToDof.end()) {
-	    std::cout << "CANNOT ROW FIND" << std::endl;
+	    std::cout << "Cannot find row index for coords: " << std::endl;
 	    rowCoords.print();
 	    exit(0);
 	  }
@@ -1664,7 +1677,7 @@ namespace AMDiS {
 	  WorldVector<double> rowCoords = dofToCoord[row];
 
 	  if (coordToDof.find(rowCoords) == coordToDof.end()) {
-	    std::cout << "CANNOT ROW FIND" << std::endl;
+	    std::cout << "Cannot find row index for coords: " << std::endl;
 	    rowCoords.print();
 	    exit(0);
 	  }
@@ -1688,6 +1701,8 @@ namespace AMDiS {
       std::cout << "File read!" << std::endl;
     }
 
+    std::cout << "Now checking ..." << std::endl;
+
     std::cout.precision(15);
 
     for (int i = 0; i < nComponents; i++) {
@@ -1710,10 +1725,13 @@ namespace AMDiS {
 	  // And so we check the difference between the value read from file and the
 	  // corresponding value in the matrix.
 	  if (diff > 1e-8) {
-	    std::cout << "WRONG value in matrix[" << i << "][" << j << "]!" << std::endl
-		      << "  DOF in other matrix[" << row << "][" << col << "] = " << value << std::endl
-		      << "  DOF in this matrix[" << row << "][" << col << "] = " << valueHere << std::endl;
-	    
+	    std::cout << "Wrong value in component matrix " << i 
+		      << "/" << j << ": " << std::endl
+		      << "  DOF in this matrix[" << row << "][" << col << "] = " 
+		      << valueHere << std::endl	    	    
+		      << "  DOF in other matrix[" << dofMapHereToFiles[row][0] << "][" 
+		      << dofMapHereToFiles[col][0] << "] = " 
+		      << value << std::endl;
 	    
 	    exit(0);
 	  }