diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h index 3aa48152f144677ec6a68af75fe4ae03b0e8c970..f57dac7d0af36a736c07f10de02db13140187390 100644 --- a/AMDiS/src/DOFAdmin.h +++ b/AMDiS/src/DOFAdmin.h @@ -150,7 +150,7 @@ namespace AMDiS { return name; } - /// Returns \ref nDof[i], i.e., the number of dofs for the position i. + /// Returns \ref nDof[i], i.e., the number of DOFs for the position i. inline const int getNumberOfDofs(int i) const { return nDof[i]; diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h index 2d4be25fd2bbfa94dbfdf045806b67211362f422..9df5513e2fdfdde7d14d524282b3f6e1f80bf97b 100644 --- a/AMDiS/src/DOFMatrix.h +++ b/AMDiS/src/DOFMatrix.h @@ -492,7 +492,10 @@ namespace AMDiS { int nnzPerRow; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - std::map *rankDofs; + /// Stores for the DOFs of the row FE spaces whether they are owned by the + /// rank or not. This is used to ensure that Dirichlet BC is handled + /// correctly in parallel computations. + std::map *rankDofs; #endif /// Inserter object: implemented as pointer, allocated and deallocated as needed diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc index 305ae867f752a6cabea6395e1a0cfc6cbd08701b..50626931b430bbb56a2ea6ac52fd6b403f6414ac 100644 --- a/AMDiS/src/Debug.cc +++ b/AMDiS/src/Debug.cc @@ -25,7 +25,8 @@ namespace AMDiS { namespace debug { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace) + void writeLocalElementDofs(int rank, int elIdx, + const FiniteElemSpace *feSpace) { using boost::lexical_cast; @@ -37,7 +38,8 @@ namespace AMDiS { } - void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace) + void writeDofMesh(int rank, DegreeOfFreedom dof, + const FiniteElemSpace *feSpace) { using boost::lexical_cast; @@ -50,7 +52,7 @@ namespace AMDiS { } - void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename) + void writeMesh(const FiniteElemSpace *feSpace, int rank, std::string filename) { using boost::lexical_cast; @@ -64,7 +66,7 @@ namespace AMDiS { #endif - void writeDofIndexMesh(FiniteElemSpace *feSpace) + void writeDofIndexMesh(const FiniteElemSpace *feSpace) { DOFVector tmp(feSpace, "tmp"); DOFIterator it(&tmp, USED_DOFS); @@ -74,7 +76,7 @@ namespace AMDiS { } - void colorEdgeInMesh(FiniteElemSpace *feSpace, + void colorEdgeInMesh(const FiniteElemSpace *feSpace, Element *el, int localEdgeNo, std::string filename) @@ -190,7 +192,8 @@ namespace AMDiS { } - Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof) + Element* getDofIndexElement(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof) { const BasisFunction* basFcts = feSpace->getBasisFcts(); int nBasFcts = basFcts->getNumber(); @@ -322,7 +325,7 @@ namespace AMDiS { } - void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof) + void printInfoByDof(const FiniteElemSpace *feSpace, DegreeOfFreedom dof) { FUNCNAME("debug::printInfoByDof()"); @@ -387,7 +390,7 @@ namespace AMDiS { } - void printAllDofCoords(FiniteElemSpace *feSpace) + void printAllDofCoords(const FiniteElemSpace *feSpace) { FUNCNAME("printAllDofCoords()"); @@ -405,7 +408,8 @@ namespace AMDiS { } - void getAllDofs(FiniteElemSpace *feSpace, std::set& dofs) + void getAllDofs(const FiniteElemSpace *feSpace, + std::set& dofs) { FUNCNAME("getAllDofs()"); @@ -793,7 +797,7 @@ namespace AMDiS { } - void testDofsByCoords(FiniteElemSpace *feSpace, + void testDofsByCoords(const FiniteElemSpace *feSpace, DofContainer &dofs0, DofContainer &dofs1) { FUNCNAME("debug::testDofsByCoords()"); diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h index 3467f9b7716e4cbd2ed92a6601109b0350af5cdc..eeda743f3df715a601ca1c1e21ee7f0018c94be0 100644 --- a/AMDiS/src/Debug.h +++ b/AMDiS/src/Debug.h @@ -41,9 +41,13 @@ namespace AMDiS { typedef std::map ElementIdxToDofs; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace); + void writeLocalElementDofs(int rank, + int elIdx, + const FiniteElemSpace *feSpace); - void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename); + void writeMesh(const FiniteElemSpace *feSpace, + int rank, + std::string filename); /** \brief * Writes a vtu file with the mesh, where all DOFs are set to zero, and only @@ -55,7 +59,9 @@ namespace AMDiS { * \param[in] dof Defines the DOF, which value is set to one in the mesh file. * \param[in] feSpace The FE space to be used. */ - void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace); + void writeDofMesh(int rank, + DegreeOfFreedom dof, + const FiniteElemSpace *feSpace); #endif /** \brief @@ -64,9 +70,9 @@ namespace AMDiS { * * \param[in] feSpace The FE space to be used. */ - void writeDofIndexMesh(FiniteElemSpace *feSpace); + void writeDofIndexMesh(const FiniteElemSpace *feSpace); - void colorEdgeInMesh(FiniteElemSpace *feSpace, + void colorEdgeInMesh(const FiniteElemSpace *feSpace, Element *el, int localEdgeNo, std::string filename); @@ -92,7 +98,8 @@ namespace AMDiS { Mesh *mesh, int elIndex); - Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof); + Element* getDofIndexElement(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof); Element* getLevel0ParentElement(Mesh *mesh, Element *el); @@ -108,13 +115,14 @@ namespace AMDiS { void printElementCoords(const FiniteElemSpace *feSpace, Element *el); - void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof); + void printInfoByDof(const FiniteElemSpace *feSpace, DegreeOfFreedom dof); void printMatValuesStatistics(Matrix *mat); - void printAllDofCoords(FiniteElemSpace *feSpace); + void printAllDofCoords(const FiniteElemSpace *feSpace); - void getAllDofs(FiniteElemSpace *feSpace, std::set& dofs); + void getAllDofs(const FiniteElemSpace *feSpace, + std::set& dofs); /** \brief * Creates a text file storing the value of a sparse matrix. Each line of the file @@ -200,7 +208,7 @@ namespace AMDiS { const DegreeOfFreedom* dof3, DofContainer &vec); - void testDofsByCoords(FiniteElemSpace *feSpace, + void testDofsByCoords(const FiniteElemSpace *feSpace, DofContainer &dofs0, DofContainer &dofs1); } diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index 6091573c13e4261dfd3e2645c14f0a7228898d59..21e1e92cb9e88a8656e750f29ed2b2a3069bdd49 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -625,7 +625,7 @@ namespace AMDiS { } - void Element::getAllDofs(FiniteElemSpace* feSpace, + void Element::getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) { diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index edc2656473a606817ac53eb4560be32f6a95341e..b69e16062f6a34c5908133d2ab03f8107e3447aa 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -421,7 +421,7 @@ namespace AMDiS { * which all vertex dofs are assembled. * \param[out] dofs List of dofs, where the result is stored. */ - virtual void getNodeDofs(FiniteElemSpace* feSpace, + virtual void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const = 0; @@ -436,11 +436,11 @@ namespace AMDiS { * all non vertex dofs are assembled. * \param[out] dofs All dofs are put to this dof list. */ - virtual void getHigherOrderDofs(FiniteElemSpace* feSpace, + virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const = 0; - void getAllDofs(FiniteElemSpace* feSpace, + void getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs); diff --git a/AMDiS/src/ElementDofIterator.cc b/AMDiS/src/ElementDofIterator.cc index d9a34df8d7d83a7b6e7b32f9ec6388752df860b4..3da06de1c09d2658641c3a3c001c646cea05bb83 100644 --- a/AMDiS/src/ElementDofIterator.cc +++ b/AMDiS/src/ElementDofIterator.cc @@ -36,14 +36,15 @@ namespace AMDiS { // Get geo index of vertices in the given dimension. posIndex = INDEX_OF_DIM(pos, dim); - // Get number of dofs per vertex (should be one in all cases). + // Get number of DOFs per vertex (should be one in all cases). nDofs = admin->getNumberOfDofs(posIndex); TEST_EXIT_DBG(nDofs != 0)("Mh, I've to think about this situation!\n"); - // Calculate displacement. Is used if there is more than one dof admin on the mesh. + // Calculate displacement. Is used if there is more than one DOF admin + // on the mesh. n0 = admin->getNumberOfPreDofs(posIndex); - // Get first dof index position for vertices. + // Get first DOF index position for vertices. node0 = mesh->getNode(posIndex); // Get number of vertices in this dimension. nElements = Global::getGeo(posIndex, mesh->getDim()); @@ -55,11 +56,11 @@ namespace AMDiS { bool ElementDofIterator::next() { - // First iterate over the dofs of one element (vertex, edge, face). + // First iterate over the DOFs of one element (vertex, edge, face). dofPos++; if (dofPos >= nDofs) { - // We are finished with all dofs of on element. Go to the next one. + // We are finished with all DOFs of on element. Go to the next one. dofPos = 0; elementPos++; @@ -72,33 +73,35 @@ namespace AMDiS { return false; // Increase position, i.e., go from vertices to edges to faces and search - // for the next position with dofs. + // for the next position with DOFs. do { pos++; // Get geo index posistion. posIndex = INDEX_OF_DIM(pos, dim); - // Get number of dofs in this position. + // Get number of DOFs in this position. nDofs = admin->getNumberOfDofs(posIndex); } while (nDofs == 0 && pos < dim); if (nDofs > 0 && pos <= dim) { - // We have found on more position with dofs. + // We have found on more position with DOFs. - // Get number of elements in this position, i.e, the number of vertices,. - // edges and faces in the given dimension. + // Get number of elements in this position, i.e, the number of + // vertices, edges and faces in the given dimension. nElements = Global::getGeo(posIndex, dim); - // Calculate displacement. Is used if there is more than one dof admin on the mesh. + // Calculate displacement. Is used if there is more than one DOF + // admin on the mesh. n0 = admin->getNumberOfPreDofs(posIndex); - // Get first dof index position for the geo index position. + // Get first DOF index position for the geo index position. node0 = mesh->getNode(posIndex); if (inOrder) - orderPosition = basisFcts->orderOfPositionIndices(element, posIndex, 0); + orderPosition = + basisFcts->orderOfPositionIndices(element, posIndex, 0); } else { - // That's all, we jave traversed all dofs of the mesh element. + // That's all, we jave traversed all DOFs of the mesh element. return false; } diff --git a/AMDiS/src/ElementDofIterator.h b/AMDiS/src/ElementDofIterator.h index 90600367f4c749c80be17ae25692f4179abad3c1..8603956ab03af05a2bf7c001c08944fa9071c370 100644 --- a/AMDiS/src/ElementDofIterator.h +++ b/AMDiS/src/ElementDofIterator.h @@ -30,7 +30,7 @@ namespace AMDiS { /** \brief - * This class implements an iterator to iterate over all dofs of one element + * This class implements an iterator to iterate over all DOFs of one element * independet of dimension and the degree of basis functions. * * Should be used in the following way: @@ -57,12 +57,12 @@ namespace AMDiS { /// Start a new traverse with the given element. void reset(const Element* el); - /// Go to next dof. Returns false, if there is dof anymore. + /// Go to next DOF. Returns false, if there is no DOF anymore. bool next(); bool nextStrict(); - /// Returns the dof index of the current dof. + /// Returns the DOF index of the current DOF. inline DegreeOfFreedom getDof() { if (inOrder) @@ -71,7 +71,7 @@ namespace AMDiS { return dofs[node0 + elementPos][n0 + dofPos]; } - /// Returns a pointer to the current dof. + /// Returns a pointer to the current DOF. inline const DegreeOfFreedom* getDofPtr() { if (inOrder) @@ -80,13 +80,22 @@ namespace AMDiS { return &dofs[node0 + elementPos][n0 + dofPos]; } - /// Returns \ref pos, the current position (vertex, edge, face) of the traverse. + /// Returns a pointer to the starting position of the current DOF + /// array. Makes only sence, if \ref nextStrict() is used for traverse. + inline const DegreeOfFreedom* getStartDof() + { + return dofs[node0 + elementPos]; + } + + /// Returns \ref pos, the current position (vertex, edge, face) of + /// the traverse. inline int getCurrentPos() { return pos; } - /// Returns \ref elementPos, the number of vertex, edge or face that is traversed. + /// Returns \ref elementPos, the number of vertex, edge or face that + /// is traversed. inline int getCurrentElementPos() { return elementPos; @@ -99,12 +108,12 @@ namespace AMDiS { protected: - /// The dof admin for which dofs should be traversed. + /// The DOF admin for which DOFs should be traversed. const DOFAdmin* admin; const BasisFunction* basisFcts; - /// Pointer to the dofs that should be traversed. + /// Pointer to the DOFs that should be traversed. const DegreeOfFreedom **dofs; /// Mesh on which the element is defined. @@ -126,32 +135,28 @@ namespace AMDiS { GeoIndex posIndex; /** \brief - * Number of dofs at the current traverse position. Examples: independent of - * dimension and degree of basis functions there is only one dof per vertex. - * But in 2d and with 3rd degree lagrange basis functions there are two - * dofs per edge. + * Number of DOFs at the current traverse position. Examples: independent of + * dimension and degree of basis functions there is only one DOF per vertex. + * But in 2D and with 3rd degree lagrange basis functions there are two + * DOFs per edge. */ int nDofs; - /** \brief - * Displacement of dof indices. Used of more than one dof admin is defined - * on the mesh. - */ + /// Displacement of DOF indices. Used if more than one DOF admin is defined + /// on the mesh. int n0; - /// Dof index of the first dof at this geo index position. + /// DOF index of the first DOF at this geo index position. int node0; - /** \brief - * Number of elements in the current geo position. Examples: 3 vertices in 2d, - * 1 face in 2d, 4 faces in 3d, etc. - */ + /// Number of elements in the current geo position. Examples: 3 vertices in + /// 2d, 1 face in 2d, 4 faces in 3d, etc. int nElements; /// Current element, i.e., ith vertex, edge or face, that is traversed. int elementPos; - /// Currrent dof that is traversed on the current element; + /// Currrent DOF that is traversed on the current element; int dofPos; }; } diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index 3b7c58a45b036e6a18bc7d040cd7bfeb6a454385..2c8de3adbbf60a66293238f096496145912d779b 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -173,13 +173,13 @@ namespace AMDiS { return "Line"; } - void getNodeDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const + void getNodeDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const { FUNCNAME("Line::getNodeDofs()"); ERROR_EXIT("Not yet implemented!\n"); } - void getHigherOrderDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const + void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const { FUNCNAME("Line::getHigherOrderDofs()"); ERROR_EXIT("Not yet implemented!\n"); diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 101c7db2b30af5a9ee66d77dbb681681f34f0ed5..9f50b51c923b582252f9f74fb38385a4f51c8082 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -271,21 +271,32 @@ namespace AMDiS { void Mesh::removeMacroElements(std::set& macros, - const FiniteElemSpace *feSpace) + vector& feSpaces) { FUNCNAME("Mesh::removeMacroElement()"); typedef map > DofElMap; typedef map DofPosMap; - TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n"); - TEST_EXIT(admin[0])("There is something wrong!\n"); + TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); + + // === Search for the FE space with the highest degree of polynomials. === + // === Using this FE space ensures that deleting DOFs defined on it, === + // === also DOFs of lower order FE spaces will be deleted correctly. === + + const FiniteElemSpace *feSpace = feSpaces[0]; + for (unsigned int i = 1; i < feSpaces.size(); i++) + if (feSpaces[i]->getBasisFcts()->getDegree() > + feSpace->getBasisFcts()->getDegree()) + feSpace = feSpaces[i]; + + // === Determine to all DOFs in mesh the macro elements where the DOF === // === is part of. === - // Map that stores for each dof pointer (which may have a list of dofs) - // all macro element indices that own this dof. + // Map that stores for each DOF pointer (which may have a list of DOFs) + // all macro element indices that own this DOF. DofElMap dofsOwner; DofPosMap dofsPosIndex; @@ -295,8 +306,8 @@ namespace AMDiS { while (elInfo) { elDofIter.reset(elInfo->getElement()); do { - dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement()); - dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex(); + dofsOwner[elDofIter.getStartDof()].insert(elInfo->getMacroElement()); + dofsPosIndex[elDofIter.getStartDof()] = elDofIter.getPosIndex(); } while (elDofIter.nextStrict()); elInfo = stack.traverseNext(elInfo); @@ -615,15 +626,14 @@ namespace AMDiS { ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF); for (unsigned int i = 0; i < admin.size(); i++) { - DOFAdmin *localAdmin = admin[i]; - int n = localAdmin->getNumberOfDofs(position); - int n0 = localAdmin->getNumberOfPreDofs(position); + int n = admin[i]->getNumberOfDofs(position); + int n0 = admin[i]->getNumberOfPreDofs(position); TEST_EXIT_DBG(n + n0 <= ndof) ("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof); for (int j = 0; j < n; j++) - localAdmin->freeDofIndex(dof[n0 + j]); + admin[i]->freeDofIndex(dof[n0 + j]); } delete [] dof; @@ -974,7 +984,7 @@ namespace AMDiS { } - void Mesh::getAllDofs(FiniteElemSpace *feSpace, + void Mesh::getAllDofs(const FiniteElemSpace *feSpace, std::set& allDofs) { FUNCNAME("Mesh::getAllDofs()"); @@ -1381,7 +1391,7 @@ namespace AMDiS { FiniteElemSpace *feSpace = FiniteElemSpace::provideFeSpace(localAdmin, basFcts, &testMesh, "tmp"); - DataCollector dc(feSpace); + DataCollector<> dc(feSpace); MacroWriter::writeMacro(&dc, newMacroFilename.str().c_str()); if (periodicFilename != "") diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index c7f1d729bf6683725b877797eb9410ed6da0ddac..38996ede47ad92cf828ab0f67da6003bbffda97a 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -409,12 +409,12 @@ namespace AMDiS { void addMacroElement(MacroElement* me); /* \brief - * Removes a set of macro elements from the mesh. This works only for the case, - * that there are no global or local refinements, i.e., all macro elements have - * no children. + * Removes a set of macro elements from the mesh. This works only for the + * case, that there are no global or local refinements, i.e., all macro + * elements have no children. */ void removeMacroElements(std::set& macros, - const FiniteElemSpace* feSpace); + vector& feSpaces); /// Frees the array of DOF pointers (see \ref createDofPtrs) void freeDofPtrs(DegreeOfFreedom **ptrs); @@ -515,7 +515,7 @@ namespace AMDiS { * @param[in] feSpace The FE space to be used for collecting DOFs. * @param[out] allDofs The set which is filled with all DOFs. */ - void getAllDofs(FiniteElemSpace *feSpace, + void getAllDofs(const FiniteElemSpace *feSpace, std::set& allDofs); /// Returns FILL_ANY_?D diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h index 4c487044df1ebf351ebf3ebd7aef2afeb26011a5..b2861f1db1bc20ab2dcea4331a893d8ca3428c2e 100644 --- a/AMDiS/src/ProblemStat.h +++ b/AMDiS/src/ProblemStat.h @@ -350,7 +350,7 @@ namespace AMDiS { return rhs; } - inline const DOFVector* getRhsVector(int i = 0) + inline DOFVector* getRhsVector(int i = 0) { return rhs->getDOFVector(i); } diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc index 2320e4081816c3bdd61a85e6d132c5f1c91b74c9..50f56de365c9094b4ba428de484cecb6baf11354 100644 --- a/AMDiS/src/SecondOrderTerm.cc +++ b/AMDiS/src/SecondOrderTerm.cc @@ -434,7 +434,8 @@ namespace AMDiS { FUNCNAME("MatrixGradient_SOT::weakEval()"); TEST_EXIT_DBG(f)("No function f!\n"); - TEST_EXIT_DBG(gradAtQPs)("Operator was not initialized correctly!\n"); + TEST_EXIT_DBG(num_rows(gradAtQPs)) + ("Operator was not initialized correctly!\n"); int nPoints = grdUhAtQP.size(); for (int iq = 0; iq < nPoints; iq++) { diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index 605c8fa769ecd296588ed1c00cc6a87e20a0770f..cf8985f7c649bada0c77ee545dacc04424e221bf 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -200,7 +200,7 @@ namespace AMDiS { } - void Tetrahedron::getNodeDofs(FiniteElemSpace* feSpace, + void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { @@ -222,7 +222,7 @@ namespace AMDiS { } - void Tetrahedron::getNodeDofsAtFace(FiniteElemSpace* feSpace, + void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { @@ -285,7 +285,7 @@ namespace AMDiS { } - void Tetrahedron::getNodeDofsAtEdge(FiniteElemSpace* feSpace, + void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { @@ -359,7 +359,7 @@ namespace AMDiS { } - void Tetrahedron::getHigherOrderDofs(FiniteElemSpace* feSpace, + void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index 07984d3c8091da53721522633db203710b72ed57..492c928736b2ec29cc36f734639c5b88e60e1c71 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -136,19 +136,19 @@ namespace AMDiS { return "Tetrahedron"; } - void getNodeDofs(FiniteElemSpace* feSpace, + void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; - void getNodeDofsAtFace(FiniteElemSpace* feSpace, + void getNodeDofsAtFace(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; - void getNodeDofsAtEdge(FiniteElemSpace* feSpace, + void getNodeDofsAtEdge(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; - void getHigherOrderDofs(FiniteElemSpace* feSpace, + void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc index e7ee0bc22bcb1ad33d79947e74b24c9c627eed0a..179aa908da8ea883e01d22aed250edb6a2b294b0 100644 --- a/AMDiS/src/Triangle.cc +++ b/AMDiS/src/Triangle.cc @@ -87,7 +87,7 @@ namespace AMDiS { } - void Triangle::getNodeDofs(FiniteElemSpace* feSpace, + void Triangle::getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { @@ -111,10 +111,10 @@ namespace AMDiS { child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); dofs.push_back(child[1]->getFirstChild()->getDof(2)); nextBound.ithObj = 0; - child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); } else { nextBound.ithObj = 0; - child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); dofs.push_back(child[1]->getFirstChild()->getDof(2)); nextBound.ithObj = 1; child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); @@ -164,7 +164,7 @@ namespace AMDiS { } - void Triangle::getHigherOrderDofs(FiniteElemSpace* feSpace, + void Triangle::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const { diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index 86e8c3b694e6b72db55da74d3d795c65f66a92eb..fa1658e18026c159e0b87c513665414477971c70 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -194,11 +194,11 @@ namespace AMDiS { return "Triangle"; } - void getNodeDofs(FiniteElemSpace* feSpace, + void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; - void getHigherOrderDofs(FiniteElemSpace* feSpace, + void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs) const; diff --git a/AMDiS/src/io/DataCollector.hh b/AMDiS/src/io/DataCollector.hh index 594ad36d1eeca15840ffafae544bb4d81097e6ca..f581acd1dde4984520b4c9496f6232462a89c432 100644 --- a/AMDiS/src/io/DataCollector.hh +++ b/AMDiS/src/io/DataCollector.hh @@ -541,4 +541,4 @@ namespace AMDiS { } } -#endif \ No newline at end of file +#endif diff --git a/AMDiS/src/io/VtkVectorWriter.hh b/AMDiS/src/io/VtkVectorWriter.hh index a0653b8d880893ec1db0522441116e551d3ee5ee..fe438692b1295412b1dd25e26896eb392e8be995 100644 --- a/AMDiS/src/io/VtkVectorWriter.hh +++ b/AMDiS/src/io/VtkVectorWriter.hh @@ -483,4 +483,4 @@ namespace AMDiS { } -#endif // AMDIS_VTKVECTORWRITER_HH \ No newline at end of file +#endif // AMDIS_VTKVECTORWRITER_HH diff --git a/AMDiS/src/parallel/ElementObjectData.cc b/AMDiS/src/parallel/ElementObjectData.cc index 075814aa0b7c055d8166596e9103257c9e8e25d3..7dd8f0495540048f667f6783b1b66153ed5fd58b 100644 --- a/AMDiS/src/parallel/ElementObjectData.cc +++ b/AMDiS/src/parallel/ElementObjectData.cc @@ -125,7 +125,7 @@ namespace AMDiS { } - void ElementObjects::createPeriodicData() + void ElementObjects::createPeriodicData(const FiniteElemSpace *feSpace) { FUNCNAME("ElementObjects::createPeriodicData()"); @@ -216,9 +216,9 @@ namespace AMDiS { ("Should not happen!\n"); periodicVertices[make_pair(dof0, dof3)] = - provideConnectedPeriodicBoundary(type0, type1); + provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1); periodicVertices[make_pair(dof3, dof0)] = - provideConnectedPeriodicBoundary(type0, type1); + provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1); for (unsigned int j = i + 1; j < multPeriodicDof2.size(); j++) if (multPeriodicDof2[j] == dof3) @@ -237,7 +237,7 @@ namespace AMDiS { make_pair(multPeriodicDof3[j], multPeriodicDof3[i]); if (periodicVertices.count(perDofs0) == 0) { - BoundaryType b = getNewBoundaryType(); + BoundaryType b = getNewBoundaryType(feSpace->getAdmin()); periodicVertices[perDofs0] = b; periodicVertices[perDofs1] = b; } @@ -268,7 +268,8 @@ namespace AMDiS { BoundaryType type0 = periodicEdges[make_pair(edge0, edge1)]; BoundaryType type1 = periodicEdges[make_pair(edge0, edge2)]; - BoundaryType type2 = provideConnectedPeriodicBoundary(type0, type1); + BoundaryType type2 = + provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1); periodicEdges[perEdge0] = type2; periodicEdges[perEdge1] = type2; @@ -311,7 +312,7 @@ namespace AMDiS { } - BoundaryType ElementObjects::getNewBoundaryType() + BoundaryType ElementObjects::getNewBoundaryType(DOFAdmin *admin) { FUNCNAME("ElementObjects::getNewBoundaryType()"); @@ -324,13 +325,14 @@ namespace AMDiS { newPeriodicBoundaryType--; mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = - new VertexVector(feSpace->getAdmin(), ""); + new VertexVector(admin, ""); return newPeriodicBoundaryType; } - BoundaryType ElementObjects::provideConnectedPeriodicBoundary(BoundaryType b0, + BoundaryType ElementObjects::provideConnectedPeriodicBoundary(DOFAdmin *admin, + BoundaryType b0, BoundaryType b1) { FUNCNAME("ElementObjects::provideConnectedPeriodicBoundary()"); @@ -339,13 +341,13 @@ namespace AMDiS { (b0 <= b1 ? make_pair(b0, b1) : make_pair(b1, b0)); if (bConnMap.count(bConn) == 0) { - BoundaryType newPeriodicBoundaryType = getNewBoundaryType(); + BoundaryType newPeriodicBoundaryType = getNewBoundaryType(admin); VertexVector &vecB0 = mesh->getPeriodicAssociations(b0); VertexVector &vecB1 = mesh->getPeriodicAssociations(b1); VertexVector &vecC = mesh->getPeriodicAssociations(newPeriodicBoundaryType); - DOFIteratorBase it(const_cast(feSpace->getAdmin()), USED_DOFS); + DOFIteratorBase it(const_cast(admin), USED_DOFS); for (it.reset(); !it.end(); ++it) { if (!it.isDofFree()) { @@ -417,7 +419,8 @@ namespace AMDiS { } - void ElementObjects::createReverseModeData(map &elIndexMap, + void ElementObjects::createReverseModeData(const FiniteElemSpace* feSpace, + map &elIndexMap, map &elIndexTypeMap) { FUNCNAME("ElementObjects::createReverseModeData()"); @@ -497,7 +500,8 @@ namespace AMDiS { EDGE, edges1[j].ithObject); bool reverseMode = - BoundaryObject::computeReverseMode(obj0, obj1, feSpace, edgeIt->second); + BoundaryObject::computeReverseMode(obj0, obj1, feSpace, + edgeIt->second); edgeReverseMode[make_pair(edges0[i], edges1[j])] = reverseMode; edgeReverseMode[make_pair(edges1[j], edges0[i])] = reverseMode; diff --git a/AMDiS/src/parallel/ElementObjectData.h b/AMDiS/src/parallel/ElementObjectData.h index 8cfaddad580179e8ca0c11832354aa4ab162e251..8bc930706ca1738a7e8c2265dda462cca7bd0ea0 100644 --- a/AMDiS/src/parallel/ElementObjectData.h +++ b/AMDiS/src/parallel/ElementObjectData.h @@ -101,18 +101,15 @@ namespace AMDiS { class ElementObjects { public: ElementObjects() - : feSpace(NULL), - mesh(NULL), + : mesh(NULL), iterGeoPos(CENTER) {} - /// Set the finite element space that should be used for the database (especially - /// the mesh is of interest). - void setFeSpace(FiniteElemSpace *fe) + /// Set the mesh that should be used for the database. + void setMesh(Mesh *m) { - feSpace = fe; - mesh = feSpace->getMesh(); + mesh = m; } @@ -135,7 +132,7 @@ namespace AMDiS { * box, all eight vertex nodes and always four of the 12 edges are iterectly * connected. */ - void createPeriodicData(); + void createPeriodicData(const FiniteElemSpace *feSpace); /** \brief @@ -148,7 +145,19 @@ namespace AMDiS { */ void createRankData(map& macroElementRankMap); - void createReverseModeData(map &elIndexMap, + /** \brief + * Creates on all boundaries the reverse mode flag. + * + * \param[in] feSpace An arbitrary FE space defined on the mesh. + * Is used to get the orientation of the DOFs on + * elements. + * \param[in] elIndexMap Maps an element index to the pointer to the + * element. + * \param[in] elIndexTypeMap Maps an element index to its type id (only + * relevant in 3D). + */ + void createReverseModeData(const FiniteElemSpace* feSpace, + map &elIndexMap, map &elIndexTypeMap); @@ -469,10 +478,11 @@ namespace AMDiS { faceLocalMap[elObj] = face; } - BoundaryType provideConnectedPeriodicBoundary(BoundaryType b0, - BoundaryType b1); + BoundaryType getNewBoundaryType(DOFAdmin *admin); - BoundaryType getNewBoundaryType(); + BoundaryType provideConnectedPeriodicBoundary(DOFAdmin *admin, + BoundaryType b0, + BoundaryType b1); /// Some auxiliary function to write the element object database to disk. void serialize(ostream &out, vector& elVec); @@ -487,13 +497,9 @@ namespace AMDiS { void deserialize(istream &in, map& data); private: - /// The used FE space. - FiniteElemSpace *feSpace; - /// The mesh that is used to store all its element information in the database. Mesh *mesh; - /// Maps to each vertex DOF all element objects that represent this vertex. map > vertexElements; diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc index 2413e3f9f267cb2ee869691dfa2063e6a850c584..bd776e170dc1a32b13651c420164d523751e1888 100644 --- a/AMDiS/src/parallel/InteriorBoundary.cc +++ b/AMDiS/src/parallel/InteriorBoundary.cc @@ -42,7 +42,7 @@ namespace AMDiS { bool BoundaryObject::computeReverseMode(BoundaryObject &obj0, BoundaryObject &obj1, - FiniteElemSpace *feSpace, + const FiniteElemSpace *feSpace, BoundaryType boundary) { FUNCNAME("BoundaryObject::computeReverseMode()"); diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h index dff95391436d4d22bed164336e375c5ab62f08c0..ed46d8670a26f7949ada81b515bae06927c7233f 100644 --- a/AMDiS/src/parallel/InteriorBoundary.h +++ b/AMDiS/src/parallel/InteriorBoundary.h @@ -48,7 +48,7 @@ namespace AMDiS { static bool computeReverseMode(BoundaryObject &obj0, BoundaryObject &obj1, - FiniteElemSpace *feSpace, + const FiniteElemSpace *feSpace, BoundaryType boundary); bool operator==(const BoundaryObject& other) const; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 0145162c60a9820f60b31649376e495ee7c3db6a..bc8f7efe87e2dc898aa09ca710bc2b95325e44ed 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -71,14 +71,11 @@ namespace AMDiS { : problemStat(0), initialized(false), name("parallel"), - feSpace(NULL), + feSpaces(0), mesh(NULL), refineManager(NULL), info(10), partitioner(NULL), - nRankDofs(0), - rStartDofs(0), - nOverallDofs(0), deserialized(false), writeSerializationFile(false), repartitioningAllowed(false), @@ -144,10 +141,27 @@ namespace AMDiS { TEST_EXIT(mpiSize > 1) ("Parallelization does not work with only one process!\n"); - TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n"); + TEST_EXIT(feSpaces.size() > 0) + ("No FE space has been defined for the mesh distributor!\n"); TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n"); - elObjects.setFeSpace(feSpace); + // === Sort FE spaces with respect to the degree of the basis === + // === functions. Use stuiped bubble sort for this. === + + bool doNext = false; + do { + for (unsigned int i = 0; i < feSpaces.size() - 1; i++) + if (feSpaces[i]->getBasisFcts()->getDegree() > + feSpaces[i + 1]->getBasisFcts()->getDegree()) { + const FiniteElemSpace *tmp = feSpaces[i + 1]; + feSpaces[i + 1] = feSpaces[i]; + feSpaces[i] = tmp; + doNext = true; + } + } while (doNext); + + + elObjects.setMesh(feSpaces[0]->getMesh()); // If the problem has been already read from a file, we need only to set // isRankDofs to all matrices and rhs vector and to remove periodic @@ -299,14 +313,14 @@ namespace AMDiS { ParallelDebug::testInteriorBoundary(*this); ParallelDebug::followBoundary(*this); - debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh"); + debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh"); MSG("Debug mode tests finished!\n"); #endif // === Create periodic DOF mapping, if there are periodic boundaries. === - createPeriodicMap(); + createPeriodicMap(feSpaces[0]); #if (DEBUG != 0) ParallelDebug::testPeriodicBoundary(*this); @@ -325,7 +339,7 @@ namespace AMDiS { // === Update periodic mapping, if there are periodic boundaries. === - createPeriodicMap(); + createPeriodicMap(feSpaces[0]); #if (DEBUG != 0) ParallelDebug::testPeriodicBoundary(*this); @@ -350,40 +364,45 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::addProblemStat()"); - if (feSpace != NULL) { - vector feSpaces = probStat->getFeSpaces(); - for (unsigned int i = 0; i < feSpaces.size(); i++) { -// MSG("MESH %p <-> %p BF %p <-> %p\n", -// feSpace->getMesh(), -// feSpaces[i]->getMesh(), -// feSpace->getBasisFcts(), -// feSpaces[i]->getBasisFcts()); - TEST_EXIT(feSpace == feSpaces[i]) - ("Parallelizaton is not supported for multiple FE spaces!\n"); - } - } else { - feSpace = probStat->getFeSpace(0); - mesh = feSpace->getMesh(); - info = probStat->getInfo(); - - TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1) - ("Only meshes with one DOFAdmin are supported!\n"); - TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0) - ("Wrong pre dof number for DOFAdmin!\n"); - - switch (mesh->getDim()) { - case 2: - refineManager = new RefinementManager2d(); - break; - case 3: - refineManager = new RefinementManager3d(); - break; - default: - ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim()); + // === Add FE spaces from stationary problem to mesh distributor. === + + for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) { + bool foundFeSpace = false; + for (unsigned int j = 0; j < feSpaces.size(); j++) { + if (feSpaces[j] == probStat->getFeSpaces()[i]) + foundFeSpace = true; + + TEST_EXIT(feSpaces[j]->getMesh() == probStat->getFeSpaces()[i]->getMesh()) + ("MeshDistributor does not yet support usage of different meshes!\n"); } - partitioner->setMesh(mesh); + if (!foundFeSpace) + feSpaces.push_back(probStat->getFeSpaces()[i]); + } + + TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); + + mesh = feSpaces[0]->getMesh(); + info = probStat->getInfo(); + + TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1) + ("Only meshes with one DOFAdmin are supported!\n"); + TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0) + ("Wrong pre dof number for DOFAdmin!\n"); + + switch (mesh->getDim()) { + case 2: + refineManager = new RefinementManager2d(); + break; + case 3: + refineManager = new RefinementManager3d(); + break; + default: + ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim()); } + + partitioner->setMesh(mesh); + // Create parallel serialization file writer, if needed. int writeSerialization = 0; @@ -502,40 +521,45 @@ namespace AMDiS { void MeshDistributor::synchVector(SystemVector &vec) { + FUNCNAME("MeshDistributor::synchVector()"); + int nComponents = vec.getSize(); StdMpi > stdMpi(mpiComm); - for (RankToDofContainer::iterator sendIt = sendDofs.begin(); - sendIt != sendDofs.end(); ++sendIt) { + typedef map >::iterator it_type; + + for (it_type sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt) { vector dofs; - int nSendDofs = sendIt->second.size(); - dofs.reserve(nComponents * nSendDofs); - + for (int i = 0; i < nComponents; i++) { - DOFVector *dofvec = vec.getDOFVector(i); - for (int j = 0; j < nSendDofs; j++) - dofs.push_back((*dofvec)[*((sendIt->second)[j])]); - } + DofContainer &feDofs = sendIt->second[vec.getFeSpace(i)]; + DOFVector& dofVec = *(vec.getDOFVector(i)); + + int nFeDofs = feDofs.size(); + dofs.resize(dofs.size() + nFeDofs); + + for (int j = 0; j < nFeDofs; j++) + dofs.push_back(dofVec[*(feDofs[j])]); + } stdMpi.send(sendIt->first, dofs); } - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt) - stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents); + for (it_type recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) + stdMpi.recv(recvIt->first); stdMpi.startCommunication(); - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt) { - int nRecvDofs = recvIt->second.size(); + for (it_type recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) { int counter = 0; for (int i = 0; i < nComponents; i++) { - DOFVector *dofvec = vec.getDOFVector(i); - for (int j = 0; j < nRecvDofs; j++) - (*dofvec)[*(recvIt->second)[j]] = - stdMpi.getRecvData(recvIt->first)[counter++]; + DofContainer &feDofs = recvIt->second[vec.getFeSpace(i)]; + DOFVector& dofVec = *(vec.getDOFVector(i)); + int nFeDofs = feDofs.size(); + + for (int j = 0; j < nFeDofs; j++) + dofVec[*(feDofs[j])] = stdMpi.getRecvData(recvIt->first)[counter++]; } } } @@ -634,8 +658,8 @@ namespace AMDiS { DofEdge dofEdge1 = edge1.first->getEdge(edge1.second); WorldVector c0, c1; - mesh->getDofIndexCoords(dofEdge0.first, feSpace, c0); - mesh->getDofIndexCoords(dofEdge0.second, feSpace, c1); + mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0); + mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1); MSG_DBG("FOUND EDGE: %d %d with coords %f %f %f %f %f %f\n", dofEdge0.first, dofEdge0.second, @@ -655,23 +679,20 @@ namespace AMDiS { } - void MeshDistributor::getAllBoundaryDofs(DofContainer& dofs) + void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace, + DofContainer& dofs) { FUNCNAME("MeshDistributor::getAllBoundaryDofs()"); DofContainerSet dofSet; - for (RankToDofContainer::iterator sendIt = sendDofs.begin(); - sendIt != sendDofs.end(); ++sendIt) - for (DofContainer::iterator dofIt = sendIt->second.begin(); - dofIt != sendIt->second.end(); ++dofIt) - dofSet.insert(*dofIt); + typedef map >::iterator it_type; - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt) - for (DofContainer::iterator dofIt = recvIt->second.begin(); - dofIt != recvIt->second.end(); ++dofIt) - dofSet.insert(*dofIt); + for (it_type it = sendDofs.begin(); it != sendDofs.end(); ++it) + dofSet.insert(it->second[feSpace].begin(), it->second[feSpace].end()); + + for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it) + dofSet.insert(it->second[feSpace].begin(), it->second[feSpace].end()); dofs.clear(); dofs.insert(dofs.begin(), dofSet.begin(), dofSet.end()); @@ -680,18 +701,35 @@ namespace AMDiS { void MeshDistributor::setRankDofs(ProblemStatSeq *probStat) { + FUNCNAME("MeshDistributor::setRankDofs()"); + int nComponents = probStat->getNumComponents(); for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) - if (probStat->getSystemMatrix(i, j)) - probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof); + if (probStat->getSystemMatrix(i, j)) { + const FiniteElemSpace *rowFeSpace = + probStat->getSystemMatrix(i, j)->getRowFeSpace(); + + TEST_EXIT(find(feSpaces.begin(), feSpaces.end(), rowFeSpace) != feSpaces.end()) + ("Should not happen!\n"); + + probStat->getSystemMatrix(i, j)->setRankDofs(dofFeData[rowFeSpace].isRankDof); + } + TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n"); TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i)) ("No solution vector!\n"); + TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() == + probStat->getSolution(i)->getFeSpace()) + ("Should not happen!\n"); - probStat->getRhs()->getDOFVector(i)->setRankDofs(isRankDof); - probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof); + const FiniteElemSpace *feSpace = probStat->getRhsVector(i)->getFeSpace(); + TEST_EXIT(find(feSpaces.begin(), feSpaces.end(), feSpace) != feSpaces.end()) + ("Should not happen!\n"); + + probStat->getRhsVector(i)->setRankDofs(dofFeData[feSpace].isRankDof); + probStat->getSolution(i)->setRankDofs(dofFeData[feSpace].isRankDof); } } @@ -793,7 +831,7 @@ namespace AMDiS { MeshStructure elCode; elCode.init(it->rankObj); - MeshManipulation mm(feSpace); + MeshManipulation mm(mesh); meshChanged |= mm.fitElementToMeshCode(elCode, it->neighObj); } } else { @@ -819,7 +857,7 @@ namespace AMDiS { } while (recvAllValues != 0); #if (DEBUG != 0) - debug::writeMesh(feSpace, -1, debugOutputDir + "mesh"); + debug::writeMesh(feSpaces[0], -1, debugOutputDir + "mesh"); #endif @@ -830,7 +868,7 @@ namespace AMDiS { // === Update periodic mapping, if there are periodic boundaries. === - createPeriodicMap(); + createPeriodicMap(feSpaces[0]); #if (DEBUG != 0) ParallelDebug::testPeriodicBoundary(*this); @@ -942,7 +980,7 @@ namespace AMDiS { // it->first); TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); - MeshManipulation mm(feSpace); + MeshManipulation mm(mesh); meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj); } } @@ -952,22 +990,23 @@ namespace AMDiS { } - void MeshDistributor::createBoundaryDofs(std::set &boundaryDofs) + void MeshDistributor::createBoundaryDofs(const FiniteElemSpace *feSpace, + std::set &boundaryDofs) { FUNCNAME("MeshDistributor::createBoundaryDofs()"); boundaryDofs.clear(); - for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + typedef map >::iterator it_type; + + for (it_type it = sendDofs.begin(); it != sendDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].begin(); ++dofIt) boundaryDofs.insert(**dofIt); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].begin(); ++dofIt) boundaryDofs.insert(**dofIt); } @@ -983,6 +1022,22 @@ namespace AMDiS { } + void MeshDistributor::serialize(ostream &out, + map > &data) + { + int mapSize = data.size(); + SerUtil::serialize(out, mapSize); + for (map >::iterator it = data.begin(); it != data.end(); ++it) { + int rank = it->first; + SerUtil::serialize(out, rank); + + for (map::iterator dcIt = it->second.begin(); + dcIt != it->second.end(); ++dcIt) + serialize(out, dcIt->second); + } + } + + void MeshDistributor::deserialize(istream &in, DofContainer &data, map &dofMap) { @@ -1004,20 +1059,9 @@ namespace AMDiS { } - void MeshDistributor::serialize(ostream &out, RankToDofContainer &data) - { - int mapSize = data.size(); - SerUtil::serialize(out, mapSize); - for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) { - int rank = it->first; - SerUtil::serialize(out, rank); - serialize(out, it->second); - } - } - - - void MeshDistributor::deserialize(istream &in, RankToDofContainer &data, - map &dofMap) + void MeshDistributor::deserialize(istream &in, + map > &data, + map > &dofMap) { data.clear(); @@ -1026,7 +1070,9 @@ namespace AMDiS { for (int i = 0; i < mapSize; i++) { int rank = 0; SerUtil::deserialize(in, rank); - deserialize(in, data[rank], dofMap); + + for (unsigned int j = 0; j < feSpaces.size(); j++) + deserialize(in, data[rank][feSpaces[j]], dofMap[feSpaces[j]]); } } @@ -1116,7 +1162,7 @@ namespace AMDiS { if (repartitioningCounter == 0) ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", - repartitioningCounter, feSpace); + repartitioningCounter, feSpaces[0]); #endif repartitioningCounter++; @@ -1132,7 +1178,7 @@ namespace AMDiS { // === Run mesh partitioner to calculate a new mesh partitioning. === - partitioner->setLocalGlobalDofMap(&mapLocalGlobalDofs); + // partitioner->setLocalGlobalDofMap(&mapLocalGlobalDofs); bool partitioningSucceed = partitioner->partition(elemWeights, ADAPTIVE_REPART); if (!partitioningSucceed) { MSG("Mesh partitioner created empty partition!\n"); @@ -1309,12 +1355,14 @@ namespace AMDiS { // Note that also if there are no macros to be deleted, this function will // update the number of elements, vertices, etc. of the mesh. - mesh->removeMacroElements(deleteMacroElements, feSpace); + mesh->removeMacroElements(deleteMacroElements, feSpaces); // === Remove double DOFs. === - MeshManipulation meshManipulation(feSpace); - meshManipulation.deleteDoubleDofs(newMacroEl, elObjects); + TEST_EXIT(feSpaces.size() == 1) + ("Repartitioning not yet implemented for multiple FE spaces!\n"); + MeshManipulation meshManipulation(mesh); + meshManipulation.deleteDoubleDofs(feSpaces[0], newMacroEl, elObjects); mesh->dofCompress(); partitioner->createPartitionMap(partitionMap); @@ -1324,14 +1372,14 @@ namespace AMDiS { // === Update periodic mapping, if there are periodic boundaries. === - createPeriodicMap(); + createPeriodicMap(feSpaces[0]); #if (DEBUG != 0) MSG("AMDiS runs in debug mode, so make some test ...\n"); ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", - repartitioningCounter, feSpace); + repartitioningCounter, feSpaces[0]); ParallelDebug::testAllElements(*this); ParallelDebug::testDoubleDofs(mesh); ParallelDebug::testInteriorBoundary(*this); @@ -1419,10 +1467,13 @@ namespace AMDiS { // === Create periodic data, if there are periodic boundary conditions. === - elObjects.createPeriodicData(); + TEST_EXIT(feSpaces.size() == 1) + ("Sebastian: Na, dass funktioniert auch noch nicht mit mehreren FE spaces. Du weisst schon, wen du jetzt mobben kannst :)!\n"); + elObjects.createPeriodicData(feSpaces[0]); // === Create data about the reverse modes of neighbouring elements. === - elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap); + elObjects.createReverseModeData(feSpaces[0], macroElIndexMap, + macroElIndexTypeMap); // === Create mesh element data for this rank. === elObjects.createRankData(partitionMap); @@ -1526,9 +1577,10 @@ namespace AMDiS { if (elObjects.isInRank(it->first.first, mpiRank) == false) continue; + TEST_EXIT(feSpaces.size() == 1)("Does not work for multiple fe spaces!\n"); WorldVector c0, c1; - mesh->getDofIndexCoords(it->first.first, feSpace, c0); - mesh->getDofIndexCoords(it->first.second, feSpace, c1); + mesh->getDofIndexCoords(it->first.first, feSpaces[0], c0); + mesh->getDofIndexCoords(it->first.second, feSpaces[0], c1); ElementObjectData& perDofEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank]; @@ -1750,68 +1802,87 @@ namespace AMDiS { sendDofs.clear(); recvDofs.clear(); + for (unsigned int i = 0; i < feSpaces.size(); i++) + createBoundaryDofs(feSpaces[i]); + } + + + void MeshDistributor::createBoundaryDofs(const FiniteElemSpace *feSpace) + { + FUNCNAME("MeshDistributor::createBoundaryDofs()"); + if (createBoundaryDofFlag.isSet(BOUNDARY_SUBOBJ_SORTED)) { + + // === Clear data. === for (int geo = FACE; geo >= VERTEX; geo--) - boundaryDofInfo.geoDofs[static_cast(geo)].clear(); + boundaryDofInfo[feSpace].geoDofs[static_cast(geo)].clear(); + // === Create send DOFs. === for (int geo = FACE; geo >= VERTEX; geo--) { for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { if (it->rankObj.subObj == geo) { DofContainer dofs; it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs); - DofContainer& tmp = sendDofs[it.getRank()]; + DofContainer& tmp = sendDofs[it.getRank()][feSpace]; tmp.insert(tmp.end(), dofs.begin(), dofs.end()); if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_SEND_DOFS)) - boundaryDofInfo.geoDofs[static_cast(geo)].insert(dofs.begin(), dofs.end()); + boundaryDofInfo[feSpace].geoDofs[static_cast(geo)].insert(dofs.begin(), dofs.end()); } } } + // === Create recv DOFs. === for (int geo = FACE; geo >= VERTEX; geo--) { for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { if (it->rankObj.subObj == geo) { DofContainer dofs; it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs); - DofContainer& tmp = recvDofs[it.getRank()]; + DofContainer& tmp = recvDofs[it.getRank()][feSpace]; tmp.insert(tmp.end(), dofs.begin(), dofs.end()); if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_RECV_DOFS)) - boundaryDofInfo.geoDofs[static_cast(geo)].insert(dofs.begin(), dofs.end()); + boundaryDofInfo[feSpace].geoDofs[static_cast(geo)].insert(dofs.begin(), dofs.end()); } } } } else { for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) it->rankObj.el->getAllDofs(feSpace, it->rankObj, - sendDofs[it.getRank()]); + sendDofs[it.getRank()][feSpace]); for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) it->rankObj.el->getAllDofs(feSpace, it->rankObj, - recvDofs[it.getRank()]); + recvDofs[it.getRank()][feSpace]); } // === Delete all empty DOF send and recv positions === - { - RankToDofContainer::iterator it = sendDofs.begin(); - while (it != sendDofs.end()) { - if (it->second.size() == 0) - sendDofs.erase(it++); - else + typedef map >::iterator it_type; + for (it_type sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt) { + map::iterator it = + sendIt->second.begin(); + while (it != sendIt->second.end()) { + if (it->second.size() == 0) { + const FiniteElemSpace* fe = it->first; + ++it; + sendIt->second.erase(fe); + } else ++it; } } - - { - RankToDofContainer::iterator it = recvDofs.begin(); - while (it != recvDofs.end()) { - if (it->second.size() == 0) - recvDofs.erase(it++); - else + for (it_type recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) { + map::iterator it = + recvIt->second.begin(); + while (it != recvIt->second.end()) { + if (it->second.size() == 0) { + const FiniteElemSpace* fe = it->first; + ++it; + recvIt->second.erase(fe); + } else ++it; } } @@ -1828,7 +1899,7 @@ namespace AMDiS { if (partitioner->getElementInRank()[(*it)->getIndex()] == false) macrosToRemove.insert(*it); - mesh->removeMacroElements(macrosToRemove, feSpace); + mesh->removeMacroElements(macrosToRemove, feSpaces); } @@ -1836,13 +1907,52 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::updateLocalGlobalNumbering()"); - mesh->dofCompress(); - #if (DEBUG != 0) debug::ElementIdxToDofs elMap; debug::createSortedDofs(mesh, elMap); #endif + sendDofs.clear(); + recvDofs.clear(); + + for (unsigned int i = 0; i < feSpaces.size(); i++) + updateLocalGlobalNumbering(feSpaces[i]); + + + +#if (DEBUG != 0) + MSG("------------- Debug information -------------\n"); + for (unsigned int i = 0; i < feSpaces.size(); i++) { + MSG("FE space %d:\n", i); + MSG(" nRankDofs = %d\n", dofFeData[feSpaces[i]].nRankDofs); + MSG(" nOverallDofs = %d\n", dofFeData[feSpaces[i]].nOverallDofs); + MSG(" rStartDofs %d\n", dofFeData[feSpaces[i]].rStartDofs); + } + + stringstream oss; + oss << debugOutputDir << "elementIndex-" << mpiRank << ".vtu"; + debug::writeElementIndexMesh(mesh, oss.str()); + + ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat"); + + debug::testSortedDofs(mesh, elMap); + ParallelDebug::testCommonDofs(*this, true); + ParallelDebug::testGlobalIndexByCoords(*this); + +#else + int tmp = 0; + Parameters::get(name + "->write parallel debug file", tmp); + if (tmp) + ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat"); +#endif + } + + + void MeshDistributor::updateLocalGlobalNumbering(const FiniteElemSpace *feSpace) + { + FUNCNAME("MeshDistributor::updateLocalGlobalNumbering()"); + + mesh->dofCompress(); // === Get all DOFs in ranks partition. === @@ -1855,14 +1965,15 @@ namespace AMDiS { // === Traverse interior boundaries and get all DOFs on them. === - createBoundaryDofs(); + createBoundaryDofs(feSpace); // All DOFs that must be received are DOFs not owned by rank and have // therefore to be removed from the set 'rankDofs'. - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); + for (map >::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) { - for (DofContainer::iterator dofIt = recvIt->second.begin(); - dofIt != recvIt->second.end(); ++dofIt) { + DofContainer &rDofs = recvIt->second[feSpace]; + for (DofContainer::iterator dofIt = rDofs.begin(); + dofIt != rDofs.end(); ++dofIt) { DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *dofIt); if (eraseIt != rankDofs.end()) @@ -1872,14 +1983,17 @@ namespace AMDiS { // Get displacment for global rank DOF ordering and global DOF number. - nRankDofs = rankDofs.size(); - mpi::getDofNumbering(mpiComm, nRankDofs, rStartDofs, nOverallDofs); + dofFeData[feSpace].nRankDofs = rankDofs.size(); + mpi::getDofNumbering(mpiComm, + dofFeData[feSpace].nRankDofs, + dofFeData[feSpace].rStartDofs, + dofFeData[feSpace].nOverallDofs); // Stores for all rank owned DOFs a new global index. DofIndexMap rankDofsNewGlobalIndex; - for (int i = 0; i < nRankDofs; i++) - rankDofsNewGlobalIndex[rankDofs[i]] = i + rStartDofs; + for (int i = 0; i < dofFeData[feSpace].nRankDofs; i++) + rankDofsNewGlobalIndex[rankDofs[i]] = i + dofFeData[feSpace].rStartDofs; // === Send and receive new DOF indices. === @@ -1888,90 +2002,75 @@ namespace AMDiS { ParallelDebug::testDofContainerCommunication(*this, sendDofs, recvDofs); #endif - StdMpi > stdMpi(mpiComm, false); - for (RankToDofContainer::iterator sendIt = sendDofs.begin(); + StdMpi > stdMpi(mpiComm); + for (map >::iterator sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt) { + DofContainer &sDofs = sendIt->second[feSpace]; stdMpi.getSendData(sendIt->first).resize(0); - stdMpi.getSendData(sendIt->first).reserve(sendIt->second.size()); - for (DofContainer::iterator dofIt = sendIt->second.begin(); - dofIt != sendIt->second.end(); ++dofIt) + stdMpi.getSendData(sendIt->first).reserve(sDofs.size()); + for (DofContainer::iterator dofIt = sDofs.begin(); + dofIt != sDofs.end(); ++dofIt) stdMpi.getSendData(sendIt->first).push_back(rankDofsNewGlobalIndex[*dofIt]); } stdMpi.updateSendDataSize(); - stdMpi.recv(recvDofs); + + for (map >::iterator recvIt = recvDofs.begin(); + recvIt != recvDofs.end(); ++recvIt) + stdMpi.recv(recvIt->first); + stdMpi.startCommunication(); // First, we set all DOFs in ranks partition to be owend by the rank. Than, // the DOFs in ranks partition that are owned by other rank are set to false. - isRankDof.clear(); + dofFeData[feSpace].isRankDof.clear(); for (int i = 0; i < nRankAllDofs; i++) - isRankDof[i] = true; + dofFeData[feSpace].isRankDof[i] = true; - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt) { + for (map >::iterator recvIt = recvDofs.begin(); + recvIt != recvDofs.end(); ++recvIt) { + DofContainer &rDofs = recvIt->second[feSpace]; int i = 0; - for (DofContainer::iterator dofIt = recvIt->second.begin(); - dofIt != recvIt->second.end(); ++dofIt) { + for (DofContainer::iterator dofIt = rDofs.begin(); + dofIt != rDofs.end(); ++dofIt) { rankDofsNewGlobalIndex[*dofIt] = stdMpi.getRecvData(recvIt->first)[i++]; - isRankDof[**dofIt] = false; + dofFeData[feSpace].isRankDof[**dofIt] = false; } } // === Create now the local to global index and local to DOF index mappings. === - mapLocalGlobalDofs.clear(); - mapLocalDofIndex.clear(); + dofFeData[feSpace].mapLocalGlobalDofs.clear(); + dofFeData[feSpace].mapLocalDofIndex.clear(); for (DofIndexMap::iterator dofIt = rankDofsNewGlobalIndex.begin(); dofIt != rankDofsNewGlobalIndex.end(); ++dofIt) - mapLocalGlobalDofs[*(dofIt->first)] = dofIt->second; + dofFeData[feSpace].mapLocalGlobalDofs[*(dofIt->first)] = dofIt->second; - for (int i = 0; i < nRankDofs; i++) - mapLocalDofIndex[i] = *(rankDofs[i]); + for (int i = 0; i < dofFeData[feSpace].nRankDofs; i++) + dofFeData[feSpace].mapLocalDofIndex[i] = *(rankDofs[i]); // === Update dof admins due to new number of DOFs. === lastMeshChangeIndex = mesh->getChangeIndex(); - - -#if (DEBUG != 0) - MSG("------------- Debug information -------------\n"); - MSG("nRankDofs = %d\n", nRankDofs); - MSG("nOverallDofs = %d\n", nOverallDofs); - MSG("rStartDofs %d\n", rStartDofs); - - stringstream oss; - oss << debugOutputDir << "elementIndex-" << mpiRank << ".vtu"; - debug::writeElementIndexMesh(mesh, oss.str()); - - ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat"); - - debug::testSortedDofs(mesh, elMap); - ParallelDebug::testCommonDofs(*this, true); - ParallelDebug::testGlobalIndexByCoords(*this); - -#else - int tmp = 0; - Parameters::get(name + "->write parallel debug file", tmp); - if (tmp) - ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat"); -#endif } - void MeshDistributor::createPeriodicMap() + void MeshDistributor::createPeriodicMap(const FiniteElemSpace *feSpace) { FUNCNAME("MeshDistributor::createPeriodicMap()"); if (periodicBoundary.boundary.size() == 0) return; + TEST_EXIT(feSpaces.size() > 1) + ("Nanu, schon wieder Arbeit fuer den Thomas! Auch das hier muss geaendert werden fuer den Einsatz von unterschiedlichen FE Dingern!\n"); + // Clear all periodic dof mappings calculated before. We do it from scratch. periodicDof.clear(); periodicDofAssociations.clear(); @@ -2007,8 +2106,10 @@ namespace AMDiS { BoundaryType type = bound.type; for (unsigned int j = 0; j < dofs0.size(); j++) { - DegreeOfFreedom globalDof0 = mapLocalGlobalDofs[*(dofs0[j])]; - DegreeOfFreedom globalDof1 = mapLocalGlobalDofs[*(dofs1[j])]; + DegreeOfFreedom globalDof0 = + dofFeData[feSpace].mapLocalGlobalDofs[*(dofs0[j])]; + DegreeOfFreedom globalDof1 = + dofFeData[feSpace].mapLocalGlobalDofs[*(dofs1[j])]; if (periodicDofAssociations[globalDof0].count(type) == 0) { periodicDof[type][globalDof0] = globalDof1; @@ -2033,7 +2134,7 @@ namespace AMDiS { // Send the global indices to the rank on the other side. stdMpi.getSendData(it->first).reserve(dofs.size()); for (unsigned int i = 0; i < dofs.size(); i++) - stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]); + stdMpi.getSendData(it->first).push_back(dofFeData[feSpace].mapLocalGlobalDofs[*(dofs[i])]); // Receive from this rank the same number of dofs. stdMpi.recv(it->first, dofs.size()); @@ -2059,7 +2160,7 @@ namespace AMDiS { // Added the received dofs to the mapping. for (unsigned int i = 0; i < dofs.size(); i++) { - int globalDofIndex = mapLocalGlobalDofs[*(dofs[i])]; + int globalDofIndex = dofFeData[feSpace].mapLocalGlobalDofs[*(dofs[i])]; int mapGlobalDofIndex = stdMpi.getRecvData(it->first)[i]; BoundaryType type = types[i]; @@ -2092,7 +2193,8 @@ namespace AMDiS { boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs); for (unsigned int i = 0; i < dofs.size(); i++) { - DegreeOfFreedom globalDof = mapLocalGlobalDofs[*dofs[i]]; + DegreeOfFreedom globalDof = + dofFeData[feSpace].mapLocalGlobalDofs[*dofs[i]]; TEST_EXIT_DBG(periodicDofAssociations.count(globalDof)) ("Should hot happen!\n"); @@ -2166,12 +2268,13 @@ namespace AMDiS { } - DegreeOfFreedom MeshDistributor::mapGlobalToLocal(DegreeOfFreedom dof) + DegreeOfFreedom MeshDistributor::mapGlobalToLocal(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof) { FUNCNAME("MeshDistributor::mapGlobalToLocal()"); - for (DofMapping::iterator it = mapLocalGlobalDofs.begin(); - it != mapLocalGlobalDofs.end(); ++it) + for (DofMapping::iterator it = dofFeData[feSpace].mapLocalGlobalDofs.begin(); + it != dofFeData[feSpace].mapLocalGlobalDofs.end(); ++it) if (it->second == dof) return it->first; @@ -2189,8 +2292,6 @@ namespace AMDiS { SerUtil::serialize(out, elemWeights); SerUtil::serialize(out, partitionMap); - SerUtil::serialize(out, nRankDofs); - SerUtil::serialize(out, nOverallDofs); elObjects.serialize(out); @@ -2201,14 +2302,25 @@ namespace AMDiS { serialize(out, sendDofs); serialize(out, recvDofs); - SerUtil::serialize(out, mapLocalGlobalDofs); - SerUtil::serialize(out, mapLocalDofIndex); - SerUtil::serialize(out, isRankDof); + // === Serialieze FE space dependent data === + + unsigned int nFeSpace = feSpaces.size(); + SerUtil::serialize(out, nFeSpace); + + for (unsigned int i = 0; i < nFeSpace; i++) { + SerUtil::serialize(out, dofFeData[feSpaces[i]].nRankDofs); + SerUtil::serialize(out, dofFeData[feSpaces[i]].nOverallDofs); + SerUtil::serialize(out, dofFeData[feSpaces[i]].rStartDofs); + + SerUtil::serialize(out, dofFeData[feSpaces[i]].isRankDof); + SerUtil::serialize(out, dofFeData[feSpaces[i]].mapLocalGlobalDofs); + SerUtil::serialize(out, dofFeData[feSpaces[i]].mapLocalDofIndex); + } + serialize(out, periodicDof); serialize(out, periodicDofAssociations); - SerUtil::serialize(out, rStartDofs); SerUtil::serialize(out, macroElementNeighbours); int nSize = allMacroElements.size(); @@ -2229,28 +2341,28 @@ namespace AMDiS { SerUtil::deserialize(in, elemWeights); SerUtil::deserialize(in, partitionMap); - SerUtil::deserialize(in, nRankDofs); - SerUtil::deserialize(in, nOverallDofs); // Create two maps: one from from element indices to the corresponding element // pointers, and one map from Dof indices to the corresponding dof pointers. map elIndexMap; - map dofMap; - ElementDofIterator elDofIter(feSpace); - TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); - while (elInfo) { - Element *el = elInfo->getElement(); - elIndexMap[el->getIndex()] = el; - - if (el->isLeaf()) { - elDofIter.reset(el); - do { - dofMap[elDofIter.getDof()] = elDofIter.getDofPtr(); - } while (elDofIter.next()); + map > dofMap; + for (unsigned int i = 0; i < feSpaces.size(); i++) { + ElementDofIterator elDofIter(feSpaces[i]); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + Element *el = elInfo->getElement(); + elIndexMap[el->getIndex()] = el; + + if (el->isLeaf()) { + elDofIter.reset(el); + do { + dofMap[feSpaces[i]][elDofIter.getDof()] = elDofIter.getDofPtr(); + } while (elDofIter.next()); + } + + elInfo = stack.traverseNext(elInfo); } - - elInfo = stack.traverseNext(elInfo); } elObjects.deserialize(in); @@ -2262,14 +2374,28 @@ namespace AMDiS { deserialize(in, sendDofs, dofMap); deserialize(in, recvDofs, dofMap); - SerUtil::deserialize(in, mapLocalGlobalDofs); - SerUtil::deserialize(in, mapLocalDofIndex); - SerUtil::deserialize(in, isRankDof); + // === Deerialieze FE space dependent data === + + unsigned int nFeSpace = 0; + SerUtil::deserialize(in, nFeSpace); + TEST_EXIT(nFeSpace == feSpaces.size()) + ("Number of FE spaces is %d, but %d are stored in serialization file!\n", + feSpaces.size(), nFeSpace); + + for (unsigned int i = 0; i < nFeSpace; i++) { + SerUtil::deserialize(in, dofFeData[feSpaces[i]].nRankDofs); + SerUtil::deserialize(in, dofFeData[feSpaces[i]].nOverallDofs); + SerUtil::deserialize(in, dofFeData[feSpaces[i]].rStartDofs); + + SerUtil::deserialize(in, dofFeData[feSpaces[i]].isRankDof); + SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapLocalGlobalDofs); + SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapLocalDofIndex); + } + deserialize(in, periodicDof); deserialize(in, periodicDofAssociations); - SerUtil::deserialize(in, rStartDofs); SerUtil::deserialize(in, macroElementNeighbours); int nSize = 0; diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 63d3491000eb024e06ad4d64b55c959bf4bbca58..a5b24cb9b5476c8fe8a77a62b5b7573f9f3d2fb7 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -25,6 +25,7 @@ #include +#include "parallel/ElementObjectData.h" #include "parallel/ParallelTypes.h" #include "parallel/MeshPartitioner.h" #include "parallel/InteriorBoundary.h" @@ -36,7 +37,6 @@ #include "FiniteElemSpace.h" #include "Serializer.h" #include "BoundaryManager.h" -#include "ElementObjectData.h" #include "SystemVector.h" namespace AMDiS { @@ -48,8 +48,33 @@ namespace AMDiS { { map geoDofs; }; - + + struct DofData + { + /// Number of DOFs in the rank mesh. + int nRankDofs; + + /// Is the index of the first global DOF index, which is owned by the rank. + int rStartDofs; + + /// Number of DOFs in the whole domain. + int nOverallDofs; + + /** \brief + * Maps all DOFs in ranks partition to a bool value. If it is true, the DOF + * is owned by the rank. Otherwise, its an interior boundary DOF that is + * owned by another rank. + */ + DofIndexToBool isRankDof; + + /// Maps local to global dof indices. + DofMapping mapLocalGlobalDofs; + + /// Maps local dof indices to real dof indices. + DofMapping mapLocalDofIndex; + }; + class MeshDistributor { private: @@ -124,47 +149,60 @@ namespace AMDiS { return mesh; } - /// Returns \ref feSpace. - inline const FiniteElemSpace* getFeSpace() + /// Returns an FE space from \ref feSpaces. + inline const FiniteElemSpace* getFeSpace(unsigned int i = 0) { - return feSpace; + FUNCNAME("MeshDistributor::getFeSpace()"); + + TEST_EXIT_DBG(i < feSpaces.size())("Should not happen!\n"); + + return feSpaces[i]; } - + + /// Returns all FE spaces, thus \ref feSpaces. + inline vector& getFeSpaces() + { + return feSpaces; + } + /// Returns \ref nRankDOFs, the number of DOFs in the rank mesh. - inline int getNumberRankDofs() + inline int getNumberRankDofs(const FiniteElemSpace *feSpace) { - return nRankDofs; + return dofFeData[feSpace].nRankDofs; } /// Returns \ref rStartDofs, the first global DOF index owned by rank. - inline int getStartDofs() + inline int getStartDofs(const FiniteElemSpace *feSpace) { - return rStartDofs; + return dofFeData[feSpace].rStartDofs; } /// Returns \ref nOverallDofs, the global number of DOFs. - inline int getNumberOverallDofs() + inline int getNumberOverallDofs(const FiniteElemSpace *feSpace) { - return nOverallDofs; + return dofFeData[feSpace].nOverallDofs; } - inline DofMapping& getMapLocalGlobalDofs() + inline DofMapping& getMapLocalGlobalDofs(const FiniteElemSpace *feSpace) { - return mapLocalGlobalDofs; + return dofFeData[feSpace].mapLocalGlobalDofs; } /// Maps a local DOF to its global index. - inline DegreeOfFreedom mapLocalToGlobal(DegreeOfFreedom dof) + inline DegreeOfFreedom mapLocalToGlobal(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof) { - return mapLocalGlobalDofs[dof]; + return dofFeData[feSpace].mapLocalGlobalDofs[dof]; } - DegreeOfFreedom mapGlobalToLocal(DegreeOfFreedom dof); + DegreeOfFreedom mapGlobalToLocal(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof); /// Maps a local DOF to its local index. - inline DegreeOfFreedom mapLocalToDofIndex(DegreeOfFreedom dof) + inline DegreeOfFreedom mapLocalToDofIndex(const FiniteElemSpace *feSpace, + DegreeOfFreedom dof) { - return mapLocalDofIndex[dof]; + return dofFeData[feSpace].mapLocalDofIndex[dof]; } /// Returns the periodic mapping for all boundary DOFs in rank. @@ -211,19 +249,29 @@ namespace AMDiS { return (periodicDof[type].count(globalDofIndex) > 0); } + map >& getSendDofs() + { + return sendDofs; + } + + map >& getRecvDofs() + { + return recvDofs; + } + /// Return true, if the given DOF is owned by the rank. If false, the DOF /// is in rank's partition, but is owned by some other rank. - inline bool getIsRankDof(DegreeOfFreedom dof) + inline bool getIsRankDof(const FiniteElemSpace *feSpace, DegreeOfFreedom dof) { - if (isRankDof.count(dof)) - return isRankDof[dof]; + if (dofFeData[feSpace].isRankDof.count(dof)) + return dofFeData[feSpace].isRankDof[dof]; return false; } - inline DofIndexToBool& getIsRankDof() + inline DofIndexToBool& getIsRankDof(const FiniteElemSpace *feSpace) { - return isRankDof; + return dofFeData[feSpace].isRankDof; } inline long getLastMeshChangeIndex() @@ -246,19 +294,10 @@ namespace AMDiS { return mpiComm; } - inline RankToDofContainer& getSendDofs() - { - return sendDofs; - } - - inline RankToDofContainer& getRecvDofs() - { - return recvDofs; - } - /// Creates a set of all DOFs that are on interior boundaries of rank's /// domain. Thus, it creates the union of \ref sendDofs and \ref recvDofs. - void createBoundaryDofs(std::set &boundaryDofs); + void createBoundaryDofs(const FiniteElemSpace *feSpace, + std::set &boundaryDofs); // Writes all data of this object to an output stream. void serialize(ostream &out); @@ -281,28 +320,32 @@ namespace AMDiS { { StdMpi > stdMpi(mpiComm); - for (RankToDofContainer::iterator sendIt = sendDofs.begin(); + const FiniteElemSpace *fe = vec.getFeSpace(); + + typedef map >::iterator it_type; + + for (it_type sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt) { vector dofs; - int nSendDofs = sendIt->second.size(); + int nSendDofs = sendIt->second[fe].size(); dofs.reserve(nSendDofs); for (int i = 0; i < nSendDofs; i++) - dofs.push_back(vec[*((sendIt->second)[i])]); + dofs.push_back(vec[*((sendIt->second[fe])[i])]); stdMpi.send(sendIt->first, dofs); } - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); + for (it_type recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) - stdMpi.recv(recvIt->first, recvIt->second.size()); + stdMpi.recv(recvIt->first, recvIt->second[fe].size()); stdMpi.startCommunication(); - for (RankToDofContainer::iterator recvIt = recvDofs.begin(); + for (it_type recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) for (unsigned int i = 0; i < recvIt->second.size(); i++) - vec[*(recvIt->second)[i]] = stdMpi.getRecvData(recvIt->first)[i]; + vec[*(recvIt->second[fe])[i]] = stdMpi.getRecvData(recvIt->first)[i]; } /** \brief @@ -319,12 +362,13 @@ namespace AMDiS { createBoundaryDofFlag = flag; } - BoundaryDofInfo& getBoundaryDofInfo() + BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace) { - return boundaryDofInfo; + return boundaryDofInfo[feSpace]; } - void getAllBoundaryDofs(DofContainer& dofs); + void getAllBoundaryDofs(const FiniteElemSpace *feSpace, + DofContainer& dofs); public: @@ -347,20 +391,24 @@ namespace AMDiS { void createBoundaryDofs(); + void createBoundaryDofs(const FiniteElemSpace *feSpace); + /// Removes all macro elements from the mesh that are not part of ranks /// partition. void removeMacroElements(); + void updateLocalGlobalNumbering(); + /// Updates the local and global DOF numbering after the mesh has been /// changed. - void updateLocalGlobalNumbering(); + void updateLocalGlobalNumbering(const FiniteElemSpace *feSpace); /** \brief * Creates to all dofs in rank's partition that are on a periodic boundary * the mapping from dof index to the other periodic dof indices. This * information is stored in \ref periodicDof. */ - void createPeriodicMap(); + void createPeriodicMap(const FiniteElemSpace *feSpace); /** \brief * This function is called only once during the initialization when the @@ -409,16 +457,18 @@ namespace AMDiS { /// Writes a vector of dof pointers to an output stream. void serialize(ostream &out, DofContainer &data); + /// Writes a \ref RankToDofContainer to an output stream. + void serialize(ostream &out, + map > &data); + /// Reads a vector of dof pointers from an input stream. void deserialize(istream &in, DofContainer &data, map &dofMap); - /// Writes a \ref RankToDofContainer to an output stream. - void serialize(ostream &out, RankToDofContainer &data); - /// Reads a \ref RankToDofContainer from an input stream. - void deserialize(istream &in, RankToDofContainer &data, - map &dofMap); + void deserialize(istream &in, + map > &data, + map > &dofMap); /// Writes a periodic dof mapping to an output stream. void serialize(ostream &out, PeriodicDofMap &data); @@ -489,8 +539,8 @@ namespace AMDiS { /// Name of the problem (as used in the init files) string name; - /// Finite element space of the problem. - FiniteElemSpace *feSpace; + /// Finite element spaces of the problem. + vector feSpaces; /// Mesh of the problem. Mesh *mesh; @@ -519,14 +569,7 @@ namespace AMDiS { */ map partitionMap; - /// Number of DOFs in the rank mesh. - int nRankDofs; - - /// Is the index of the first global DOF index, which is owned by the rank. - int rStartDofs; - - /// Number of DOFs in the whole domain. - int nOverallDofs; + map dofFeData; /// Data structure to store all sub-objects of all elements of the /// macro mesh. @@ -564,27 +607,16 @@ namespace AMDiS { * This map contains for each rank the list of DOFs the current rank must * send to exchange solution DOFs at the interior boundaries. */ - RankToDofContainer sendDofs; + // map sendDofs; + map > sendDofs; /** \brief * This map contains on each rank the list of DOFs from which the current * rank will receive DOF values (i.e., this are all DOFs at an interior * boundary). The DOF indices are given in rank's local numbering. */ - RankToDofContainer recvDofs; - - /// Maps local to global dof indices. - DofMapping mapLocalGlobalDofs; - - /// Maps local dof indices to real dof indices. - DofMapping mapLocalDofIndex; - - /** \brief - * Maps all DOFs in ranks partition to a bool value. If it is true, the DOF - * is owned by the rank. Otherwise, its an interior boundary DOF that is - * owned by another rank. - */ - DofIndexToBool isRankDof; + // map recvDofs; + map > recvDofs; /** \brief * If periodic boundaries are used, this map stores, for each periodic @@ -655,7 +687,7 @@ namespace AMDiS { Flag createBoundaryDofFlag; - BoundaryDofInfo boundaryDofInfo; + map boundaryDofInfo; public: /// The boundary DOFs are sorted by subobject entities, i.e., first all diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc index 256cfa2fd5bda66c027f30a51317e64c63491ffd..d992c72aa95580aa5415dbcf32d9454a3b8b5f10 100644 --- a/AMDiS/src/parallel/MeshManipulation.cc +++ b/AMDiS/src/parallel/MeshManipulation.cc @@ -19,9 +19,8 @@ namespace AMDiS { - MeshManipulation::MeshManipulation(FiniteElemSpace *space) - : feSpace(space), - mesh(space->getMesh()) + MeshManipulation::MeshManipulation(Mesh *m) + : mesh(m) { switch (mesh->getDim()) { case 2: @@ -42,7 +41,8 @@ namespace AMDiS { } - void MeshManipulation::deleteDoubleDofs(std::set& newMacroEl, + void MeshManipulation::deleteDoubleDofs(const FiniteElemSpace *feSpace, + std::set& newMacroEl, ElementObjects &objects) { FUNCNAME("MeshManipulation::deleteDoubleDofs()"); @@ -248,10 +248,10 @@ namespace AMDiS { // s0 and s1 are the number of the edge/face in both child of the element, // which contain the edge/face the function has to traverse through. If the // edge/face is not contained in one of the children, s0 or s1 is -1. - int s0 = - boundEl.el->getSubObjOfChild(0, boundEl.subObj, boundEl.ithObj, boundEl.elType); - int s1 = - boundEl.el->getSubObjOfChild(1, boundEl.subObj, boundEl.ithObj, boundEl.elType); + int s0 = boundEl.el->getSubObjOfChild(0, boundEl.subObj, + boundEl.ithObj, boundEl.elType); + int s1 = boundEl.el->getSubObjOfChild(1, boundEl.subObj, + boundEl.ithObj, boundEl.elType); TEST_EXIT_DBG(s0 != -1 || s1 != -1)("This should not happen!\n"); diff --git a/AMDiS/src/parallel/MeshManipulation.h b/AMDiS/src/parallel/MeshManipulation.h index b74c7e06ebf289046bdba7b3fe8b03685eeea657..cce08a18b33df59c3f0aa9b88f4d1f6e6125b332 100644 --- a/AMDiS/src/parallel/MeshManipulation.h +++ b/AMDiS/src/parallel/MeshManipulation.h @@ -33,11 +33,13 @@ namespace AMDiS { class MeshManipulation { public: - MeshManipulation(FiniteElemSpace *space); + MeshManipulation(Mesh *m); ~MeshManipulation(); - void deleteDoubleDofs(std::set& newMacroEl, ElementObjects &elObj); + void deleteDoubleDofs(const FiniteElemSpace *feSpace, + std::set& newMacroEl, + ElementObjects &elObj); /** \brief * Starts the procedure to fit a given edge/face of an element with a mesh @@ -83,8 +85,6 @@ namespace AMDiS { private: - FiniteElemSpace *feSpace; - Mesh *mesh; RefinementManager *refineManager; diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc index 7dd5c52dcaaf395d64bd998299bf4abeaed8130d..7eeb7b84da1c492c87881754eec3001cc7c8841e 100644 --- a/AMDiS/src/parallel/ParallelDebug.cc +++ b/AMDiS/src/parallel/ParallelDebug.cc @@ -147,13 +147,15 @@ namespace AMDiS { return; + TEST_EXIT(pdb.feSpaces.size() == 1)("Shoudl not happen!\n"); + // === 1. check: All periodic DOFs should have at least a correct number === // === of periodic associations. === for (map >::iterator it = pdb.periodicDofAssociations.begin(); it != pdb.periodicDofAssociations.end(); ++it) { WorldVector c; - pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, c); + pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c); int nAssoc = it->second.size(); TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7)) ("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", @@ -259,11 +261,11 @@ namespace AMDiS { continue; DofContainer dofs; - boundIt->rankObj.el->getAllDofs(pdb.feSpace, boundIt->rankObj, dofs); + boundIt->rankObj.el->getAllDofs(pdb.feSpaces[0], boundIt->rankObj, dofs); for (unsigned int i = 0; i < dofs.size(); i++) { WorldVector c; - pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpace, c); + pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpaces[0], c); sendCoords[it->first].push_back(c); rankToDofType[it->first].push_back(boundIt->type); } @@ -320,6 +322,8 @@ namespace AMDiS { FUNCNAME("ParallelDebug::testCommonDofs()"); clock_t first = clock(); + // Get FE space with basis functions of the highest degree + const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1]; int testCommonDofs = 1; Parameters::get("dbg->test common dofs", testCommonDofs); @@ -342,19 +346,18 @@ namespace AMDiS { // this rank expectes to receive from. RankToCoords recvCoords; - DOFVector > coords(pdb.feSpace, "dofCorrds"); - pdb.mesh->getDofIndexCoords(pdb.feSpace, coords); + DOFVector > coords(feSpace, "dofCorrds"); + pdb.mesh->getDofIndexCoords(feSpace, coords); - for (RankToDofContainer::iterator it = pdb.sendDofs.begin(); - it != pdb.sendDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + typedef map >::iterator it_type; + for (it_type it = pdb.sendDofs.begin(); it != pdb.sendDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) sendCoords[it->first].push_back(coords[**dofIt]); - for (RankToDofContainer::iterator it = pdb.recvDofs.begin(); - it != pdb.recvDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) recvCoords[it->first].push_back(coords[**dofIt]); vector sendSize(pdb.mpiSize, 0); @@ -443,7 +446,7 @@ namespace AMDiS { oss << ")"; MSG("%s\n", oss.str().c_str()); - debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i])); + debug::printInfoByDof(feSpace, *(pdb.recvDofs[it->first][feSpace][i])); } ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank); foundError = 1; @@ -461,29 +464,31 @@ namespace AMDiS { { FUNCNAME("ParallelDebug::testGlobalIndexByCoords()"); - DOFVector > coords(pdb.feSpace, "tmp"); - pdb.mesh->getDofIndexCoords(pdb.feSpace, coords); + // Get FE space with basis functions of the highest degree + const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1]; + + DOFVector > coords(feSpace, "tmp"); + pdb.mesh->getDofIndexCoords(feSpace, coords); typedef map, int> CoordsIndexMap; CoordsIndexMap coordsToIndex; DOFIterator > it(&coords, USED_DOFS); for (it.reset(); !it.end(); ++it) - coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()]; + coordsToIndex[(*it)] = + pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()]; StdMpi stdMpi(pdb.mpiComm, true); - for (RankToDofContainer::iterator it = pdb.sendDofs.begin(); - it != pdb.sendDofs.end(); ++it) + typedef map >::iterator it_type; + for (it_type it = pdb.sendDofs.begin(); it != pdb.sendDofs.end(); ++it) stdMpi.send(it->first, coordsToIndex); - for (RankToDofContainer::iterator it = pdb.recvDofs.begin(); - it != pdb.recvDofs.end(); ++it) + for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it) stdMpi.recv(it->first); stdMpi.startCommunication(); int foundError = 0; - for (RankToDofContainer::iterator it = pdb.recvDofs.begin(); - it != pdb.recvDofs.end(); ++it) { + for (it_type it = pdb.recvDofs.begin(); it != pdb.recvDofs.end(); ++it) { CoordsIndexMap& otherCoords = stdMpi.getRecvData(it->first); for (CoordsIndexMap::iterator coordsIt = otherCoords.begin(); @@ -553,27 +558,35 @@ namespace AMDiS { void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, - RankToDofContainer &sendDofs, - RankToDofContainer &recvDofs) + map > &sendDofs, + map > &recvDofs) { FUNCNAME("ParallelDebug::testDofContainerCommunication()"); + typedef map >::iterator it_type; + map sendNumber; - for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) - sendNumber[it->first] = it->second.size(); + for (it_type it = sendDofs.begin(); it != sendDofs.end(); ++it) + for (map::iterator dcIt = it->second.begin(); dcIt != it->second.end(); ++dcIt) + sendNumber[it->first] += dcIt->second.size(); + + map recvNumber; + for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it) + for (map::iterator dcIt = it->second.begin(); dcIt != it->second.end(); ++dcIt) + recvNumber[it->first] += dcIt->second.size(); StdMpi stdMpi(pdb.mpiComm); stdMpi.send(sendNumber); - for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) + for (it_type it = recvDofs.begin(); it != recvDofs.end(); ++it) stdMpi.recv(it->first); stdMpi.startCommunication(); int foundError = 0; for (map::iterator it = stdMpi.getRecvData().begin(); it != stdMpi.getRecvData().end(); ++it) { - if (it->second != static_cast(recvDofs[it->first].size())) { + if (it->second != recvNumber[it->first]) { ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", - recvDofs[it->first].size(), it->first, it->second); + recvNumber[it->first], it->first, it->second); foundError = 1; } } @@ -618,33 +631,34 @@ namespace AMDiS { void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank) { if (rank == -1 || pdb.mpiRank == rank) { + const FiniteElemSpace *feSpace = pdb.feSpaces[0]; + cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl; - for (DofMapping::iterator it = pdb.mapLocalGlobalDofs.begin(); - it != pdb.mapLocalGlobalDofs.end(); it++) { + for (DofMapping::iterator it = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin(); + it != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); it++) { DegreeOfFreedom localdof = -1; - if (pdb.mapLocalDofIndex.count(it->first) > 0) - localdof = pdb.mapLocalDofIndex[it->first]; + if (pdb.dofFeData[feSpace].mapLocalDofIndex.count(it->first) > 0) + localdof = pdb.dofFeData[feSpace].mapLocalDofIndex[it->first]; cout << "DOF " << it->first << " " << it->second << " " << localdof << endl; WorldVector coords; - pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, coords); + pdb.mesh->getDofIndexCoords(it->first, feSpace, coords); coords.print(); - for (RankToDofContainer::iterator rankit = pdb.sendDofs.begin(); - rankit != pdb.sendDofs.end(); ++rankit) { - for (DofContainer::iterator dofit = rankit->second.begin(); - dofit != rankit->second.end(); ++dofit) + typedef map >::iterator it_type; + for (it_type rankit = pdb.sendDofs.begin(); rankit != pdb.sendDofs.end(); ++rankit) { + for (DofContainer::iterator dofit = rankit->second[feSpace].begin(); + dofit != rankit->second[feSpace].end(); ++dofit) if (**dofit == it->first) cout << "SEND DOF TO " << rankit->first << endl; } - for (RankToDofContainer::iterator rankit = pdb.recvDofs.begin(); - rankit != pdb.recvDofs.end(); ++rankit) { - for (DofContainer::iterator dofit = rankit->second.begin(); - dofit != rankit->second.end(); ++dofit) + for (it_type rankit = pdb.recvDofs.begin(); rankit != pdb.recvDofs.end(); ++rankit) { + for (DofContainer::iterator dofit = rankit->second[feSpace].begin(); + dofit != rankit->second[feSpace].end(); ++dofit) if (**dofit == it->first) cout << "RECV DOF FROM " << rankit->first << endl; } @@ -662,6 +676,8 @@ namespace AMDiS { ERROR_EXIT("Function must be rewritten!\n"); #if 0 + const FiniteElemSpace* feSpace = pdb.feSpaces[0]; + typedef map DofMapping; typedef map > PeriodicDofMap; @@ -677,15 +693,15 @@ namespace AMDiS { cout << endl; DegreeOfFreedom localdof = -1; - for (DofMapping::iterator dofIt = pdb.mapLocalGlobalDofs.begin(); - dofIt != pdb.mapLocalGlobalDofs.end(); ++dofIt) + for (DofMapping::iterator dofIt = pdb.dofFeData[feSpace].mapLocalGlobalDofs.begin(); + dofIt != pdb.dofFeData[feSpace].mapLocalGlobalDofs.end(); ++dofIt) if (dofIt->second == it->first) localdof = dofIt->first; TEST_EXIT(localdof != -1)("There is something wrong!\n"); WorldVector coords; - pdb.mesh->getDofIndexCoords(localdof, pdb.feSpace, coords); + pdb.mesh->getDofIndexCoords(localdof, feSpace, coords); coords.print(); } } @@ -706,7 +722,7 @@ namespace AMDiS { dofit != rankDofs.end(); ++dofit) { cout << " " << **dofit << endl; WorldVector coords; - pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords); + pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords); coords.print(); } @@ -715,7 +731,7 @@ namespace AMDiS { dofit != rankAllDofs.end(); ++dofit) { cout << " " << **dofit << endl; WorldVector coords; - pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords); + pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords); coords.print(); } } @@ -763,21 +779,23 @@ namespace AMDiS { { FUNCNAME("ParallelDebug::writeCoordsFile()"); + const FiniteElemSpace *feSpace = pdb.feSpaces[0]; + stringstream filename; filename << prefix << "-" << pdb.mpiRank << "." << postfix; - DOFVector > coords(pdb.feSpace, "tmp"); - pdb.mesh->getDofIndexCoords(pdb.feSpace, coords); + DOFVector > coords(feSpace, "tmp"); + pdb.mesh->getDofIndexCoords(feSpace, coords); typedef map > ElDofMap; ElDofMap elDofMap; TraverseStack stack; - const BasisFunction *basisFcts = pdb.feSpace->getBasisFcts(); + const BasisFunction *basisFcts = feSpace->getBasisFcts(); vector localIndices(basisFcts->getNumber()); ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { basisFcts->getLocalIndices(elInfo->getElement(), - pdb.feSpace->getAdmin(), localIndices); + feSpace->getAdmin(), localIndices); elDofMap[elInfo->getElement()->getIndex()] = localIndices; elInfo = stack.traverseNext(elInfo); } @@ -792,8 +810,8 @@ namespace AMDiS { DOFIterator > it(&coords, USED_DOFS); for (it.reset(); !it.end(); ++it) { file << it.getDOFIndex() << " " - << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " " - << pdb.getIsRankDof(it.getDOFIndex()); + << pdb.dofFeData[feSpace].mapLocalGlobalDofs[it.getDOFIndex()] << " " + << pdb.getIsRankDof(feSpace, it.getDOFIndex()); for (int i = 0; i < pdb.mesh->getDim(); i++) file << " " << (*it)[i]; file << "\n"; @@ -838,7 +856,7 @@ namespace AMDiS { void ParallelDebug::writePartitioningFile(string filename, int counter, - FiniteElemSpace *feSpace) + const FiniteElemSpace *feSpace) { FUNCNAME("ParallelDebug::writePartitioningFile()"); @@ -878,13 +896,13 @@ namespace AMDiS { if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex)) debug::writeLocalElementDofs(pdb.mpiRank, it->rankObj.elIndex, - pdb.feSpace); + pdb.feSpaces[0]); for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it) if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex)) debug::writeLocalElementDofs(pdb.mpiRank, it->rankObj.elIndex, - pdb.feSpace); + pdb.feSpaces[0]); } diff --git a/AMDiS/src/parallel/ParallelDebug.h b/AMDiS/src/parallel/ParallelDebug.h index 8b1a9aab38ac28f4270a39d37953ec2573abca69..e6069cc79ab3ee150add4026efbb70f5cb30026b 100644 --- a/AMDiS/src/parallel/ParallelDebug.h +++ b/AMDiS/src/parallel/ParallelDebug.h @@ -90,8 +90,8 @@ namespace AMDiS { * \param[in] recvDofs The map of all DOFs the rank will receive. */ static void testDofContainerCommunication(MeshDistributor &pdb, - RankToDofContainer &sendDofs, - RankToDofContainer &recvDofs); + map > &sendDofs, + map > &recvDofs); /// Tests if there are multiple DOFs in mesh with the same coords. static void testDoubleDofs(Mesh *mesh); @@ -159,7 +159,7 @@ namespace AMDiS { */ static void writePartitioningFile(std::string filename, int counter, - FiniteElemSpace *feSpace); + const FiniteElemSpace *feSpace); static bool followThisBound(int rankElIndex, int neighElIndex); diff --git a/AMDiS/src/parallel/ParallelMapper.h b/AMDiS/src/parallel/ParallelMapper.h index 5d17e11434ed26fb5cc917e6d790a92c710aa067..a9fa1ec7e89d5db4eaf4784ad492adb7aa70d037 100644 --- a/AMDiS/src/parallel/ParallelMapper.h +++ b/AMDiS/src/parallel/ParallelMapper.h @@ -3,45 +3,51 @@ #include "parallel/MeshDistributor.h" namespace AMDiS { - //todo: request from a problem -> distributor mapper? - template< typename size_type > - class ParallelMapper_base { - MeshDistributor& md; - /// the current row in the problem system - size_type r; - /// the current column in the problem system - size_type c; - int nComponents; + //todo: request from a problem -> distributor mapper? + template< typename size_type > + class ParallelMapper_base { + MeshDistributor& md; + /// the current row in the problem system + size_type r; + /// the current column in the problem system + size_type c; + int nComponents; + + public: + + ParallelMapper_base(MeshDistributor& md, int nComponents): + md(md),r(0),c(0),nComponents(nComponents) + { + TEST_EXIT_DBG(nComponents > 0)("the system must have at least one component"); + } + size_type row(size_type r_) + { + // TODO: JUST TO COMPILE, mapLocalToGlobal is now FE space dependent! + size_type ret = 0; + // size_type ret = md.mapLocalToGlobal(r_)*nComponents+r; + return ret; + } + size_type col(size_type c_) + { + // TODO: JUST TO COMPILE, mapLocalToGlobal is now FE space dependent! + size_type ret = 0; + // size_type ret = md.mapLocalToGlobal(c_)*nComponents+c; + return ret ; + } + + inline MeshDistributor& meshDistributor() { return md; } + inline int getNumComponents() { return nComponents; } + + inline void setRow(size_type r_) { r = r_ ;} + inline void setCol(size_type c_) { c = c_; } + // TODO: JUST TO COMPILE, getNumberOverallDofs() is now FE space dependent! + inline size_type getNumRows() { return nComponents; } + inline size_type getNumCols() { return nComponents; } - public: - - ParallelMapper_base(MeshDistributor& md, int nComponents): - md(md),r(0),c(0),nComponents(nComponents) - { - TEST_EXIT_DBG(nComponents > 0)("the system must have at least one component"); - } - size_type row(size_type r_) - { - size_type ret = md.mapLocalToGlobal(r_)*nComponents+r; - //size_type ret = md.mapLocalToGlobal(r_)+r*md.getNumberOverallDofs() ; - return ret; - } - size_type col(size_type c_) - { - size_type ret = md.mapLocalToGlobal(c_)*nComponents+c; - //size_type ret = md.mapLocalToGlobal(c_)+c*md.getNumberOverallDofs() ; - return ret ; - } - - inline MeshDistributor& meshDistributor() { return md; } - inline int getNumComponents() { return nComponents; } - - inline void setRow(size_type r_) { r = r_ ;} - inline void setCol(size_type c_) { c = c_; } - inline size_type getNumRows() { return md.getNumberOverallDofs() * nComponents; } - inline size_type getNumCols() { return md.getNumberOverallDofs() * nComponents; } - }; - - typedef ParallelMapper_base< unsigned int > ParallelMapper; +/* inline size_type getNumRows() { return md.getNumberOverallDofs() * nComponents; } */ +/* inline size_type getNumCols() { return md.getNumberOverallDofs() * nComponents; } */ + }; + + typedef ParallelMapper_base< unsigned int > ParallelMapper; } #endif diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index dcbc0b4d1e66dbe7733759a9f42e187e68087dbf..1a58c810da13614a1af2f5308e1dc1de5ce47c8a 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -229,6 +229,8 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createPrimals()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + // === Define all vertices on the interior boundaries of the macro mesh === // === to be primal variables. === @@ -236,7 +238,7 @@ namespace AMDiS { DofIndexSet primals; DofContainerSet& vertices = - meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX]; + meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX]; TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n"); for (DofContainerSet::iterator it = vertices.begin(); it != vertices.end(); ++it) @@ -249,7 +251,7 @@ namespace AMDiS { globalPrimalIndex.clear(); nRankPrimals = 0; for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) - if (meshDistributor->getIsRankDof(*it)) { + if (meshDistributor->getIsRankDof(feSpace, *it)) { globalPrimalIndex[*it] = nRankPrimals; nRankPrimals++; } @@ -277,25 +279,25 @@ namespace AMDiS { // === Communicate primal's global index from ranks that own the === // === primals to ranks that contain this primals but are not owning === // === them. === + + typedef map >::iterator it_type; StdMpi > stdMpi(meshDistributor->getMpiComm()); - RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); - for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (it_type it = meshDistributor->getSendDofs().begin(); + it != meshDistributor->getSendDofs().end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (globalPrimalIndex.count(**dofIt)) stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]); stdMpi.updateSendDataSize(); - RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { bool recvFromRank = false; - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (primals.count(**dofIt) && - meshDistributor->getIsRankDof(**dofIt) == false) { + meshDistributor->getIsRankDof(feSpace, **dofIt) == false) { recvFromRank = true; break; } @@ -305,13 +307,13 @@ namespace AMDiS { } stdMpi.startCommunication(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { int i = 0; - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) { + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) { if (primals.count(**dofIt) && - meshDistributor->getIsRankDof(**dofIt) == false) + meshDistributor->getIsRankDof(feSpace, **dofIt) == false) globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++]; } } @@ -325,17 +327,20 @@ namespace AMDiS { void PetscSolverFeti::createDuals() { FUNCNAME("PetscSolverFeti::createDuals()"); + + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); // === Create for each dual node that is owned by the rank, the set === // === of ranks that contain this node (denoted by W(x_j)). === boundaryDofRanks.clear(); - RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); - for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) { - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) { + typedef map >::iterator it_type; + + for (it_type it = meshDistributor->getSendDofs().begin(); + it != meshDistributor->getSendDofs().end(); ++it) { + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) { // If DOF is not primal, i.e., its a dual node if (globalPrimalIndex.count(**dofIt) == 0) { boundaryDofRanks[**dofIt].insert(mpiRank); @@ -349,21 +354,20 @@ namespace AMDiS { // === ranks that also have this node. === StdMpi > > stdMpi(meshDistributor->getMpiComm()); - for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (it_type it = meshDistributor->getSendDofs().begin(); + it != meshDistributor->getSendDofs().end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (globalPrimalIndex.count(**dofIt) == 0) stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]); stdMpi.updateSendDataSize(); - RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { bool recvFromRank = false; - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (globalPrimalIndex.count(**dofIt) == 0) { recvFromRank = true; break; @@ -374,13 +378,13 @@ namespace AMDiS { } stdMpi.startCommunication(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { int i = 0; - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (globalPrimalIndex.count(**dofIt) == 0) - boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++]; + boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++]; } @@ -389,14 +393,14 @@ namespace AMDiS { duals.clear(); globalDualIndex.clear(); - int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs(); + int nRankAllDofs = feSpace->getAdmin()->getUsedDofs(); nRankB = nRankAllDofs - globalPrimalIndex.size(); nOverallB = 0; rStartB = 0; mpi::getDofNumbering(meshDistributor->getMpiComm(), nRankB, rStartB, nOverallB); DofContainer allBoundaryDofs; - meshDistributor->getAllBoundaryDofs(allBoundaryDofs); + meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs); int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size(); int nRankDuals = 0; @@ -421,12 +425,14 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createLagrange()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + // === Reserve for each dual node, on the rank that owns this node, the === // === appropriate number of Lagrange constraints. === nRankLagrange = 0; for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) { - if (meshDistributor->getIsRankDof(*it)) { + if (meshDistributor->getIsRankDof(feSpace, *it)) { dofFirstLagrange[*it] = nRankLagrange; int degree = boundaryDofRanks[*it].size(); nRankLagrange += (degree * (degree - 1)) / 2; @@ -444,7 +450,7 @@ namespace AMDiS { nRankLagrange, rStartLagrange, nOverallLagrange); for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) - if (meshDistributor->getIsRankDof(*it)) + if (meshDistributor->getIsRankDof(feSpace, *it)) dofFirstLagrange[*it] += rStartLagrange; MSG("nRankLagrange = %d nOverallLagrange = %d\n", @@ -454,11 +460,12 @@ namespace AMDiS { // === Communicate dofFirstLagrange to all other ranks. === StdMpi > stdMpi(meshDistributor->getMpiComm()); - RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); - for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) { + typedef map >::iterator it_type; + + for (it_type it = meshDistributor->getSendDofs().begin(); + it != meshDistributor->getSendDofs().end(); ++it) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) { if (globalPrimalIndex.count(**dofIt) == 0) { TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n"); stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]); @@ -466,12 +473,11 @@ namespace AMDiS { } stdMpi.updateSendDataSize(); - RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { bool recvData = false; - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) if (globalPrimalIndex.count(**dofIt) == 0) { recvData = true; break; @@ -483,12 +489,13 @@ namespace AMDiS { stdMpi.startCommunication(); - for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { + for (it_type it = meshDistributor->getRecvDofs().begin(); + it != meshDistributor->getRecvDofs().end(); ++it) { int counter = 0; - for (unsigned int i = 0; i < it->second.size(); i++) - if (globalPrimalIndex.count(*(it->second[i])) == 0) - dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++]; + DofContainer &dc = it->second[feSpace]; + for (unsigned int i = 0; i < dc.size(); i++) + if (globalPrimalIndex.count(*(dc[i])) == 0) + dofFirstLagrange[*(dc[i])] = stdMpi.getRecvData(it->first)[counter++]; } } @@ -497,8 +504,9 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createIndeB()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); localIndexB.clear(); - DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin(); + DOFAdmin* admin = feSpace->getAdmin(); // === To ensure that all interior node on each rank are listen first in === // === the global index of all B nodes, insert all interior nodes first, === diff --git a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc index 6d9cf30807f19ac265fb277fb2bca05b47206a22..94a0546d1b4162bbd099f781072067ab469a1944 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc @@ -24,9 +24,10 @@ namespace AMDiS { TEST_EXIT_DBG(mat)("No DOF matrix defined!\n"); double wtime = MPI::Wtime(); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); nComponents = mat->getNumRows(); - int nRankRows = meshDistributor->getNumberRankDofs(); - int nOverallRows = meshDistributor->getNumberOverallDofs(); + int nRankRows = meshDistributor->getNumberRankDofs(feSpace); + int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace); #if (DEBUG != 0) MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime); @@ -81,8 +82,9 @@ namespace AMDiS { TEST_EXIT_DBG(vec)("NO DOF vector defined!\n"); nComponents = vec->getSize(); - int nRankRows = meshDistributor->getNumberRankDofs(); - int nOverallRows = meshDistributor->getNumberOverallDofs(); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + int nRankRows = meshDistributor->getNumberRankDofs(feSpace); + int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace); nestVec.resize(nComponents); @@ -108,6 +110,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalBlockMatrix::solvePetscMatrix()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); VecDuplicate(petscRhsVec, &petscSolVec); // PETSc. @@ -120,12 +123,12 @@ namespace AMDiS { Vec tmp; VecNestGetSubVec(petscSolVec, i, &tmp); - int nRankDofs = meshDistributor->getNumberRankDofs(); + int nRankDofs = meshDistributor->getNumberRankDofs(feSpace); PetscScalar *vecPointer; VecGetArray(tmp, &vecPointer); for (int j = 0; j < nRankDofs; j++) - dofvec[meshDistributor->mapLocalToDofIndex(j)] = vecPointer[j]; + dofvec[meshDistributor->mapLocalToDofIndex(feSpace, j)] = vecPointer[j]; VecRestoreArray(tmp, &vecPointer); } @@ -165,6 +168,8 @@ namespace AMDiS { TEST_EXIT(mat)("No DOFMatrix!\n"); TEST_EXIT(petscMat)("No PETSc matrix!\n"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; @@ -187,7 +192,7 @@ namespace AMDiS { cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { // Global index of the current row DOF. - int rowIndex = meshDistributor->mapLocalToGlobal(*cursor); + int rowIndex = meshDistributor->mapLocalToGlobal(feSpace, *cursor); cols.clear(); values.clear(); @@ -195,7 +200,7 @@ namespace AMDiS { for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Global index of the current column index. - int colIndex = meshDistributor->mapLocalToGlobal(col(*icursor)); + int colIndex = meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)); // Ignore all zero entries, expect it is a diagonal entry. if (value(*icursor) == 0.0 && rowIndex != colIndex) @@ -217,10 +222,12 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalBlockMatrix::setDofVector()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + // Traverse all used DOFs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { - int index = meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex()); + int index = meshDistributor->mapLocalToGlobal(feSpace, dofIt.getDOFIndex()); double value = *dofIt; VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index 874fe7dd4589a088ee6c0d39da7274b01ae88c06..b9c6e831f590b67e0b1256ad8cf78edbc0acc8de 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -20,13 +20,16 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n"); TEST_EXIT_DBG(mat)("No DOF matrix defined!\n"); double wtime = MPI::Wtime(); int nComponents = mat->getNumRows(); - int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; - int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents; + int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents; + int nOverallRows = + meshDistributor->getNumberOverallDofs(feSpace) * nComponents; // === Create PETSc vector (solution and a temporary vector). === @@ -66,9 +69,9 @@ namespace AMDiS { #if (DEBUG != 0) int a, b; MatGetOwnershipRange(petscMatrix, &a, &b); - TEST_EXIT(a == meshDistributor->getStartDofs() * nComponents) + TEST_EXIT(a == meshDistributor->getStartDofs(feSpace) * nComponents) ("Wrong matrix ownership range!\n"); - TEST_EXIT(b == meshDistributor->getStartDofs() * nComponents + nRankRows) + TEST_EXIT(b == meshDistributor->getStartDofs(feSpace) * nComponents + nRankRows) ("Wrong matrix ownership range!\n"); #endif @@ -108,12 +111,14 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + TEST_EXIT_DBG(vec)("No DOF vector defined!\n"); TEST_EXIT_DBG(meshDistributor)("No mesh distributor defined!\n"); int nComponents = vec->getSize(); - int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; - int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents; + int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents; + int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace) * nComponents; VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscRhsVec); @@ -132,6 +137,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); int nComponents = vec.getSize(); // === Set old solution to be initiual guess for PETSc solver. === @@ -149,14 +155,14 @@ namespace AMDiS { KSPSolve(solver, petscRhsVec, petscSolVec); // === Transfere values from PETSc's solution vectors to the DOF vectors. === - int nRankDofs = meshDistributor->getNumberRankDofs(); + int nRankDofs = meshDistributor->getNumberRankDofs(feSpace); PetscScalar *vecPointer; VecGetArray(petscSolVec, &vecPointer); for (int i = 0; i < nComponents; i++) { DOFVector &dofvec = *(vec.getDOFVector(i)); for (int j = 0; j < nRankDofs; j++) - dofvec[meshDistributor->mapLocalToDofIndex(j)] = + dofvec[meshDistributor->mapLocalToDofIndex(feSpace, j)] = vecPointer[j * nComponents + i]; } @@ -194,6 +200,8 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + TEST_EXIT(mat)("No DOFMatrix!\n"); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; @@ -221,7 +229,7 @@ namespace AMDiS { cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { // Global index of the current row DOF. - int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor); + int globalRowDof = meshDistributor->mapLocalToGlobal(feSpace, *cursor); // Test if the current row DOF is a periodic DOF. bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof); @@ -239,7 +247,8 @@ namespace AMDiS { icursor != icend; ++icursor) { // Global index of the current column index. - int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor)); + int globalColDof = + meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)); // Test if the current col dof is a periodic dof. bool periodicCol = meshDistributor->isPeriodicDof(globalColDof); // Calculate PETSc col index. @@ -316,7 +325,7 @@ namespace AMDiS { icursor != icend; ++icursor) { // Global index of the current column index. - int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor)); + int globalColDof = meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)); // Ignore all zero entries, expect it is a diagonal entry. if (value(*icursor) == 0.0 && globalRowDof != globalColDof) @@ -418,15 +427,17 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::setDofVector()"); + const FiniteElemSpace *feSpace; + // Traverse all used DOFs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { - if (rankOnly && !meshDistributor->getIsRankDof(dofIt.getDOFIndex())) + if (rankOnly && !meshDistributor->getIsRankDof(feSpace, dofIt.getDOFIndex())) continue; // Calculate global row index of the DOF. DegreeOfFreedom globalRowDof = - meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex()); + meshDistributor->mapLocalToGlobal(feSpace, dofIt.getDOFIndex()); // Calculate PETSc index of the row DOF. int index = globalRowDof * dispMult + dispAdd; @@ -456,8 +467,9 @@ namespace AMDiS { TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n"); TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n"); + const FiniteElemSpace* feSpace = meshDistributor->getFeSpace(0); int nComponents = mat->getNumRows(); - int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; + int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents; d_nnz = new int[nRankRows]; o_nnz = new int[nRankRows]; for (int i = 0; i < nRankRows; i++) { @@ -483,18 +495,19 @@ namespace AMDiS { // First, create for all ranks we send data to MatrixNnzEntry object with 0 entries. - for (RankToDofContainer::iterator it = meshDistributor->getRecvDofs().begin(); + typedef map >::iterator it_type; + for (it_type it = meshDistributor->getRecvDofs().begin(); it != meshDistributor->getRecvDofs().end(); ++it) { sendMatrixEntry[it->first].resize(0); - for (DofContainer::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) + for (DofContainer::iterator dofIt = it->second[feSpace].begin(); + dofIt != it->second[feSpace].end(); ++dofIt) sendDofToRank[**dofIt] = it->first; } std::set recvFromRank; - for (RankToDofContainer::iterator it = meshDistributor->getSendDofs().begin(); + for (it_type it = meshDistributor->getSendDofs().begin(); it != meshDistributor->getSendDofs().end(); ++it) recvFromRank.insert(it->first); @@ -515,12 +528,12 @@ namespace AMDiS { for (cursor_type cursor = begin(bmat), cend = end(bmat); cursor != cend; ++cursor) { - int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor); + int globalRowDof = meshDistributor->mapLocalToGlobal(feSpace, *cursor); vector rows; rows.push_back(globalRowDof); vector rowRank; - if (meshDistributor->getIsRankDof(*cursor)) { + if (meshDistributor->getIsRankDof(feSpace, *cursor)) { rowRank.push_back(meshDistributor->getMpiRank()); } else { // Find out who is the member of this DOF. @@ -534,22 +547,22 @@ namespace AMDiS { int petscRowIdx = globalRowDof * nComponents + i; - if (meshDistributor->getIsRankDof(*cursor)) { + if (meshDistributor->getIsRankDof(feSpace, *cursor)) { // === The current row DOF is a rank dof, so create the corresponding === // === nnz values directly on rank's nnz data. === // This is the local row index of the local PETSc matrix. int localPetscRowIdx = - petscRowIdx - meshDistributor->getStartDofs() * nComponents; + petscRowIdx - meshDistributor->getStartDofs(feSpace) * nComponents; TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows) ("Should not happen! \n Debug info: localRowIdx = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d nCompontens = %d nRankRows = %d\n", *cursor, - meshDistributor->mapLocalToGlobal(*cursor), + meshDistributor->mapLocalToGlobal(feSpace, *cursor), petscRowIdx, localPetscRowIdx, - meshDistributor->getStartDofs(), + meshDistributor->getStartDofs(feSpace), nComponents, nRankRows); @@ -558,13 +571,13 @@ namespace AMDiS { for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { int petscColIdx = - meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j; + meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)) * nComponents + j; if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) { // The row DOF is a rank DOF, if also the column is a rank DOF, // increment the d_nnz values for this row, otherwise the o_nnz value. - if (petscColIdx >= meshDistributor->getStartDofs() * nComponents && - petscColIdx < meshDistributor->getStartDofs() * nComponents + nRankRows) + if (petscColIdx >= meshDistributor->getStartDofs(feSpace) * nComponents && + petscColIdx < meshDistributor->getStartDofs(feSpace) * nComponents + nRankRows) d_nnz[localPetscRowIdx]++; else o_nnz[localPetscRowIdx]++; @@ -583,7 +596,7 @@ namespace AMDiS { icend = end(cursor); icursor != icend; ++icursor) { if (value(*icursor) != 0.0) { int petscColIdx = - meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j; + meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)) * nComponents + j; sendMatrixEntry[sendToRank]. push_back(make_pair(petscRowIdx, petscColIdx)); @@ -615,14 +628,14 @@ namespace AMDiS { int r = it->second[i].first; int c = it->second[i].second; - int localRowIdx = r - meshDistributor->getStartDofs() * nComponents; + int localRowIdx = r - meshDistributor->getStartDofs(feSpace) * nComponents; TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows) ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n", r, localRowIdx, nRankRows, it->first); - if (c < meshDistributor->getStartDofs() * nComponents || - c >= meshDistributor->getStartDofs() * nComponents + nRankRows) + if (c < meshDistributor->getStartDofs(feSpace) * nComponents || + c >= meshDistributor->getStartDofs(feSpace) * nComponents + nRankRows) o_nnz[localRowIdx]++; else d_nnz[localRowIdx]++; diff --git a/AMDiS/src/parallel/PetscSolverSchur.cc b/AMDiS/src/parallel/PetscSolverSchur.cc index 2409c18ac8927e9c994d30c8058d4fcc2558f17d..c0dd58d34ff558e29ecdfd4b7c184309b606ffa9 100644 --- a/AMDiS/src/parallel/PetscSolverSchur.cc +++ b/AMDiS/src/parallel/PetscSolverSchur.cc @@ -25,20 +25,21 @@ namespace AMDiS { TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n"); MPI::Intracomm& mpiComm = meshDistributor->getMpiComm(); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); typedef map RankToDofContainer; typedef map DofIndexToBool; boundaryDofs.clear(); std::set boundaryLocalDofs; - RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); - for (RankToDofContainer::iterator rankIt = sendDofs.begin(); - rankIt != sendDofs.end(); ++rankIt) { - for (DofContainer::iterator dofIt = rankIt->second.begin(); - dofIt != rankIt->second.end(); ++dofIt) { + typedef map >::iterator it_type; + for (it_type rankIt = meshDistributor->getSendDofs().begin(); + rankIt != meshDistributor->getSendDofs().end(); ++rankIt) { + for (DofContainer::iterator dofIt = rankIt->second[feSpace].begin(); + dofIt != rankIt->second[feSpace].end(); ++dofIt) { boundaryLocalDofs.insert(**dofIt); - boundaryDofs.insert(meshDistributor->mapLocalToGlobal(**dofIt)); + boundaryDofs.insert(meshDistributor->mapLocalToGlobal(feSpace, **dofIt)); } } @@ -49,8 +50,10 @@ namespace AMDiS { rStartBoundaryDofs, nOverallBoundaryDofs); - DofContainerSet& edgeDofs = meshDistributor->getBoundaryDofInfo().geoDofs[EDGE]; - DofContainerSet& vertexDofs = meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX]; + DofContainerSet& edgeDofs = + meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[EDGE]; + DofContainerSet& vertexDofs = + meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX]; int nEdgeDofs = edgeDofs.size(); int nVertexDofs = vertexDofs.size(); @@ -75,14 +78,14 @@ namespace AMDiS { int counter = rStartEdgeDofs; for (DofContainerSet::iterator it = edgeDofs.begin(); it != edgeDofs.end(); ++it) - mapGlobalBoundaryDof[meshDistributor->mapLocalToGlobal(**it)] = + mapGlobalBoundaryDof[meshDistributor->mapLocalToGlobal(feSpace, **it)] = counter++; } { int counter = nOverallEdgeDofs + rStartVertexDofs; for (DofContainerSet::iterator it = vertexDofs.begin(); it != vertexDofs.end(); ++it) - mapGlobalBoundaryDof[meshDistributor->mapLocalToGlobal(**it)] = + mapGlobalBoundaryDof[meshDistributor->mapLocalToGlobal(feSpace, **it)] = counter++; } #else @@ -97,21 +100,20 @@ namespace AMDiS { std::set otherBoundaryLocalDofs; - RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); - for (RankToDofContainer::iterator rankIt = recvDofs.begin(); - rankIt != recvDofs.end(); ++rankIt) - for (DofContainer::iterator dofIt = rankIt->second.begin(); - dofIt != rankIt->second.end(); ++dofIt) + for (it_type rankIt = meshDistributor->getRecvDofs().begin(); + rankIt != meshDistributor->getRecvDofs().end(); ++rankIt) + for (DofContainer::iterator dofIt = rankIt->second[feSpace].begin(); + dofIt != rankIt->second[feSpace].end(); ++dofIt) otherBoundaryLocalDofs.insert(**dofIt); interiorDofs.clear(); - DofIndexToBool& isRankDof = meshDistributor->getIsRankDof(); + DofIndexToBool& isRankDof = meshDistributor->getIsRankDof(feSpace); for (DofIndexToBool::iterator dofIt = isRankDof.begin(); dofIt != isRankDof.end(); ++dofIt) { if (dofIt->second && boundaryLocalDofs.count(dofIt->first) == 0 && otherBoundaryLocalDofs.count(dofIt->first) == 0) - interiorDofs.insert(meshDistributor->mapLocalToGlobal(dofIt->first)); + interiorDofs.insert(meshDistributor->mapLocalToGlobal(feSpace, dofIt->first)); } nInteriorDofs = interiorDofs.size(); @@ -130,14 +132,14 @@ namespace AMDiS { TEST_EXIT_DBG(nInteriorDofs > 0)("Should not happen!\n"); - StdMpi > stdMpi(mpiComm, false); - for (RankToDofContainer::iterator sendIt = meshDistributor->getSendDofs().begin(); + StdMpi > stdMpi(mpiComm); + for (it_type sendIt = meshDistributor->getSendDofs().begin(); sendIt != meshDistributor->getSendDofs().end(); ++sendIt) { stdMpi.getSendData(sendIt->first).resize(0); stdMpi.getSendData(sendIt->first).reserve(sendIt->second.size()); - for (DofContainer::iterator dofIt = sendIt->second.begin(); - dofIt != sendIt->second.end(); ++dofIt) { - int globalSendDof = meshDistributor->mapLocalToGlobal(**dofIt); + for (DofContainer::iterator dofIt = sendIt->second[feSpace].begin(); + dofIt != sendIt->second[feSpace].end(); ++dofIt) { + int globalSendDof = meshDistributor->mapLocalToGlobal(feSpace, **dofIt); TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalSendDof)) ("No mapping for boundary DOF %d!\n", globalSendDof); @@ -147,16 +149,20 @@ namespace AMDiS { } stdMpi.updateSendDataSize(); - stdMpi.recv(meshDistributor->getRecvDofs()); + for (it_type recvIt = meshDistributor->getRecvDofs().begin(); + recvIt != meshDistributor->getRecvDofs().end(); ++recvIt) + stdMpi.recv(recvIt->first); + stdMpi.startCommunication(); - for (RankToDofContainer::iterator recvIt = meshDistributor->getRecvDofs().begin(); + for (it_type recvIt = meshDistributor->getRecvDofs().begin(); recvIt != meshDistributor->getRecvDofs().end(); ++recvIt) { int i = 0; - for (DofContainer::iterator dofIt = recvIt->second.begin(); - dofIt != recvIt->second.end(); ++dofIt) { - int globalRecvDof = meshDistributor->mapLocalToGlobal(**dofIt); + for (DofContainer::iterator dofIt = recvIt->second[feSpace].begin(); + dofIt != recvIt->second[feSpace].end(); ++dofIt) { + int globalRecvDof = meshDistributor->mapLocalToGlobal(feSpace, **dofIt); mapGlobalBoundaryDof[globalRecvDof] = stdMpi.getRecvData(recvIt->first)[i++]; + boundaryDofs.insert(globalRecvDof); } } @@ -179,7 +185,8 @@ namespace AMDiS { void PetscSolverSchur::fillPetscMatrix(Matrix *mat) { FUNCNAME("PetscSolverSchur::fillPetscMatrix()"); - + + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); int nComponents = mat->getNumRows(); updateDofData(nComponents); @@ -243,8 +250,8 @@ namespace AMDiS { MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY); - int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; - int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents; + int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents; + int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace) * nComponents; VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscSolVec); VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscTmpVec); @@ -255,9 +262,10 @@ namespace AMDiS { { FUNCNAME("PetscSolverSchur::fillPetscRhs()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); int nComponents = vec->getSize(); - int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; - int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents; + int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents; + int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace) * nComponents; VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscRhsVec); @@ -274,6 +282,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverSchur::solvePetscMatrix()"); + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); int nComponents = vec.getSize(); KSPCreate(PETSC_COMM_WORLD, &solver); @@ -300,7 +309,7 @@ namespace AMDiS { DOFVector::Iterator dofIt(vec.getDOFVector(i), USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { DegreeOfFreedom globalRowDof = - meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex()); + meshDistributor->mapLocalToGlobal(feSpace, dofIt.getDOFIndex()); if (boundaryDofs.count(globalRowDof)) { int index = (mapGlobalBoundaryDof[globalRowDof] - rStartBoundaryDofs + nInteriorDofs) * (i + 1); @@ -349,6 +358,7 @@ namespace AMDiS { TEST_EXIT(mat)("No DOFMatrix!\n"); + const FiniteElemSpace* feSpace = meshDistributor->getFeSpace(0); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; @@ -370,7 +380,7 @@ namespace AMDiS { cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { // Global index of the current row DOF. - int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor); + int globalRowDof = meshDistributor->mapLocalToGlobal(feSpace, *cursor); colsBoundary.clear(); colsInterior.clear(); @@ -379,7 +389,7 @@ namespace AMDiS { for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { - int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor)); + int globalColDof = meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)); if (boundaryDofs.count(globalColDof)) { TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalColDof)) @@ -434,14 +444,16 @@ namespace AMDiS { void PetscSolverSchur::setDofVector(Vec& petscVec, DOFVector* vec, int dispMult, int dispAdd, bool rankOnly) { + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); + DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { - if (rankOnly && !meshDistributor->getIsRankDof(dofIt.getDOFIndex())) + if (rankOnly && !meshDistributor->getIsRankDof(feSpace, dofIt.getDOFIndex())) continue; // Calculate global row index of the DOF. DegreeOfFreedom globalRowDof = - meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex()); + meshDistributor->mapLocalToGlobal(feSpace, dofIt.getDOFIndex()); double value = *dofIt;