Bugfix in parallel code.

AMDiS/src/InteriorBoundary.cc

@@ -18,7 +18,9 @@ namespace AMDiS {
     
     switch (feSpace->getMesh()->getDim()) {
     case 2:
-      ERROR_EXIT("Not yet implemented!\n");
+      if (ithObj == 2 || ithObj != otherBound.ithObj)
+	otherMode = true;
+
       break;
     case 3:
       if (ithObj == 2 || ithObj == 3) {

AMDiS/src/InteriorBoundary.h

@@ -116,6 +116,82 @@ namespace AMDiS {
   public:
     typedef std::map<int, std::vector<AtomicBoundary> > RankToBoundMap;
 
+    /// Iterator for the interior boundary object.
+    class iterator {      
+    public:
+      iterator(InteriorBoundary &b)
+	: bound(b)
+      {
+	reset();
+      }
+
+      /// Set the iterator to the first position.
+      void reset()
+      {
+	mapIt = bound.boundary.begin();
+	nextNonempty();
+
+	if (mapIt != bound.boundary.end())
+	  vecIt = mapIt->second.begin();	
+      }
+
+      /// Test if iterator is at the final position.
+      bool end() const
+      {
+	return (mapIt == bound.boundary.end());
+      }
+
+      /// Move iterator to the next position.
+      void operator++()
+      {
+	++vecIt;
+	if (vecIt == mapIt->second.end()) {
+	  ++mapIt;
+	  nextNonempty();
+
+	  if (mapIt != bound.boundary.end())
+	    vecIt = mapIt->second.begin();	  
+	}	
+      }
+
+      void nextRank()
+      {
+	++mapIt;
+	nextNonempty();
+
+	if (mapIt != bound.boundary.end())
+	  vecIt = mapIt->second.begin();	
+      }
+
+      int getRank() 
+      {
+	return mapIt->first;
+      }
+
+      AtomicBoundary& getBound()
+      {
+	return *vecIt;
+      }
+
+    protected:
+
+      inline void nextNonempty()
+      {
+ 	while (mapIt->second.size() == 0) {
+ 	  ++mapIt;
+	  if (mapIt == bound.boundary.end())
+	    return;
+	}
+      }
+
+    protected:
+      RankToBoundMap::iterator mapIt;
+
+      std::vector<AtomicBoundary>::iterator vecIt;
+
+      InteriorBoundary &bound;
+    };
+
   public:
     InteriorBoundary() {}
 

AMDiS/src/ParallelDomainBase.cc

@@ -730,6 +730,8 @@ namespace AMDiS {
   {
     FUNCNAME("ParallelDomainBase::checkMeshChange()");
 
+    int dim = mesh->getDim();
+
     // === If mesh has not been changed on all ranks, return. ===
 
     int recvAllValues = 0;
@@ -748,19 +750,15 @@ namespace AMDiS {
       // To check the interior boundaries, the ownership of the boundaries is not 
       // important. Therefore, we add all boundaries to one boundary container.
       RankToBoundMap allBound;
-      for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); 
-	   it != myIntBoundary.boundary.end(); ++it)
-	for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
-	     bit != it->second.end(); ++bit)
-	  if (bit->rankObj.subObj == FACE)
-	    allBound[it->first].push_back(*bit);
-
-      for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); 
-	   it != otherIntBoundary.boundary.end(); ++it)
-	for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
-	     bit != it->second.end(); ++bit)
-	  if (bit->rankObj.subObj == FACE)
-	    allBound[it->first].push_back(*bit);
+
+      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
+	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
+	  allBound[it.getRank()].push_back(it.getBound());
+
+      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
+	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
+	  allBound[it.getRank()].push_back(it.getBound());
+
 
       // === Check the boundaries and adapt mesh if necessary. ===
 
@@ -771,7 +769,6 @@ namespace AMDiS {
       int sendValue = static_cast<int>(!meshChanged);
       recvAllValues = 0;
       mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
-      MSG("LOOP\n");
     } while (recvAllValues != 0);
 
 
@@ -805,8 +802,6 @@ namespace AMDiS {
   {
     FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");
 
-    MSG("CHECK_AND_ADAPT 1!\n");
-
     // === Create mesh structure codes for all ranks boundary elements. ===
        
     std::map<int, MeshCodeVec> sendCodes;
@@ -820,21 +815,14 @@ namespace AMDiS {
       }
     }
 
-    MSG("CHECK_AND_ADAPT 2!\n");
-
     StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
     stdMpi.send(sendCodes);
     stdMpi.recv(allBound);
     stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
  
-    MSG("CHECK_AND_ADAPT 3!\n");
-
-    clock_t first = clock();
-
     // === Compare received mesh structure codes. ===
     
     bool meshFitTogether = true;
-    int superCounter = 0;
 
     for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
       
@@ -849,8 +837,6 @@ namespace AMDiS {
 	
 	if (elCode.getCode() != recvCodes[i].getCode()) {
 	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
-	  superCounter++;
-	  clock_t fff = clock();
 
 	  int refCounter = 0;
 	  bool b = fitElementToMeshCode(recvCodes[i], 
@@ -858,12 +844,6 @@ namespace AMDiS {
 					boundIt->rankObj.ithObj, 
 					boundIt->rankObj.elType, refCounter);  
 
-	  MSG("SIZE OF ELINFO3d = %d\n", sizeof(ElInfo3d));
-	  MSG("Code length %d with %d refs needed %.5f seconds\n", 
-	      recvCodes[i].getNumElements(), 
-	      refCounter,
-	      TIME_USED(fff, clock()));
-
 	  if (b)
 	    meshFitTogether = false;
 	}
@@ -872,9 +852,6 @@ namespace AMDiS {
       }
     }
 
-    MSG("time for %d needed %.5f seconds\n", superCounter, TIME_USED(first, clock()));
-    MSG("CHECK_AND_ADAPT 4!\n");
-
     return meshFitTogether;
   }
 
@@ -1370,77 +1347,69 @@ namespace AMDiS {
     std::set<DegreeOfFreedom> rankBoundaryDofs;
 
     // First, traverse the rank owned elements af the interior boundaries.
-    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
-	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {
-      for (std::vector<AtomicBoundary>::iterator boundIt = rankIt->second.begin(); 
-	   boundIt != rankIt->second.end(); ++boundIt) {
-	Element *el = boundIt->rankObj.el;
-	basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
-
-	if (dim == 3) {
-	  for (int j = 0; j < 3; j++) {
-	    int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j);
-	    DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
-	    DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
-	    GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
-	    
-	    // If the edge at rank's interior boundarie has a higher owner rank, than
-	    // we have to remove this edge from the corresponding boundary element.
-	    // Otherwise, it is part of the interior boundary and we add it to the set
-	    // rankBoundaryEdges.
-	    if (edgeOwner[edge] > mpiRank)
-	      boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
-	    else
-	      rankBoundaryEdges.insert(edge);
-	  }
-	}
+    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
+      Element *el = it.getBound().rankObj.el;
+      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
       
+      if (dim == 3) {
 	for (int j = 0; j < 3; j++) {
-	  int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j);
-	  DegreeOfFreedom dof = localIndices[dofNo];
-
-	  if (dofOwner[dof] > mpiRank)
-	    boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
+	  int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j);
+	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
+	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
+	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
+	  
+	  // If the edge at rank's interior boundarie has a higher owner rank, than
+	  // we have to remove this edge from the corresponding boundary element.
+	  // Otherwise, it is part of the interior boundary and we add it to the set
+	  // rankBoundaryEdges.
+	  if (edgeOwner[edge] > mpiRank)
+	    it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
 	  else
-	    rankBoundaryDofs.insert(dof);
+	    rankBoundaryEdges.insert(edge);
 	}
       }
+      
+      for (int j = 0; j < 3; j++) {
+	int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j);
+	DegreeOfFreedom dof = localIndices[dofNo];
+	
+	if (dofOwner[dof] > mpiRank)
+	  it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
+	else
+	  rankBoundaryDofs.insert(dof);
+      }    
     }
 
-    // Now, do the same with all other elements at the interio boundaries.
-    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
-	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
-      for (std::vector<AtomicBoundary>::iterator boundIt = rankIt->second.begin(); 
-	   boundIt != rankIt->second.end(); ++boundIt) {
-	Element *el = boundIt->rankObj.el;
-	basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
-
+    // Now, do the same with all other elements at the interior boundaries.
+    for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
+      Element *el = it.getBound().rankObj.el;
+      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
+      
+      if (dim == 3) {
 	for (int j = 0; j < 3; j++) {
-	  int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j);
+	  int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j);
 	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
 	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
 	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
-
-	  if (edgeOwner[edge] > rankIt->first)
-	    boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
-	  else
-	    rankBoundaryEdges.insert(edge);
-
-	}
-
-	for (int j = 0; j < 3; j++) {
-	  int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j);
-	  DegreeOfFreedom dof = localIndices[dofNo];
 	  
-	  if (dofOwner[dof] > rankIt->first)
-	    boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
+	  if (edgeOwner[edge] > it.getRank())
+	    it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
 	  else
-	    rankBoundaryDofs.insert(dof);
-	}	
+	    rankBoundaryEdges.insert(edge);	    
+	}
       }
+      
+      for (int j = 0; j < 3; j++) {
+	int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j);
+	DegreeOfFreedom dof = localIndices[dofNo];
+	
+	if (dofOwner[dof] > it.getRank())
+	  it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
+	else
+	  rankBoundaryDofs.insert(dof);
+      }	      
     }
 
-
     // === Create the new interior boundaries consisting only of edges. This    ===
     // === boundaries are created on that ranks, which do not own the boundary  ===
     // === but are on the other side of the edge. Than, theses ranks inform the ===
@@ -1894,45 +1863,29 @@ namespace AMDiS {
       }
     }
 
-    
-    for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
-	 it != myIntBoundary.boundary.end(); ++it) {    
-     
-      DofContainer &dofsToSend = sendDofs[it->first];
-
-      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
-	   boundIt != it->second.end(); ++boundIt) {
-	DofContainer dofs;	
-	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
-   	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
-
-	for (int i = 0; i < static_cast<int>(dofs.size()); i++)
-	  dofsToSend.push_back(dofs[i]);
-      }
-    } 
-
-    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
-	 it != otherIntBoundary.boundary.end(); ++it) {
-
-      DofContainer &dofsToRecv = recvDofs[it->first];
-
-      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
-	   boundIt != it->second.end(); ++boundIt) {
-	DofContainer dofs;
-   	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
-	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
-
-	for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
-	  DofContainer::iterator eraseIt = 
-	    find(rankDofs.begin(), rankDofs.end(), dofs[i]);
-	  if (eraseIt != rankDofs.end()) 
-	    rankDofs.erase(eraseIt);	 
-  
-	  dofsToRecv.push_back(dofs[i]);
-	}
-      }
+    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
+      DofContainer dofs;	
+      it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs);
+      it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().rankObj, dofs);
+      for (int i = 0; i < static_cast<int>(dofs.size()); i++)
+	sendDofs[it.getRank()].push_back(dofs[i]);
     }
 
+    for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
+      DofContainer dofs;
+      it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().rankObj, dofs);
+      it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs);
+      
+      for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
+	DofContainer::iterator eraseIt = 
+	  find(rankDofs.begin(), rankDofs.end(), dofs[i]);
+	if (eraseIt != rankDofs.end()) 
+	  rankDofs.erase(eraseIt);	 
+	
+	recvDofs[it.getRank()].push_back(dofs[i]);
+      }      
+    }    
+
     nRankDofs = rankDofs.size();
 
     // === Get starting position for global rank dof ordering. ====