diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index ca2b0dfff6f935366d64e49ced293594aa601ab8..96d23346f3351ec0ba230b71bd9734b9cc6cd032 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -1687,6 +1687,7 @@ namespace AMDiS { MSG("| number of levels: %d\n", nLevels); MSG("| number of FE spaces: %d\n", uniqueFeSpaces.size()); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) { MSG("| FE space = %d (pointer adr %p):\n", i, uniqueFeSpaces[i]); MSG("| nRankDofs = %d\n", dofMap[uniqueFeSpaces[i]].nRankDofs); diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc index b567fa797d682cbf658bf2c6859daab87e649df8..c42017ce29c262747026d88a22988ad764e1173e 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.cc +++ b/AMDiS/src/parallel/ParallelDofMapping.cc @@ -42,8 +42,7 @@ namespace AMDiS { : levelData(ld), dofComm(NULL), feSpace(NULL), - needGlobalMapping(false), - isNonLocal(false) + globalMapping(false) { clear(); } @@ -80,11 +79,9 @@ namespace AMDiS { // === If required, compute also the global indices. === - if (needGlobalMapping) { - computeGlobalMapping(); - - if (isNonLocal) - computeNonLocalIndices(); + if (globalMapping) { + computeGlobalMapping(); + computeNonLocalIndices(); } } @@ -170,67 +167,50 @@ namespace AMDiS { } - void ComponentDataEqFeSpace::init(vector<const FiniteElemSpace*> &f0, - vector<const FiniteElemSpace*> &f1, - bool isNonLocal, - MeshLevelData &levelData) + void FeSpaceData::init(vector<const FiniteElemSpace*> &f0, + vector<const FiniteElemSpace*> &f1, + bool globalMapping, + MeshLevelData &levelData) { - feSpaces = f0; - feSpacesUnique = f1; - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) { - addFeSpace(*it, levelData); - componentData[*it].setNeedGlobalMapping(isNonLocal); - componentData[*it].setNonLocal(isNonLocal); + componentSpaces = f0; + feSpaces = f1; + for (vector<const FiniteElemSpace*>::iterator it = feSpaces.begin(); + it != feSpaces.end(); ++it) { + if (componentData.count(*it)) + componentData.find(*it)->second.clear(); + else + componentData.insert(make_pair(*it, ComponentDofMap(&levelData))); + + componentData[*it].setFeSpace(*it); + componentData[*it].setGlobalMapping(globalMapping); } } - void ComponentDataEqFeSpace::addFeSpace(const FiniteElemSpace* feSpace, - MeshLevelData &levelData) + void ComponentData::init(vector<const FiniteElemSpace*> &f0, + vector<const FiniteElemSpace*> &f1, + bool globalMapping, + MeshLevelData &levelData) { - if (componentData.count(feSpace)) - componentData.find(feSpace)->second.clear(); - else - componentData.insert(make_pair(feSpace, ComponentDofMap(&levelData))); - - componentData.find(feSpace)->second.setFeSpace(feSpace); - } - - - void ComponentDataDiffFeSpace::init(vector<const FiniteElemSpace*> &f0, - vector<const FiniteElemSpace*> &f1, - bool isNonLocal, - MeshLevelData &levelData) - { - feSpaces = f0; - feSpacesUnique = f1; - - for (unsigned int component = 0; component < feSpaces.size(); component++) { - addComponent(component, feSpaces[component], levelData); - componentData[component].setNeedGlobalMapping(isNonLocal); - componentData[component].setNonLocal(isNonLocal); + componentSpaces = f0; + feSpaces = f1; + + for (unsigned int component = 0; component < componentSpaces.size(); component++) { + if (componentData.count(component)) + componentData.find(component)->second.clear(); + else + componentData.insert(make_pair(component, ComponentDofMap(&levelData))); + + componentData[component].setFeSpace(componentSpaces[component]); + componentData[component].setGlobalMapping(globalMapping); } } - void ComponentDataDiffFeSpace::addComponent(unsigned int component, - const FiniteElemSpace* feSpace, - MeshLevelData &levelData) - { - if (componentData.count(component)) - componentData.find(component)->second.clear(); - else - componentData.insert(make_pair(component, ComponentDofMap(&levelData))); - - componentData.find(component)->second.setFeSpace(feSpace); - } - - ParallelDofMapping::ParallelDofMapping(DofMappingMode mode) : levelData(NULL), dofComm(NULL), - isNonLocal(true), + globalMapping(true), needMatIndexFromGlobal(false), nRankDofs(1), nLocalDofs(1), @@ -239,10 +219,10 @@ namespace AMDiS { { switch (mode) { case COMPONENT_WISE: - data = new ComponentDataDiffFeSpace(); + data = new ComponentData(); break; case FESPACE_WISE: - data = new ComponentDataEqFeSpace(); + data = new FeSpaceData(); break; } @@ -261,9 +241,10 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::init()"); levelData = &ldata; - isNonLocal = b; + globalMapping = b; + componentSpaces = fe; - data->init(fe, uniqueFe, isNonLocal, ldata); + data->init(fe, uniqueFe, globalMapping, ldata); } @@ -386,29 +367,28 @@ 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++) { + for (unsigned int component = 0; component < componentSpaces.size(); + component++) { // Traverse all DOFs of the FE space and create for all rank owned DOFs // a matrix index. - DofMap& dofMap = (*data)[i].getMap(); + DofMap& dofMap = (*data)[component].getMap(); for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { - if ((*data)[i].isRankDof(it->first)) { + if ((*data)[component].isRankDof(it->first)) { int globalMatIndex = it->second.local + offset; if (globalIndex) - dofToMatIndex.add(i, it->second.global, globalMatIndex); + dofToMatIndex.add(component, it->second.global, globalMatIndex); else - dofToMatIndex.add(i, it->first, globalMatIndex); + dofToMatIndex.add(component, it->first, globalMatIndex); } } // Increase the offset for the next FE space by the number of DOFs owned // by the rank in the current FE space. - offset += (*data)[i].nRankDofs; + offset += (*data)[component].nRankDofs; // If there are no non local DOFs, continue with the next FE space. - if (!isNonLocal) + if (!globalMapping) continue; TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n"); @@ -417,16 +397,17 @@ namespace AMDiS { // === interior boundaries. === StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm); - for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); + for (DofComm::Iterator it(dofComm->getSendDofs(), 0, componentSpaces[component]); !it.end(); it.nextRank()) { vector<DegreeOfFreedom> sendGlobalDofs; for (; !it.endDofIter(); it.nextDof()) if (dofMap.count(it.getDofIndex())) if (globalIndex) - sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global)); + sendGlobalDofs.push_back(dofToMatIndex.get(component, + dofMap[it.getDofIndex()].global)); else - sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex())); + sendGlobalDofs.push_back(dofToMatIndex.get(component, it.getDofIndex())); int rank = it.getRank(); if (meshLevel > 0) @@ -435,7 +416,7 @@ namespace AMDiS { stdMpi.send(rank, sendGlobalDofs); } - for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); + for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, componentSpaces[component]); !it.end(); it.nextRank()) { int rank = it.getRank(); if (meshLevel > 0) @@ -446,7 +427,7 @@ namespace AMDiS { stdMpi.startCommunication(); - for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); + for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, componentSpaces[component]); !it.end(); it.nextRank()) { int rank = it.getRank(); if (meshLevel > 0) @@ -457,9 +438,9 @@ namespace AMDiS { if (dofMap.count(it.getDofIndex())) { DegreeOfFreedom d = stdMpi.getRecvData(rank)[counter++]; if (globalIndex) - dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d); + dofToMatIndex.add(component, dofMap[it.getDofIndex()].global, d); else - dofToMatIndex.add(i, it.getDofIndex(), d); + dofToMatIndex.add(component, it.getDofIndex(), d); } } } @@ -469,7 +450,7 @@ namespace AMDiS { stdMpi.clear(); - for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]); + for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); !it.end(); it.nextRank()) { vector<DegreeOfFreedom> sendGlobalDofs; @@ -477,9 +458,10 @@ namespace AMDiS { if (dofMap.count(it.getDofIndex())) { if (globalIndex) { sendGlobalDofs.push_back(dofMap[it.getDofIndex()].global); - sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global)); + sendGlobalDofs.push_back(dofToMatIndex.get(component, + dofMap[it.getDofIndex()].global)); } else { - sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex())); + sendGlobalDofs.push_back(dofToMatIndex.get(component, it.getDofIndex())); } } @@ -493,7 +475,7 @@ namespace AMDiS { stdMpi.startCommunication(); - for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]); + for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); !it.end(); it.nextRank()) { int rank = it.getRank(); @@ -508,12 +490,13 @@ namespace AMDiS { TEST_EXIT_DBG(counter + 2 <= stdMpi.getRecvData(it.getRank()).size()) ("Should not happen!\n"); - dofToMatIndex.add(i, + dofToMatIndex.add(component, stdMpi.getRecvData(it.getRank())[counter], stdMpi.getRecvData(it.getRank())[counter + 1]); counter += 2; } else { - dofToMatIndex.add(i, it.getDofIndex(), + dofToMatIndex.add(component, + it.getDofIndex(), stdMpi.getRecvData(it.getRank())[counter++]); } } diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h index bad2f7371fa7ef7c2c4d0b8c12bcd1f41f32dc50..f38230954479676a5005bf6e0bca376726169c94 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.h +++ b/AMDiS/src/parallel/ParallelDofMapping.h @@ -222,17 +222,10 @@ namespace AMDiS { feSpace = fe; } - /// Informs the mapping whether the mapping will include DOFs that are not - /// owned by the rank. - void setNonLocal(bool b) - { - isNonLocal = b; - } - /// Informs the mapping whether a global index must be computed. - void setNeedGlobalMapping(bool b) + void setGlobalMapping(bool b) { - needGlobalMapping = b; + globalMapping = b; } /// Sets the DOF communicator. @@ -276,12 +269,7 @@ namespace AMDiS { boost::container::flat_set<DegreeOfFreedom> nonRankDofs; /// If true, a global index mapping will be computed for all DOFs. - bool needGlobalMapping; - - /// Is true if there are DOFs in at least one subdomain that are not owned - /// by the rank. If the value is false, each rank contains only DOFs that - /// are also owned by this rank. - bool isNonLocal; + bool globalMapping; public: /// @@ -289,6 +277,8 @@ namespace AMDiS { }; + /// Abstract iterator interface to access containrs containing values of + /// type \ref ComponentDofMap. class ComponentIterator { public: virtual ComponentDofMap& operator*() = 0; @@ -303,240 +293,249 @@ namespace AMDiS { }; + /// Abstract interface to acces DOF mapping data for each component. Allows + /// to hide specific implementations, which allow, e.g., to efficiently map + /// all components having the same FE space to the same DOF mapping. class ComponentDataInterface { public: + /// Access via component number virtual ComponentDofMap& operator[](int compNumber) = 0; + /// Acess via FE space pointer virtual ComponentDofMap& operator[](const FiniteElemSpace *feSpace) = 0; + /// Checks whether the DOF mapping is defined for a specific + /// component number. virtual bool isDefinedFor(int compNumber) const = 0; + /// Returns iterator which iterates over all DOF mappings. virtual ComponentIterator& getIteratorData() = 0; + /// Returns iterator which iterates over the DOF mappings of all + /// components. If the data is defined for each FE space and more than + /// one commponent is defined on the same FE space, the iterative will + /// also iterative multple times over the same DOF mapping object. virtual ComponentIterator& getIteratorComponent() = 0; - virtual void init(vector<const FiniteElemSpace*> &f0, - vector<const FiniteElemSpace*> &f1, - bool isNonLocal, + /** \brief + * Initialization of the object. + * + * \param[in] componentSpaces Set of the FE spaces for each component. + * \param[in] feSpaces Set of all different FE spaces. + * \param[in] globalMapping Mapping is parallel (true) or only + * local (false). + * \param[in] levelData Data for multi level method. + */ + virtual void init(vector<const FiniteElemSpace*> &componentSpaces, + vector<const FiniteElemSpace*> &feSpaces, + bool globalMapping, MeshLevelData &levelData) = 0; - vector<const FiniteElemSpace*>& getFeSpaces() - { - return feSpaces; - } - protected: /// The FE spaces for all components. - vector<const FiniteElemSpace*> feSpaces; + vector<const FiniteElemSpace*> componentSpaces; /// The set of all FE spaces. It uniquly contains all different FE spaces /// from \ref feSpaces. - vector<const FiniteElemSpace*> feSpacesUnique; + vector<const FiniteElemSpace*> feSpaces; }; - class ComponentDataEqFeSpace : public ComponentDataInterface + /// This class concretizes the interface class \ref ComponentDataInterface. A + /// DOF mapping is implemented for each component. + class ComponentData : public ComponentDataInterface { public: - ComponentDataEqFeSpace() - : iterData(this), - iterComponent(this) + ComponentData() + : iter(this) {} + /// Returns DOF mapping for a given component number. ComponentDofMap& operator[](int compNumber) { - const FiniteElemSpace *feSpace = feSpaces[compNumber]; - return componentData.find(feSpace)->second; - } + TEST_EXIT_DBG(componentData.count(compNumber)) + ("No data for component %d!\n", compNumber); - ComponentDofMap& operator[](const FiniteElemSpace *feSpace) - { - TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");; - - return componentData.find(feSpace)->second; + return componentData.find(compNumber)->second; } - bool isDefinedFor(int compNumber) const + /// Just to implement the corresponding virtual function in \ref + /// ComponentDataInterface. Of course it does not work as we have data for + /// each component. Thus there may be different mappings for the same + /// FE space. + ComponentDofMap& operator[](const FiniteElemSpace *feSpace) { - const FiniteElemSpace *feSpace = feSpaces[compNumber]; - return static_cast<bool>(componentData.count(feSpace)); + ERROR_EXIT("FE Space acces is not possible for component wise defined DOF mappings\n"); } - + + /// Return data iterator. ComponentIterator& getIteratorData() { - iterData.reset(); - return iterData; + iter.reset(); + return iter; } + /// Return component iterator. ComponentIterator& getIteratorComponent() { - iterComponent.reset(); - return iterComponent; + iter.reset(); + return iter; } + /// Checks whether the DOF mapping is defined for a specific + /// component number. + bool isDefinedFor(int compNumber) const + { + return (static_cast<unsigned int>(compNumber) < componentData.size()); + } + + /// Initialization void init(vector<const FiniteElemSpace*> &f0, vector<const FiniteElemSpace*> &f1, - bool isNonLocal, + bool globalMapping, MeshLevelData &levelData); - protected: - void addFeSpace(const FiniteElemSpace* feSpace, - MeshLevelData &levelData); - - class IteratorData : public ComponentIterator { + /// Iterator class to iterate over all parallel DOF mappings. + class Iterator : public ComponentIterator { public: - IteratorData(ComponentDataEqFeSpace *d) - : data(d) + Iterator(ComponentData *d) + : data(d), + componentCounter(-1) {} ComponentDofMap& operator*() { - (*data)[*it]; + (*data)[componentCounter]; } ComponentDofMap* operator->() { - &((*data)[*it]); + &((*data)[componentCounter]); } + bool end() { - return (it == data->feSpacesUnique.end()); + return (it == data->componentSpaces.end()); } void next() { ++it; + ++componentCounter; + + if (it == data->componentSpaces.end()) + componentCounter = -1; } void reset() { - it = data->feSpacesUnique.begin(); + it = data->componentSpaces.begin(); + componentCounter = 0; } protected: - ComponentDataEqFeSpace *data; + /// Pointer to data class over which the iterator must iterate. + ComponentData *data; + /// Internal iterator of the internal data from \ref ComponentData. vector<const FiniteElemSpace*>::iterator it; - }; - - class IteratorComponent : public ComponentIterator { - public: - IteratorComponent(ComponentDataEqFeSpace *d) - : data(d) - {} - - ComponentDofMap& operator*() - { - (*data)[*it]; - } - - ComponentDofMap* operator->() - { - &((*data)[*it]); - } - - bool end() - { - return (it == data->feSpaces.end()); - } - void next() - { - ++it; - } - - void reset() - { - it = data->feSpaces.begin(); - } + /// Component number of current iteration. + int componentCounter; + }; - protected: - ComponentDataEqFeSpace *data; - vector<const FiniteElemSpace*>::iterator it; - }; - + /// Data mapping from component numbers to DOF mapping objects. + map<unsigned int, ComponentDofMap> componentData; - map<const FiniteElemSpace*, ComponentDofMap> componentData; + /// Iterator object. + Iterator iter; - IteratorData iterData; - IteratorComponent iterComponent; - friend class IteratorData; - - friend class IteratorComponent; + friend class Iterator; }; - class ComponentDataDiffFeSpace : public ComponentDataInterface + + /// This class concretizes the interface class \ref ComponentDataInterface. A + /// DOF mapping is implemented for each finite element space. Thus, different + /// components sharing the same FE space are handled by the same DOF mapping. + class FeSpaceData : public ComponentDataInterface { public: - ComponentDataDiffFeSpace() - : iter(this) + FeSpaceData() + : iterData(this), + iterComponent(this) {} + /// Returns DOF mapping for a given component number. ComponentDofMap& operator[](int compNumber) { - TEST_EXIT_DBG(componentData.count(compNumber))("No data for component %d!\n", compNumber); - - return componentData.find(compNumber)->second; + const FiniteElemSpace *feSpace = componentSpaces[compNumber]; + return componentData.find(feSpace)->second; } + /// Returns DOF mapping for a given FE space. ComponentDofMap& operator[](const FiniteElemSpace *feSpace) { - ERROR_EXIT("FE Space acces is not possible for component wise defined DOF mappings\n"); + TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");; + + return componentData.find(feSpace)->second; } - + + /// Checks whether the DOF mapping is defined for a specific + /// component number. + bool isDefinedFor(int compNumber) const + { + const FiniteElemSpace *feSpace = componentSpaces[compNumber]; + return static_cast<bool>(componentData.count(feSpace)); + } + + /// Return data iterator. ComponentIterator& getIteratorData() { - iter.reset(); - return iter; + iterData.reset(); + return iterData; } + /// Return component iterator. ComponentIterator& getIteratorComponent() { - iter.reset(); - return iter; - } - - bool isDefinedFor(int compNumber) const - { - return (static_cast<unsigned int>(compNumber) < componentData.size()); + iterComponent.reset(); + return iterComponent; } + /// Initialization void init(vector<const FiniteElemSpace*> &f0, vector<const FiniteElemSpace*> &f1, - bool isNonLocal, + bool globalMapping, MeshLevelData &levelData); + protected: - void addComponent(unsigned int component, - const FiniteElemSpace* feSpace, - MeshLevelData &levelData); - class Iterator : public ComponentIterator { + /// Iterator class to iterate over all parallel DOF mappings. + class IteratorData : public ComponentIterator { public: - Iterator(ComponentDataDiffFeSpace *d) - : data(d), - componentCounter(-1) + IteratorData(FeSpaceData *d) + : data(d) {} ComponentDofMap& operator*() { - (*data)[componentCounter]; + (*data)[*it]; } ComponentDofMap* operator->() { - &((*data)[componentCounter]); + &((*data)[*it]); } - bool end() { return (it == data->feSpaces.end()); @@ -545,33 +544,72 @@ namespace AMDiS { void next() { ++it; - ++componentCounter; - - if (it == data->feSpaces.end()) - componentCounter = -1; } void reset() { it = data->feSpaces.begin(); - componentCounter = 0; } protected: - ComponentDataDiffFeSpace *data; + FeSpaceData *data; vector<const FiniteElemSpace*>::iterator it; + }; - int componentCounter; + + /// Iterator class to iterate over all component DOF mappings. + class IteratorComponent : public ComponentIterator { + public: + IteratorComponent(FeSpaceData *d) + : data(d) + {} + + ComponentDofMap& operator*() + { + (*data)[*it]; + } + + ComponentDofMap* operator->() + { + &((*data)[*it]); + } + + bool end() + { + return (it == data->componentSpaces.end()); + } + + void next() + { + ++it; + } + + void reset() + { + it = data->componentSpaces.begin(); + } + + protected: + FeSpaceData *data; + + vector<const FiniteElemSpace*>::iterator it; }; + - map<unsigned int, ComponentDofMap> componentData; + map<const FiniteElemSpace*, ComponentDofMap> componentData; - Iterator iter; + IteratorData iterData; - friend class Iterator; + IteratorComponent iterComponent; + + friend class IteratorData; + + friend class IteratorComponent; }; + + /// Used to specify whether a parallel DOF mapping is defined for each /// specific component or for each FE space. enum DofMappingMode { @@ -587,9 +625,16 @@ namespace AMDiS { class ParallelDofMapping { public: + /** \brief + * Constructur for parallel DOF mapping. + * + * \param[in] mode Defines if DOF mapping is defined either per + * component or per FE space. + */ ParallelDofMapping(DofMappingMode mode); - /** \brief Initialize the parallel DOF mapping. + /** \brief + * Initialize the parallel DOF mapping. * * \param[in] m MPI communicator. * \param[in] fe The FE spaces of all components of the @@ -597,22 +642,22 @@ namespace AMDiS { * \param[in] uniqueFe Unique list of FE spaces. Thus, two * arbitrary elements of this list are always * different. - * \param[in] isNonLocal If true, at least one rank's mapping con- + * \param[in] globalMapping If true, at least one rank's mapping con- * taines DOFs that are not owend by the rank. */ void init(MeshLevelData& mld, vector<const FiniteElemSpace*> &fe, vector<const FiniteElemSpace*> &uniqueFe, - bool isNonLocal = true); + bool globalMapping = true); /// In the case of having only one FE space, this init function can be used. void init(MeshLevelData& mld, const FiniteElemSpace *feSpace, - bool isNonLocal = true) + bool globalMapping = true) { vector<const FiniteElemSpace*> feSpaces; feSpaces.push_back(feSpace); - init(mld, feSpaces, feSpaces, isNonLocal); + init(mld, feSpaces, feSpaces, globalMapping); } void setMpiComm(MPI::Intracomm &m, int l); @@ -659,6 +704,8 @@ namespace AMDiS { return (*data)[feSpace]; } + /// Checks whether the DOF mapping is defined for a specific + /// component number. inline bool isDefinedFor(int compNumber) const { return data->isDefinedFor(compNumber); @@ -667,7 +714,7 @@ namespace AMDiS { /// Returns the number of solution components the mapping is defined on. inline int getNumberOfComponents() const { - return static_cast<int>(data->getFeSpaces().size()); + return static_cast<int>(componentSpaces.size()); } /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank. @@ -727,12 +774,6 @@ namespace AMDiS { return dofToMatIndex.get(ithComponent, d) - rStartDofs; } - // inline const FiniteElemSpace* getFeSpace(int i = 0) - // { - // TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n"); - // return feSpacesUnique[i]; - // } - // Writes all data of this object to an output stream. void serialize(ostream &out) { @@ -798,7 +839,7 @@ namespace AMDiS { /// Is true if there are DOFs in at least one subdomain that are not owned /// by the rank. If the value is false, each rank contains only DOFs that /// are also owned by this rank. - bool isNonLocal; + bool globalMapping; /// If matrix indices should be computed, this variable defines if the /// mapping from DOF indices to matrix row indices is defined on local @@ -829,6 +870,9 @@ namespace AMDiS { /// Mapping from global DOF indices to global matrix indices under /// consideration of possibly multiple components. DofToMatIndex dofToMatIndex; + + /// FE spaces of all components + vector<const FiniteElemSpace*> componentSpaces; }; } diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index d3ede6baab4813a23f796f601b9a0b732c9bd45a..0e224116aac59ec72f11adf53f8e04a76d22ffaf 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -142,26 +142,12 @@ namespace AMDiS { FUNCNAME("PetscSolverFeti::createDirichletData()"); int nComponents = mat.getSize(); - for (int i = 0; i < nComponents; i++) { - DOFMatrix* dofMat = mat[i][i]; + for (int component = 0; component < nComponents; component++) { + DOFMatrix* dofMat = mat[component][component]; if (!dofMat) continue; - - const FiniteElemSpace *feSpace = dofMat->getRowFeSpace(); - std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows(); - if (dirichletRows.count(feSpace)) { - // === Run test if Dirichlet rows are all the same in all components === - // === of one FE space. === - TEST_EXIT(dirichletRows[feSpace].size() == dRows.size()) - ("Wrong number of dirichlet rows!\n"); - - for (std::set<DegreeOfFreedom>::iterator it0 = dirichletRows[feSpace].begin(), - it1 = dRows.begin(); it1 != dRows.end(); ++it0, ++it1) { - TEST_EXIT(*it0 == *it1)("Wrong DOFs %d %d!\n", *it0, *it1); - } - } else { - dirichletRows[feSpace] = dRows; - } + + dirichletRows[component] = dofMat->getDirichletRows(); } } @@ -317,7 +303,7 @@ namespace AMDiS { for (DofContainerSet::iterator it = vertices.begin(); it != vertices.end(); ++it) { - if (dirichletRows[feSpace].count(**it)) + if (dirichletRows[component].count(**it)) continue; if (meshLevel == 0) { @@ -363,7 +349,7 @@ namespace AMDiS { for (DofContainer::iterator it = allBoundaryDofs.begin(); it != allBoundaryDofs.end(); ++it) { - if (dirichletRows[feSpace].count(**it)) + if (dirichletRows[component].count(**it)) continue; if (isPrimal(component, **it)) @@ -392,7 +378,7 @@ namespace AMDiS { for (DofContainer::iterator it = allBoundaryDofs.begin(); it != allBoundaryDofs.end(); ++it) { - if (dirichletRows[feSpace].count(**it)) + if (dirichletRows[component].count(**it)) continue; if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) @@ -572,7 +558,7 @@ namespace AMDiS { isPrimal(component, i) || isDual(component, i) || isInterface(component, i) || - dirichletRows[feSpace].count(i)) + dirichletRows[component].count(i)) continue; if (meshLevel == 0) { @@ -763,7 +749,7 @@ namespace AMDiS { for (DofContainer::iterator dit = edgeDofs.begin(); dit != edgeDofs.end(); ++dit) { - if (dirichletRows[feSpaces[0]].count(**dit) == 0) { + if (dirichletRows[0].count(**dit) == 0) { dirichletOnlyEdge = false; break; } @@ -803,7 +789,7 @@ namespace AMDiS { TEST_EXIT_DBG(isPrimal(component, **it) == false) ("Should not be primal!\n"); - if (dirichletRows[feSpaces[component]].count(**it)) + if (dirichletRows[component].count(**it)) continue; int col = lagrangeMap.getMatIndex(component, **it); @@ -1575,8 +1561,8 @@ namespace AMDiS { printDirichlet); if (printDirichlet) { int nComponents = mat->getSize(); - for (int i = 0; i < nComponents; i++) { - DOFMatrix *seqMat = (*mat)[i][i]; + for (int component = 0; component < nComponents; component++) { + DOFMatrix *seqMat = (*mat)[component][component]; if (!seqMat) continue; @@ -1585,9 +1571,9 @@ namespace AMDiS { std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows(); for (std::set<DegreeOfFreedom>::iterator it = dirichletRows.begin(); it != dirichletRows.end(); ++it) { - if (localDofMap[i].isSet(*it)) { + if (localDofMap[component].isSet(*it)) { MSG("Dirichlet dof %d in component %d with interior mat index %d\n", - *it, i, localDofMap.getMatIndex(i, *it)); + *it, component, localDofMap.getMatIndex(component, *it)); } } } @@ -1934,7 +1920,7 @@ namespace AMDiS { for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) { - if (dirichletRows[rowFeSpace].count(*cursor)) + if (dirichletRows[rowComponent].count(*cursor)) continue; if (isPrimal(rowComponent, *cursor)) @@ -1951,7 +1937,7 @@ namespace AMDiS { for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) { - if (dirichletRows[colFeSpace].count(col(*icursor))) + if (dirichletRows[colComponent].count(col(*icursor))) continue; if (isPrimal(colComponent, col(*icursor))) @@ -2016,7 +2002,7 @@ namespace AMDiS { // Traverse all columns. for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) { - if (dirichletRows[colFeSpace].count(col(*icursor))) + if (dirichletRows[colComponent].count(col(*icursor))) continue; if (!isDual(colComponent, col(*icursor))) diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index 6b5b74daeebdb1abaa6d11f48b31b72e9527fa9a..49b61e6f26b8d8018792ad82aeb87678f1c007a8 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -350,7 +350,9 @@ namespace AMDiS { int pressureComponent; - map<const FiniteElemSpace*, std::set<DegreeOfFreedom> > dirichletRows; + /// Maps from component number to set of DOFs which are Dirichlet DOfs in + /// this component. + map<int, std::set<DegreeOfFreedom> > dirichletRows; friend class PetscSolverFetiDebug; };