diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h index 6f00a4c00999d211f5ec7dfe309c2a38e0b21b27..557623b7f0912d53e83bc65cdc1140497b1cd338 100644 --- a/AMDiS/src/OEMSolver.h +++ b/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; diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index b8b8249d7458e87aa4f3e40d1da14fa265da6d09..fdcfef122f84be08f66689f0275e3de641961a49 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/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); }