diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index fcb4dc51cc8edd08d29a286d5b0c5c39094959e7..14c1ab5b64710be236e3130e3bbc726ba4e7dcd8 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -132,23 +132,35 @@ namespace AMDiS { return mesh; } -#if 0 /// Returns an FE space from \ref feSpaces. inline const FiniteElemSpace* getFeSpace(unsigned int i = 0) { FUNCNAME("MeshDistributor::getFeSpace()"); - TEST_EXIT_DBG(i < feSpaces.size())("Should not happen!\n"); + TEST_EXIT_DBG(i < uniqueFeSpaces.size())("Should not happen!\n"); - return feSpaces[i]; + return uniqueFeSpaces[i]; } /// Returns all FE spaces, thus \ref feSpaces. inline vector<const FiniteElemSpace*>& getFeSpaces() { - return feSpaces; + return uniqueFeSpaces; + } + + inline const FiniteElemSpace* getComponentFeSpace(unsigned int i = 0) + { + FUNCNAME("MeshDistributor::getFeSpace()"); + + TEST_EXIT_DBG(i < componentSpaces.size())("Should not happen!\n"); + + return componentSpaces[i]; + } + + inline vector<const FiniteElemSpace*>& getComponentFeSpaces() + { + return componentSpaces; } -#endif /// Returns the DOF mapping object, \ref dofMap. inline ParallelDofMapping& getDofMap() @@ -488,8 +500,8 @@ namespace AMDiS { string name; /// Finite element spaces of the problem. - //vector<const FiniteElemSpace*> feSpaces; vector<const FiniteElemSpace*> componentSpaces; + vector<const FiniteElemSpace*> uniqueFeSpaces; diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc index d84f33780bfcce07d880afac19285c6e468fed24..d4f24f10e2055ca0a5280cc51051acdd346027d3 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.cc +++ b/AMDiS/src/parallel/ParallelDofMapping.cc @@ -191,7 +191,7 @@ namespace AMDiS { TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n"); for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next()) - data->clear(); + it->clear(); nRankDofs = -1; nLocalDofs = -1; @@ -322,6 +322,8 @@ namespace AMDiS { // === Create the matrix indices for all component FE spaces. === + vector<const FiniteElemSpace*> &feSpaces = data->getFeSpaces(); + for (unsigned int i = 0; i < feSpaces.size(); i++) { // Traverse all DOFs of the FE space and create for all rank owned DOFs @@ -462,12 +464,6 @@ namespace AMDiS { { FUNCNAME("ParallelDofMapping::createIndexSet()"); - TEST_EXIT_DBG(firstComponent + nComponents <= feSpaces.size()) - ("Should not happen!\n"); - - TEST_EXIT_DBG(data.count(feSpaces[firstComponent])) - ("No data for FE space at address %p!\n", feSpaces[firstComponent]); - int firstRankDof = -1; ComponentDofMap &compMap = (*data)[firstComponent]; DofMap &dofMap = compMap.getMap(); diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h index 9a8e08d58fb303ebbcc73a5070e2eed3a08e9567..3a86c89aaca8b305dfc06e1d9bd65041e39d2a20 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.h +++ b/AMDiS/src/parallel/ParallelDofMapping.h @@ -295,6 +295,8 @@ namespace AMDiS { public: virtual ComponentDofMap& operator*() = 0; + virtual ComponentDofMap* operator->() = 0; + virtual bool end() = 0; virtual void next() = 0; @@ -319,6 +321,11 @@ namespace AMDiS { vector<const FiniteElemSpace*> &f1, bool isNonLocal) = 0; + vector<const FiniteElemSpace*>& getFeSpaces() + { + return feSpaces; + } + protected: /// The FE spaces for all components. vector<const FiniteElemSpace*> feSpaces; @@ -398,6 +405,10 @@ namespace AMDiS { { } + ComponentDofMap* operator->() + { + } + bool end() { } @@ -416,6 +427,11 @@ namespace AMDiS { ComponentDofMap& operator*() { } + + ComponentDofMap* operator->() + { + } + bool end() { @@ -478,6 +494,11 @@ namespace AMDiS { { } + ComponentDofMap* operator->() + { + } + + bool end() { } @@ -496,6 +517,11 @@ namespace AMDiS { ComponentDofMap& operator*() { } + + ComponentDofMap* operator->() + { + } + bool end() { @@ -618,11 +644,10 @@ namespace AMDiS { } /// Returns the number of solution components the mapping is defined on. - // inline int getNumberOfComponents() const - // { - // ERROR_EXIT("WRITE SOMETHING MEANINGFUL!\n"); - // // return static_cast<int>(feSpaces.size()); - // } + inline int getNumberOfComponents() const + { + return static_cast<int>(data->getFeSpaces().size()); + } /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank. inline int getRankDofs() diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index ccae7d27d6029522a8be81e61b830d4ff23f5092..e8a7ac4631cd513d3167e2525d79a657b8bc4d79 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -200,16 +200,16 @@ namespace AMDiS { interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0); } - for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { - const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); + for (unsigned int i = 0; i < meshDistributor->getComponentFeSpaces().size(); i++) { + const FiniteElemSpace *feSpace = meshDistributor->getComponentFeSpace(i); - createPrimals(feSpace); + createPrimals(i, feSpace); - createDuals(feSpace); + createDuals(i, feSpace); - createInterfaceNodes(feSpace); + createInterfaceNodes(i, feSpace); - createIndexB(feSpace); + createIndexB(i, feSpace); } primalDofMap.update(); @@ -221,10 +221,10 @@ namespace AMDiS { if (stokesMode) interfaceDofMap.update(); - for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { - const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); - createLagrange(feSpace); - createAugmentedLagrange(feSpace); + for (unsigned int i = 0; i < meshDistributor->getComponentFeSpaces().size(); i++) { + const FiniteElemSpace *feSpace = meshDistributor->getComponentFeSpace(i); + createLagrange(i, feSpace); + createAugmentedLagrange(i, feSpace); } lagrangeMap.update(); @@ -295,7 +295,8 @@ namespace AMDiS { } - void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createPrimals(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createPrimals()"); @@ -338,13 +339,14 @@ namespace AMDiS { for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) if (meshDistributor->getDofMap()[feSpace].isRankDof(*it)) { - primalDofMap[feSpace].insertRankDof(*it); + primalDofMap[component].insertRankDof(*it); } else - primalDofMap[feSpace].insertNonRankDof(*it); + primalDofMap[component].insertNonRankDof(*it); } - void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createDuals(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createDuals()"); @@ -361,20 +363,21 @@ namespace AMDiS { if (dirichletRows[feSpace].count(**it)) continue; - if (isPrimal(feSpace, **it)) + if (isPrimal(component, **it)) continue; if (meshLevel == 0) { - dualDofMap[feSpace].insertRankDof(**it); + dualDofMap[component].insertRankDof(**it); } else { if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it)) - dualDofMap[feSpace].insertRankDof(**it); + dualDofMap[component].insertRankDof(**it); } } } - void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createInterfaceNodes(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createInterfaceNodes()"); @@ -390,14 +393,15 @@ namespace AMDiS { continue; if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) - interfaceDofMap[feSpace].insertRankDof(**it); + interfaceDofMap[component].insertRankDof(**it); else - interfaceDofMap[feSpace].insertNonRankDof(**it); + interfaceDofMap[component].insertNonRankDof(**it); } } - void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createLagrange(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createLagrange()"); @@ -442,12 +446,12 @@ namespace AMDiS { meshLevel, feSpace); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) - if (!isPrimal(feSpace, it.getDofIndex())) + if (!isPrimal(component, it.getDofIndex())) if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1) sdRankDofs[it.getRank()].insert(it.getDofIndex()); } - if (dualDofMap[feSpace].nLocalDofs == 0) + if (dualDofMap[component].nLocalDofs == 0) return; @@ -459,7 +463,7 @@ namespace AMDiS { meshLevel, feSpace); !it.end(); it.nextRank()) { for (; !it.endDofIter(); it.nextDof()) { - if (!isPrimal(feSpace, it.getDofIndex())) { + if (!isPrimal(component, it.getDofIndex())) { boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank); if (meshLevel == 0 || @@ -478,7 +482,7 @@ namespace AMDiS { for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) - if (!isPrimal(feSpace, it.getDofIndex())) + if (!isPrimal(component, it.getDofIndex())) if (meshLevel == 0 || (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex()))) stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]); @@ -489,7 +493,7 @@ namespace AMDiS { !it.end(); it.nextRank()) { bool recvFromRank = false; for (; !it.endDofIter(); it.nextDof()) { - if (!isPrimal(feSpace, it.getDofIndex())) { + if (!isPrimal(component, it.getDofIndex())) { if (meshLevel == 0 || (meshLevel > 0 && meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) { @@ -509,14 +513,14 @@ namespace AMDiS { !it.end(); it.nextRank()) { int i = 0; for (; !it.endDofIter(); it.nextDof()) - if (!isPrimal(feSpace, it.getDofIndex())) + if (!isPrimal(component, it.getDofIndex())) if (meshLevel == 0 || (meshLevel > 0 && meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) boundaryDofRanks[feSpace][it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[i++]; else - lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex()); + lagrangeMap[component].insertNonRankDof(it.getDofIndex()); } @@ -527,18 +531,19 @@ namespace AMDiS { DofMap& dualMap = dualDofMap[feSpace].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) { - lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange); + lagrangeMap[component].insertRankDof(it->first, nRankLagrange); int degree = boundaryDofRanks[feSpace][it->first].size(); nRankLagrange += (degree * (degree - 1)) / 2; } else { - lagrangeMap[feSpace].insertNonRankDof(it->first); + lagrangeMap[component].insertNonRankDof(it->first); } } - lagrangeMap[feSpace].nRankDofs = nRankLagrange; + lagrangeMap[component].nRankDofs = nRankLagrange; } - void PetscSolverFeti::createAugmentedLagrange(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createAugmentedLagrange(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createAugmentedLagrange()"); @@ -547,7 +552,8 @@ namespace AMDiS { } - void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace) + void PetscSolverFeti::createIndexB(int component, + const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createIndexB()"); @@ -560,14 +566,14 @@ namespace AMDiS { int nLocalInterior = 0; for (int i = 0; i < admin->getUsedSize(); i++) { if (admin->isDofFree(i) || - isPrimal(feSpace, i) || - isDual(feSpace, i) || - isInterface(feSpace, i) || + isPrimal(component, i) || + isDual(component, i) || + isInterface(component, i) || dirichletRows[feSpace].count(i)) continue; if (meshLevel == 0) { - localDofMap[feSpace].insertRankDof(i, nLocalInterior); + localDofMap[component].insertRankDof(i, nLocalInterior); if (fetiPreconditioner == FETI_DIRICHLET) interiorDofMap[feSpace].insertRankDof(i, nLocalInterior); @@ -575,9 +581,9 @@ namespace AMDiS { nLocalInterior++; } else { if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i)) - localDofMap[feSpace].insertRankDof(i); + localDofMap[component].insertRankDof(i); else - localDofMap[feSpace].insertNonRankDof(i); + localDofMap[component].insertNonRankDof(i); TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE) ("Not yet implemnted!\n"); @@ -589,12 +595,12 @@ namespace AMDiS { for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin(); it != dualDofMap[feSpace].getMap().end(); ++it) { if (meshLevel == 0) { - localDofMap[feSpace].insertRankDof(it->first); + localDofMap[component].insertRankDof(it->first); } else { if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first)) - localDofMap[feSpace].insertRankDof(it->first); + localDofMap[component].insertRankDof(it->first); else - localDofMap[feSpace].insertNonRankDof(it->first); + localDofMap[component].insertNonRankDof(it->first); } } } @@ -785,19 +791,19 @@ namespace AMDiS { int rowCounter = rStartEdges; for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); edgeIt != allEdges.end(); ++edgeIt) { - for (int i = 0; i < feSpaces.size(); i++) { + for (int component = 0; component < feSpaces.size(); component++) { DofContainer edgeDofs; - edgeIt->el->getAllDofs(feSpaces[i], *edgeIt, edgeDofs); + edgeIt->el->getAllDofs(feSpaces[component], *edgeIt, edgeDofs); for (DofContainer::iterator it = edgeDofs.begin(); it != edgeDofs.end(); ++it) { - TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false) + TEST_EXIT_DBG(isPrimal(component, **it) == false) ("Should not be primal!\n"); - if (dirichletRows[feSpaces[i]].count(**it)) + if (dirichletRows[feSpaces[component]].count(**it)) continue; - int col = lagrangeMap.getMatIndex(i, **it); + int col = lagrangeMap.getMatIndex(component, **it); double value = 1.0; MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES); } @@ -1927,7 +1933,7 @@ namespace AMDiS { if (dirichletRows[rowFeSpace].count(*cursor)) continue; - if (isPrimal(rowFeSpace, *cursor)) + if (isPrimal(rowComponent, *cursor)) continue; switch (fetiPreconditioner) { @@ -1944,11 +1950,11 @@ namespace AMDiS { if (dirichletRows[colFeSpace].count(col(*icursor))) continue; - if (isPrimal(colFeSpace, col(*icursor))) + if (isPrimal(colComponent, col(*icursor))) continue; - if (!isDual(rowFeSpace, *cursor)) { - if (!isDual(colFeSpace, col(*icursor))) { + if (!isDual(rowComponent, *cursor)) { + if (!isDual(colComponent, col(*icursor))) { int colIndex = interiorDofMap.getLocalMatIndex(colComponent, col(*icursor)); colsLocal.push_back(colIndex); @@ -1960,7 +1966,7 @@ namespace AMDiS { valuesLocalOther.push_back(value(*icursor)); } } else { - if (!isDual(colFeSpace, col(*icursor))) { + if (!isDual(colComponent, col(*icursor))) { int colIndex = interiorDofMap.getLocalMatIndex(colComponent, col(*icursor)); colsLocalOther.push_back(colIndex); @@ -1974,7 +1980,7 @@ namespace AMDiS { } } // for each nnz in row - if (!isDual(rowFeSpace, *cursor)) { + if (!isDual(rowComponent, *cursor)) { int rowIndex = interiorDofMap.getLocalMatIndex(rowComponent, *cursor); MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(), @@ -1997,7 +2003,7 @@ namespace AMDiS { case FETI_LUMPED: - if (!isDual(rowFeSpace, *cursor)) + if (!isDual(rowComponent, *cursor)) continue; colsLocal.clear(); @@ -2009,7 +2015,7 @@ namespace AMDiS { if (dirichletRows[colFeSpace].count(col(*icursor))) continue; - if (!isDual(colFeSpace, col(*icursor))) + if (!isDual(colComponent, col(*icursor))) continue; int colIndex = diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index 6daa0cae6630b1d140cad976bdd787b90ef67be0..6b5b74daeebdb1abaa6d11f48b31b72e9527fa9a 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -109,23 +109,29 @@ namespace AMDiS { /// Defines which boundary nodes are primal. Creates global index of /// the primal variables. - void createPrimals(const FiniteElemSpace *feSpace); + void createPrimals(int component, + const FiniteElemSpace *feSpace); /// Defines the set of dual variables and creates the global index of /// dual variables. - void createDuals(const FiniteElemSpace *feSpace); + void createDuals(int component, + const FiniteElemSpace *feSpace); /// - void createInterfaceNodes(const FiniteElemSpace *feSpace); + void createInterfaceNodes(int component, + const FiniteElemSpace *feSpace); /// Create Lagrange multiplier variables corresponding to the dual /// variables. - void createLagrange(const FiniteElemSpace *feSpace); + void createLagrange(int component, + const FiniteElemSpace *feSpace); - void createAugmentedLagrange(const FiniteElemSpace *feSpace); + void createAugmentedLagrange(int component, + const FiniteElemSpace *feSpace); /// Creates a global index of the B variables. - void createIndexB(const FiniteElemSpace *feSpace); + void createIndexB(int component, + const FiniteElemSpace *feSpace); /// Creates the Lagrange multiplier constraints and assembles them /// to \ref mat_lagrange. diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index 5df5ff9738d4fc021fefae5271da066985c527cb..d1ed8b8b6bfaf2bb8a1677588565977d0e5fc722 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -150,8 +150,7 @@ namespace AMDiS { for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) { - bool isRowCoarse = - isCoarseSpace(rowComponent, feSpaces[rowComponent], *cursor); + bool isRowCoarse = isCoarseSpace(rowComponent, *cursor); cols.clear(); colsOther.clear(); @@ -162,8 +161,7 @@ namespace AMDiS { for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) { - bool isColCoarse = - isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor)); + bool isColCoarse = isCoarseSpace(colComponent, col(*icursor)); if (isColCoarse == false) if ((*interiorMap)[dofMat->getColFeSpace()].isSet(col(*icursor)) == false) @@ -786,7 +784,7 @@ namespace AMDiS { if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex())) continue; - if (isCoarseSpace(nRowVec, feSpace, dofIt.getDOFIndex())) { + if (isCoarseSpace(nRowVec, dofIt.getDOFIndex())) { TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n"); int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());