diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc index a007b8b7400665e31eb3e2959a197e98aac5456e..c19e885c53625454023c1da2a8025951833f56c6 100644 --- a/AMDiS/src/DOFAdmin.cc +++ b/AMDiS/src/DOFAdmin.cc @@ -108,8 +108,8 @@ namespace AMDiS { { FUNCNAME("DOFAdmin::freeDofIndex()"); - TEST_EXIT_DBG(usedCount > 0)("no dofs in use\n"); - TEST_EXIT_DBG(dof >= 0 && dof < size)("invalid dof index %d\n",dof); + TEST_EXIT_DBG(usedCount > 0)("No DOFs in use!\n"); + TEST_EXIT_DBG(dof >= 0 && dof < size)("Invalid DOF index %d!\n", dof); std::list<DOFIndexedBase*>::iterator di; std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end(); @@ -277,8 +277,9 @@ namespace AMDiS { for (it.reset(); !it.end(); ++it) newDofIndex[it.getDOFIndex()] = 1; + // create a MONOTONE compress int n = 0, last = 0; - for (int i = 0; i < size; i++) { /* create a MONOTONE compress */ + for (int i = 0; i < size; i++) { if (newDofIndex[i] == 1) { newDofIndex[i] = n++; last = i; diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index bc9c922169ce309eec8d6fc6979e20d7184f2d14..07633a87c61a2fa880d83590c57f91245c7afeaa 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -631,12 +631,26 @@ namespace AMDiS { void Element::getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr) + bool baseDofPtr, + vector<GeoIndex>* dofGeoIndex) { getNodeDofs(feSpace, bound, dofs, baseDofPtr); + if (dofGeoIndex != NULL) { + // In the case dofGeoIndex should be filled, set all node DOFs to be + // vertex DOFs. + dofGeoIndex->resize(dofs.size()); + for (unsigned int i = 0; i < dofs.size(); i++) + (*dofGeoIndex)[i] = VERTEX; + } + if (feSpace->getBasisFcts()->getDegree() > 1) - getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr); + getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr, dofGeoIndex); + + if (dofGeoIndex) { + TEST_EXIT_DBG(dofs.size() == dofGeoIndex->size()) + ("Arrays do not fit together: %d %d\n", dofs.size(), dofGeoIndex->size()); + } } } diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index c5aff3e12d6b58b33102e979691e2671d6a589d9..18683c5a3f863388fa7df8d52b7466c086799b18 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -427,18 +427,24 @@ namespace AMDiS { * dof* [\ref dof] of the element is inserted. If * false, &(dof[.][n0]) is put to the result vector, * with n0 beging the number of predofs. + * \param[out] dofGeoIndex Optional, the function can store to each DOF in + * the DofContainer dofs the geometric index, thus + * identifing the DOF to be a vertex, edge, face or + * center DOF. */ virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr = false) const = 0; + bool baseDofPtr = false, + vector<GeoIndex>* dofGeoIndex = NULL) const = 0; /// Combines \ref getNodeDofs and \ref getHigherOrderDofs to one function. /// See parameter description there. void getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr = false); + bool baseDofPtr = false, + vector<GeoIndex>* dofGeoIndex = NULL); /** \} */ diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index 9a6a1c989230387e67b3dded9f8836e1c2ad7e38..661db9f71e7e1bf67194bd1c79270b45a9fe7341 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -181,7 +181,7 @@ namespace AMDiS { } void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, - DofContainer&, bool) const + DofContainer&, bool, vector<GeoIndex>*) const { FUNCNAME("Line::getHigherOrderDofs()"); ERROR_EXIT("Not yet implemented!\n"); diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index 4da2debf869dabe444976cbd63906490cf3bd815..669e68bccb1c822ddd6384be5cc666fd3031b72c 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -376,10 +376,13 @@ namespace AMDiS { void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr) const + bool baseDofPtr, + vector<GeoIndex>* dofGeoIndex) const { FUNCNAME("Tetrahedron::getHigherOrderDofs()"); + TEST_EXIT(dofGeoIndex != NULL)("Not yet implemented!\n"); + switch (bound.subObj) { case VERTEX: return; @@ -403,14 +406,18 @@ namespace AMDiS { if (bound.reverseMode) { if (nextBound1.ithObj >= 0) - child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, baseDofPtr); + child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, + baseDofPtr, dofGeoIndex); if (nextBound0.ithObj >= 0) - child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, baseDofPtr); + child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, + baseDofPtr, dofGeoIndex); } else { if (nextBound0.ithObj >= 0) - child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, baseDofPtr); + child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, + baseDofPtr, dofGeoIndex); if (nextBound1.ithObj >= 0) - child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, baseDofPtr); + child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, + baseDofPtr, dofGeoIndex); } } else { // Either the edge is not contained in further refined children, or @@ -451,12 +458,14 @@ namespace AMDiS { if (childFace0 != -1) { nextBound.ithObj = childFace0; - child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, baseDofPtr); + child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } if (childFace1 != -1) { nextBound.ithObj = childFace1; - child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, baseDofPtr); + child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } } else { ElementDofIterator elDofIter(feSpace, true); diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index 089eb842a3b8b4d52e3a4e74b660ab2ebe42509c..65e52af0a196c7c1626af1ad1e2a6bbf9037593e 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -154,7 +154,8 @@ namespace AMDiS { void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr = false) const; + bool baseDofPtr = false, + vector<GeoIndex>* dofGeoIndex = NULL) const; void prepareNextBound(BoundaryObject &bound, int ithChild) const; diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc index 74c3a9b9400385dd66049c91028663b0bde62cd3..7c7afe0d65d77686215faf5e40aee02c8afac7ef 100644 --- a/AMDiS/src/Triangle.cc +++ b/AMDiS/src/Triangle.cc @@ -115,16 +115,20 @@ namespace AMDiS { if (bound.reverseMode) { nextBound.ithObj = 1; - child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 0; - child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); } else { nextBound.ithObj = 0; - child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 1; - child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); } } } @@ -136,16 +140,20 @@ namespace AMDiS { if (bound.reverseMode) { nextBound.ithObj = 1; - child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 0; - child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); } else { nextBound.ithObj = 0; - child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 1; - child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs, + baseDofPtr); } } } @@ -156,16 +164,16 @@ namespace AMDiS { if (bound.reverseMode) { nextBound.ithObj = 1; - child[1]->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 0; - child[0]->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); } else { nextBound.ithObj = 0; - child[0]->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 1; - child[1]->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); } } break; @@ -178,7 +186,8 @@ namespace AMDiS { void Triangle::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr) const + bool baseDofPtr, + vector<GeoIndex>* dofGeoIndex) const { FUNCNAME("Triange::getHigherOrderDofs()"); @@ -194,7 +203,8 @@ namespace AMDiS { case 0: if (child[1]) { nextBound.ithObj = 2; - child[1]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } else { addThisEdge = true; } @@ -203,7 +213,8 @@ namespace AMDiS { case 1: if (child[0]) { nextBound.ithObj = 2; - child[0]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } else { addThisEdge = true; } @@ -213,14 +224,18 @@ namespace AMDiS { if (child[0]) { if (bound.reverseMode) { nextBound.ithObj = 1; - child[1]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); nextBound.ithObj = 0; - child[0]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } else { nextBound.ithObj = 0; - child[0]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); nextBound.ithObj = 1; - child[1]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, + baseDofPtr, dofGeoIndex); } } else { addThisEdge = true; @@ -238,24 +253,31 @@ namespace AMDiS { if (baseDofPtr) { do { - if (elDofIter.getCurrentPos() == 1 && + if (elDofIter.getPosIndex() == EDGE && elDofIter.getCurrentElementPos() == bound.ithObj) addDofs.push_back(elDofIter.getBaseDof()); } while (elDofIter.nextStrict()); } else { do { - if (elDofIter.getCurrentPos() == 1 && + if (elDofIter.getPosIndex() == EDGE && elDofIter.getCurrentElementPos() == bound.ithObj) - addDofs.push_back(elDofIter.getDofPtr()); + addDofs.push_back(elDofIter.getDofPtr()); } while (elDofIter.next()); } - if (bound.reverseMode) - for (int i = addDofs.size() - 1; i >= 0; i--) + if (bound.reverseMode) { + for (int i = addDofs.size() - 1; i >= 0; i--) { dofs.push_back(addDofs[i]); - else - for (unsigned int i = 0; i < addDofs.size(); i++) - dofs.push_back(addDofs[i]); + if (dofGeoIndex != NULL) + dofGeoIndex->push_back(EDGE); + } + } else { + for (unsigned int i = 0; i < addDofs.size(); i++) { + dofs.push_back(addDofs[i]); + if (dofGeoIndex != NULL) + dofGeoIndex->push_back(EDGE); + } + } } } diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index 85da4d6287b4cfbf69158d9bca1f90c446fe13ee..24e7cccb9d7e10da07cc6909cc244fabb7589f68 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -202,7 +202,8 @@ namespace AMDiS { void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, DofContainer& dofs, - bool baseDofPtr = false) const; + bool baseDofPtr = false, + vector<GeoIndex>* dofGeoIndex = NULL) const; protected: /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index c658b126cdad765ee1b840a793578b5505605d6b..2d67e07a0df5ff1f30572cc8483ddc219fd78e05 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -1086,16 +1086,11 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::repartitionMesh()"); - TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1) - ("Only meshes with one DOFAdmin are supported!\n"); - - // === First we check if the rank with the maximum number of DOFs has at === // === least 20% more DOFs than the rank with the minimum number of DOFs. === // === In this case, the mesh will be repartition. === int repartitioning = 0; - vector<int> nDofsInRank(mpiSize); int nDofs = mesh->getDofAdmin(0).getUsedDofs(); mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0); @@ -1121,7 +1116,6 @@ namespace AMDiS { mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0); } - if (repartitioning == 0) return; @@ -1199,6 +1193,8 @@ namespace AMDiS { // === Add new macro elements to mesh. === + TEST_EXIT(mesh->getGeo(FACE) || mesh->getGeo(CENTER)) + ("Mh, dass macht dann doch noch etwas Arbeit f�r den Thomas!\n"); for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); it != newMacroEl.end(); ++it) { @@ -1244,7 +1240,6 @@ namespace AMDiS { } } - StdMpi<MeshCodeVec> stdMpi(mpiComm, true); stdMpi.send(sendCodes); for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin(); @@ -1252,7 +1247,6 @@ namespace AMDiS { stdMpi.recv(it->first); stdMpi.startCommunication(); - if (interchangeVectors.size() == 0) WARNING("There are no interchange vectors defined!\n"); @@ -1368,9 +1362,9 @@ namespace AMDiS { MSG("Debug mode tests finished!\n"); #endif - // === In 3D we have to make some test, if the resulting mesh is valid. If === - // === it is not valid, there is no possiblity yet to fix this problem, just === - // === exit with an error message. === + // === In 3D we have to make some test, if the resulting mesh is valid. === + // === If it is not valid, there is no possiblity yet to fix this === + // === problem, just exit with an error message. === check3dValidMesh(); @@ -1516,9 +1510,11 @@ namespace AMDiS { AtomicBoundary& b = rankIntBoundary.getNewAtomic(it2->first); b = bound; if (geoIndex == EDGE) - b.neighObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, it2->second); + b.neighObj.reverseMode = + elObjects.getEdgeReverseMode(rankBoundEl, it2->second); if (geoIndex == FACE) - b.neighObj.reverseMode = elObjects.getFaceReverseMode(rankBoundEl, it2->second); + b.neighObj.reverseMode = + elObjects.getFaceReverseMode(rankBoundEl, it2->second); } } else { @@ -1538,9 +1534,11 @@ namespace AMDiS { AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner); b = bound; if (geoIndex == EDGE) - b.rankObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, ownerBoundEl); + b.rankObj.reverseMode = + elObjects.getEdgeReverseMode(rankBoundEl, ownerBoundEl); if (geoIndex == FACE) - b.rankObj.reverseMode = elObjects.getFaceReverseMode(rankBoundEl, ownerBoundEl); + b.rankObj.reverseMode = + elObjects.getFaceReverseMode(rankBoundEl, ownerBoundEl); } } } diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc index 0e75b3633a36bf1ccc3b799aa8cd25ed449ccfad..5583cd6b52a5ae272a09b6bb9d7c9890de2b9e9e 100644 --- a/AMDiS/src/parallel/MeshManipulation.cc +++ b/AMDiS/src/parallel/MeshManipulation.cc @@ -131,7 +131,7 @@ namespace AMDiS { if (macrosProcessed.count(elIt->elIndex) == 1) { TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex)) ("Should not happen!\n"); - + Element *el0 = (*it)->getElement(); Element *el1 = macroIndexMap[elIt->elIndex]->getElement(); @@ -141,17 +141,24 @@ namespace AMDiS { BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false); DofContainer dofs0, dofs1; - - el0->getAllDofs(feSpace, b0, dofs0); - el1->getAllDofs(feSpace, b1, dofs1); + vector<GeoIndex> dofGeoIndex0, dofGeoIndex1; + el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0); + el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1); #if (DEBUG != 0) - debug::testDofsByCoords(feSpace, dofs0, dofs1); + if (feSpaces.size() == 1) + debug::testDofsByCoords(feSpace, dofs0, dofs1); + else + TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n"); #endif for (unsigned int i = 0; i < dofs0.size(); i++) { mapDelDofs[dofs0[i]] = dofs1[i]; - dofPosIndex[dofs0[i]] = EDGE; + dofPosIndex[dofs0[i]] = dofGeoIndex0[i]; + + TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i]) + ("Should not happen: %d %d\n", + dofGeoIndex0[i], dofGeoIndex1[i]); } break; @@ -182,10 +189,9 @@ namespace AMDiS { BoundaryObject b0(el0, 0, FACE, i, reverseMode); BoundaryObject b1(el1, 0, FACE, elIt->ithObject, false); - DofContainer dofs0, dofs1; - - el0->getAllDofs(feSpace, b0, dofs0); - el1->getAllDofs(feSpace, b1, dofs1); + DofContainer dofs0, dofs1; + el0->getAllDofs(feSpace, b0, dofs0, true); + el1->getAllDofs(feSpace, b1, dofs1, true); #if (DEBUG != 0) debug::testDofsByCoords(feSpace, dofs0, dofs1); @@ -240,7 +246,8 @@ namespace AMDiS { for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin(); it != mapDelDofs.end(); ++it) - mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX); + mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), + dofPosIndex[it->first]); }