diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index ee43d7d25c55b8376911d78e44066c5d6f7b7b71..6c4d4927813f94d7b4fafd7394d63e1e6a07c88c 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -259,31 +259,28 @@ namespace AMDiS { } MSG("FETI-DP data created on mesh level %d\n", meshLevel); - for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { + for (unsigned int i = 0; i < meshDistributor->getComponentFeSpaces().size(); i++) { const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); - MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace); + + MSG("FETI-DP data for %d-ith component (FE space %p):\n", i, feSpace); if (feSpace == pressureFeSpace) { MSG(" nRankInterface = %d nOverallInterface = %d\n", - interfaceDofMap[feSpace].nRankDofs, - interfaceDofMap[feSpace].nOverallDofs); + interfaceDofMap[i].nRankDofs, + interfaceDofMap[i].nOverallDofs); } else { MSG(" nRankPrimals = %d nLocalPrimals = %d nOverallPrimals = %d\n", - primalDofMap[feSpace].nRankDofs, - primalDofMap[feSpace].nLocalDofs, - primalDofMap[feSpace].nOverallDofs); + primalDofMap[i].nRankDofs, + primalDofMap[i].nLocalDofs, + primalDofMap[i].nOverallDofs); MSG(" nRankDuals = %d nOverallDuals = %d\n", - dualDofMap[feSpace].nRankDofs, - dualDofMap[feSpace].nOverallDofs); + dualDofMap[i].nRankDofs, + dualDofMap[i].nOverallDofs); MSG(" nRankLagrange = %d nOverallLagrange = %d\n", - lagrangeMap[feSpace].nRankDofs, - lagrangeMap[feSpace].nOverallDofs); - -// TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == -// static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs())) -// ("Should not happen!\n"); + lagrangeMap[i].nRankDofs, + lagrangeMap[i].nOverallDofs); } } @@ -534,7 +531,7 @@ namespace AMDiS { // === appropriate number of Lagrange constraints. === int nRankLagrange = 0; - DofMap& dualMap = dualDofMap[feSpace].getMap(); + DofMap& dualMap = dualDofMap[component].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) { lagrangeMap[component].insertRankDof(it->first, nRankLagrange); @@ -598,8 +595,8 @@ namespace AMDiS { // === And finally, add the global indicies of all dual nodes. === - for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin(); - it != dualDofMap[feSpace].getMap().end(); ++it) { + for (DofMap::iterator it = dualDofMap[component].getMap().begin(); + it != dualDofMap[component].getMap().end(); ++it) { if (meshLevel == 0) { localDofMap[component].insertRankDof(it->first); } else { @@ -638,19 +635,19 @@ namespace AMDiS { // === m == r, than the rank sets -1.0 for the corresponding === // === constraint. === - for (unsigned int k = 0; k < feSpaces.size(); k++) { - DofMap &dualMap = dualDofMap[feSpaces[k]].getMap(); + for (unsigned int component = 0; component < feSpaces.size(); component++) { + DofMap &dualMap = dualDofMap[component].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { - TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first)) + TEST_EXIT_DBG(boundaryDofRanks[feSpaces[component]].count(it->first)) ("Should not happen!\n"); // Global index of the first Lagrange constriant for this node. - int index = lagrangeMap.getMatIndex(k, it->first); + int index = lagrangeMap.getMatIndex(component, it->first); // Copy set of all ranks that contain this dual node. - vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), - boundaryDofRanks[feSpaces[k]][it->first].end()); + vector<int> W(boundaryDofRanks[feSpaces[component]][it->first].begin(), + boundaryDofRanks[feSpaces[component]][it->first].end()); // Number of ranks that contain this dual node. int degree = W.size(); @@ -662,7 +659,7 @@ namespace AMDiS { if (W[i] == mpiRank || W[j] == mpiRank) { MatSetValue(mat_lagrange, index + counter, - localDofMap.getMatIndex(k, it->first) + rStartInterior, + localDofMap.getMatIndex(component, it->first) + rStartInterior, (W[i] == mpiRank ? 1.0 : -1.0), INSERT_VALUES); } @@ -673,7 +670,7 @@ namespace AMDiS { // === Create scaling factors for scaling the lagrange matrix, which === // === is required for FETI-DP preconditioners. === - if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) { + if (meshDistributor->getDofMap()[feSpaces[component]].isRankDof(it->first)) { int nConstraints = (degree * (degree - 1)) / 2; for (int i = 0; i < nConstraints; i++) { VecSetValue(vec_scale_lagrange, @@ -1228,12 +1225,12 @@ namespace AMDiS { TEST_EXIT_DBG(meshLevel == 0) ("Should not happen, check usage of localDofMap!\n"); - for (unsigned int i = 0; i < feSpaces.size(); i++) { - DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); + for (unsigned int component = 0; component < feSpaces.size(); component++) { + DofMap &dualMap = dualDofMap[component].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { DegreeOfFreedom d = it->first; - int matIndexLocal = localDofMap.getLocalMatIndex(i, d); - int matIndexDual = dualDofMap.getLocalMatIndex(i, d); + int matIndexLocal = localDofMap.getLocalMatIndex(component, d); + int matIndexDual = dualDofMap.getLocalMatIndex(component, d); fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual; } } @@ -1271,15 +1268,15 @@ namespace AMDiS { MatGetVecs(mat_duals_duals, PETSC_NULL, &(lumpedData->tmp_vec_duals1)); - for (unsigned int i = 0; i < feSpaces.size(); i++) { - if (stokesMode && i == pressureComponent) + for (unsigned int component = 0; component < feSpaces.size(); component++) { + if (stokesMode && component == pressureComponent) continue; - DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); + DofMap &dualMap = dualDofMap[component].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { DegreeOfFreedom d = it->first; - int matIndexLocal = localDofMap.getLocalMatIndex(i, d); - int matIndexDual = dualDofMap.getLocalMatIndex(i, d); + int matIndexLocal = localDofMap.getLocalMatIndex(component, d); + int matIndexDual = dualDofMap.getLocalMatIndex(component, d); lumpedData->localToDualMap[matIndexLocal] = matIndexDual; } } @@ -1588,7 +1585,7 @@ namespace AMDiS { std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows(); for (std::set<DegreeOfFreedom>::iterator it = dirichletRows.begin(); it != dirichletRows.end(); ++it) { - if (localDofMap[feSpace].isSet(*it)) { + if (localDofMap[i].isSet(*it)) { MSG("Dirichlet dof %d in component %d with interior mat index %d\n", *it, i, localDofMap.getMatIndex(i, *it)); } @@ -1640,7 +1637,7 @@ namespace AMDiS { // === Create scatter to get solutions of all primal nodes that are === // === contained in rank's domain. === - vector<const FiniteElemSpace*> feSpaces = vec.getComponentFeSpaces(); + int nComponents = vec.getSize(); vector<PetscInt> globalIsIndex, localIsIndex; globalIsIndex.reserve(primalDofMap.getLocalDofs()); @@ -1648,10 +1645,10 @@ namespace AMDiS { { int cnt = 0; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - DofMap& dofMap = primalDofMap[feSpaces[i]].getMap(); + for (unsigned int component = 0; component < nComponents; component++) { + DofMap& dofMap = primalDofMap[component].getMap(); for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) { - globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first)); + globalIsIndex.push_back(primalDofMap.getMatIndex(component, it->first)); localIsIndex.push_back(cnt++); } } @@ -1692,26 +1689,27 @@ namespace AMDiS { // === And copy from PETSc local vectors to the DOF vectors. === + vector<const FiniteElemSpace*> feSpaces = vec.getComponentFeSpaces(); int cnt = 0; - for (int i = 0; i < vec.getSize(); i++) { - DOFVector<double>& dofVec = *(vec.getDOFVector(i)); + for (int component = 0; component < nComponents; component++) { + DOFVector<double>& dofVec = *(vec.getDOFVector(component)); - for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap().begin(); - it != localDofMap[feSpaces[i]].getMap().end(); ++it) { + for (DofMap::iterator it = localDofMap[component].getMap().begin(); + it != localDofMap[component].getMap().end(); ++it) { if (meshLevel == 0) { - int petscIndex = localDofMap.getLocalMatIndex(i, it->first); + int petscIndex = localDofMap.getLocalMatIndex(component, it->first); dofVec[it->first] = localSolB[petscIndex]; } else { - if (meshDistributor->getDofMapSd()[feSpaces[i]].isRankDof(it->first)) { - int petscIndex = localDofMap.getLocalMatIndex(i, it->first); + if (meshDistributor->getDofMapSd()[feSpaces[component]].isRankDof(it->first)) { + int petscIndex = localDofMap.getLocalMatIndex(component, it->first); TEST_EXIT(petscIndex < bsize)("Should not happen!\n"); dofVec[it->first] = localSolB[petscIndex]; } } } - for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap().begin(); - it != primalDofMap[feSpaces[i]].getMap().end(); ++it) + for (DofMap::iterator it = primalDofMap[component].getMap().begin(); + it != primalDofMap[component].getMap().end(); ++it) dofVec[it->first] = localSolPrimal[cnt++]; }