Small improvements in parallel code.

AMDiS/src/DualTraverse.cc

@@ -208,14 +208,9 @@ namespace AMDiS {
     if (largeElInfo->getLevel() == smallElInfo->getLevel())
       return largeFace;
 
-    int refinementPathLength = smallElInfo->getRefinementPathLength();
-    unsigned long refinementPath = smallElInfo->getRefinementPath();
-
-    TEST_EXIT_DBG(refinementPathLength != 0)
+    TEST_EXIT_DBG(smallElInfo->getRefinementPathLength() != 0)
       ("Refinement path is empty. This should not happen!\n");
 
-    std::cout << "REFINEMENTPATHLENGTH = " << refinementPathLength << std::endl;
-
     return -1;
   }
 }

AMDiS/src/Traverse.cc

@@ -39,7 +39,7 @@ namespace AMDiS {
     traverse_fill_flag = fill_flag; 
 
     TEST_EXIT_DBG(mesh)("No mesh!\n");
-    TEST_EXIT_DBG(traverse_mesh->getMacroElements().size() > 0)("Mesh is empty!\n");
+    TEST_EXIT(traverse_mesh->getMacroElements().size() > 0)("Mesh is empty!\n");
 
     if (stack_size < 1)
       enlargeTraverseStack();

AMDiS/src/parallel/ElementObjectData.cc

@@ -422,10 +422,14 @@ namespace AMDiS {
   {
     FUNCNAME("ElementObjects::createReverseModeData()");
 
-    // In 2D, all reverse modes are always true!
+    // === In 2D, all reverse modes are always true! ===
+
     if (mesh->getDim() == 2)
       return;
 
+
+    // === First, create reverse modes for all "directly" neighbouring elements. ===
+
     for (map<DofEdge, vector<ElementObjectData> >::iterator edgeIt = edgeElements.begin();
 	 edgeIt != edgeElements.end(); ++edgeIt) {
 
@@ -473,6 +477,54 @@ namespace AMDiS {
 	}
       }
     }
+
+
+    // === And create reverse modes for periodic neighbouring elements. ===
+
+    for (PerBoundMap<DofEdge>::iterator edgeIt = periodicEdges.begin();
+	 edgeIt != periodicEdges.end(); ++edgeIt) {
+      vector<ElementObjectData> &edges0 = edgeElements[edgeIt->first.first];
+      vector<ElementObjectData> &edges1 = edgeElements[edgeIt->first.second];
+
+      for (unsigned int i = 0; i < edges0.size(); i++) {
+	BoundaryObject obj0(elIndexMap[edges0[i].elIndex], 
+			    elIndexTypeMap[edges0[i].elIndex], 
+			    EDGE, edges0[i].ithObject);
+
+	for (unsigned int j = 0; j < edges1.size(); j++) {
+	  BoundaryObject obj1(elIndexMap[edges1[j].elIndex], 
+			      elIndexTypeMap[edges1[j].elIndex], 
+			      EDGE, edges1[j].ithObject);
+
+	  bool reverseMode = 
+	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, edgeIt->second);
+
+	  edgeReverseMode[make_pair(edges0[i], edges1[j])] = reverseMode;
+	  edgeReverseMode[make_pair(edges1[j], edges0[i])] = reverseMode;
+	}
+      }
+    }
+
+    for (PerBoundMap<DofFace>::iterator faceIt = periodicFaces.begin();
+	 faceIt != periodicFaces.end(); ++faceIt) {
+      vector<ElementObjectData> &faces0 = faceElements[faceIt->first.first];
+      vector<ElementObjectData> &faces1 = faceElements[faceIt->first.second];
+
+      TEST_EXIT_DBG(faces0.size() == faces1.size() == 1)("Should not happen!\n");
+
+      BoundaryObject obj0(elIndexMap[faces0[0].elIndex], 
+			  elIndexTypeMap[faces0[0].elIndex], 
+			  FACE, faces0[0].ithObject);      
+      BoundaryObject obj1(elIndexMap[faces1[0].elIndex], 
+			  elIndexTypeMap[faces1[0].elIndex], 
+			  FACE, faces1[0].ithObject);
+      
+      bool reverseMode = 
+	BoundaryObject::computeReverseMode(obj0, obj1, feSpace, faceIt->second);
+      
+      faceReverseMode[make_pair(faces0[0], faces1[0])] = reverseMode;
+      faceReverseMode[make_pair(faces1[0], faces0[0])] = reverseMode;
+    }  
   }
 
 

AMDiS/src/parallel/InteriorBoundary.cc

@@ -92,16 +92,6 @@ namespace AMDiS {
   }
 
 
-  void BoundaryObject::setReverseMode(BoundaryObject &otherBound, 
-				      FiniteElemSpace *feSpace,
-				      BoundaryType boundary)
-  {
-    FUNCNAME("BoundaryObject::setReverseMode()");
-
-    otherBound.reverseMode = computeReverseMode(*this, otherBound, feSpace, boundary);
-  }
-
-
   bool BoundaryObject::operator==(const BoundaryObject& other) const
   {
     return (other.elIndex == elIndex && 

AMDiS/src/parallel/InteriorBoundary.h

@@ -59,10 +59,6 @@ namespace AMDiS {
 				   FiniteElemSpace *feSpace,
 				   BoundaryType boundary);
 
-    void setReverseMode(BoundaryObject &otherBound, 
-			FiniteElemSpace *feSpace,
-			BoundaryType boundary);
-
     bool operator==(const BoundaryObject& other) const;
     
     bool operator!=(const BoundaryObject& other) const;

AMDiS/src/parallel/MeshDistributor.cc

@@ -119,9 +119,9 @@ namespace AMDiS {
     TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
     TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
 
-    int a;
-    char *b;
-    float zoltanVersion;
+    int a = 0;
+    char *b = NULL;
+    float zoltanVersion = 0.0;
     Zoltan_Initialize(a, &b, &zoltanVersion);
 
     elObjects.setFeSpace(feSpace);
@@ -165,12 +165,9 @@ namespace AMDiS {
     setInitialElementWeights();
 
     // and now partition the mesh    
-    oldPartitionMap = partitioner->getPartitionMap();
-
     bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
     TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
-
-    partitionMap = partitioner->getPartitionMap();
+    partitioner->createPartitionMap(partitionMap);
 
 
 #if (DEBUG != 0)
@@ -742,6 +739,25 @@ namespace AMDiS {
       repartitionMesh();
       nMeshChangesAfterLastRepartitioning = 0;
     }
+
+    // === Print imbalance factor. ===
+
+    vector<int> nDofsInRank(mpiSize);
+    int nDofs = mesh->getDofAdmin(0).getUsedDofs();
+    mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0);
+
+    if (mpiRank == 0) {
+      int nOverallDofs = 0;
+      int maxDofs = numeric_limits<int>::min();
+      for (int i = 0; i < mpiSize; i++) {
+	nOverallDofs += nDofsInRank[i];
+	maxDofs = std::max(maxDofs, nDofsInRank[i]);
+      }      
+      int avrgDofs = nOverallDofs / mpiSize;
+      double imbalance = (static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;
+
+      MSG("Imbalancing factor: %.1f\%\n", imbalance);
+    }
   }
 
   
@@ -935,7 +951,6 @@ namespace AMDiS {
       return;
 
     double timePoint = MPI::Wtime();
-    repartitioningCounter++;
 
 #if (DEBUG != 0)
     ParallelDebug::testDoubleDofs(mesh);
@@ -944,6 +959,7 @@ namespace AMDiS {
       ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", 
 					   repartitioningCounter, feSpace);
 #endif
+    repartitioningCounter++;
 
     // === Create new element weights. ===
 
@@ -964,8 +980,6 @@ namespace AMDiS {
       MSG("Mesh partitioner created empty partition!\n");
       return;
     }
-    oldPartitionMap = partitionMap;
-
 
     // === Create map that maps macro element indices to pointers to the ===
     // === macro elements.                                               ===
@@ -1051,6 +1065,9 @@ namespace AMDiS {
       stdMpi.recv(it->first);
     stdMpi.startCommunication();
 
+    if (interchangeVectors.size() == 0)
+      WARNING("There are no interchange vectors defined!\n");
+
     StdMpi<vector<vector<double> > > stdMpi2(mpiComm, true);
     stdMpi2.send(sendValues);
     for (std::map<int, std::vector<int> >::iterator it = partitioner->getRecvElements().begin();
@@ -1058,6 +1075,7 @@ namespace AMDiS {
       stdMpi2.recv(it->first);
     stdMpi2.startCommunication();
 
+
     // === Adapte the new macro elements due to the received mesh structure codes. ===
 
     for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
@@ -1132,17 +1150,15 @@ namespace AMDiS {
     // update the number of elements, vertices, etc. of the mesh.
     mesh->removeMacroElements(deleteMacroElements, feSpace);
 
-
     // === Remove double DOFs. ===
 
     MeshManipulation meshManipulation(feSpace);
     meshManipulation.deleteDoubleDofs(newMacroEl, elObjects);
 
     mesh->dofCompress();
-    partitionMap = partitioner->getPartitionMap();
+    partitioner->createPartitionMap(partitionMap);
 
     updateInteriorBoundaryInfo();
-
     updateLocalGlobalNumbering();
 
 #if (DEBUG != 0)
@@ -1215,11 +1231,11 @@ namespace AMDiS {
     }
 
 
-    elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap);
-
     // === Create periodic data, if there are periodic boundary conditions. ===
     elObjects.createPeriodicData();
 
+    // === Create data about the reverse modes of neighbouring elements. ===
+    elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap);
 
     // === Create mesh element data for this rank. ===
     elObjects.createRankData(partitionMap);
@@ -1283,7 +1299,6 @@ namespace AMDiS {
 	      
 	      AtomicBoundary& b = myIntBoundary.getNewAtomic(it2->first);
 	      b = bound;
-	      //	      b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
 	      if (geoIndex == EDGE)
 		b.neighObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, it2->second);
 	      if (geoIndex == FACE)
@@ -1306,9 +1321,6 @@ namespace AMDiS {
 	    
 	    AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner);
 	    b = bound;	    
-
-	    //	    b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
-	    
 	    if (geoIndex == EDGE)
 	      b.rankObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
 	    if (geoIndex == FACE)
@@ -1391,9 +1403,9 @@ namespace AMDiS {
 	b = bound;
      
 	if (mpiRank > otherElementRank)
-	  b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+	  b.neighObj.reverseMode = elObjects.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
 	else
-	  b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+	  b.rankObj.reverseMode = elObjects.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
       }
     }
 
@@ -1408,26 +1420,26 @@ namespace AMDiS {
       TEST_EXIT_DBG(elObjects.getElements(it->first.second).size() == 1)
  	("Should not happen!\n");
 
-      ElementObjectData& perEdgeEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
+      ElementObjectData& perFaceEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
 
       for (map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
  	   elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
       
 	int otherElementRank = elIt->first;
-	ElementObjectData& perEdgeEl1 = elIt->second;
+	ElementObjectData& perFaceEl1 = elIt->second;
 
 	AtomicBoundary bound;	    	    
-	bound.rankObj.el = macroElIndexMap[perEdgeEl0.elIndex];
-	bound.rankObj.elIndex = perEdgeEl0.elIndex;
-	bound.rankObj.elType = macroElIndexTypeMap[perEdgeEl0.elIndex];
+	bound.rankObj.el = macroElIndexMap[perFaceEl0.elIndex];
+	bound.rankObj.elIndex = perFaceEl0.elIndex;
+	bound.rankObj.elType = macroElIndexTypeMap[perFaceEl0.elIndex];
 	bound.rankObj.subObj = FACE;
-	bound.rankObj.ithObj = perEdgeEl0.ithObject;
+	bound.rankObj.ithObj = perFaceEl0.ithObject;
 	
-	bound.neighObj.el = macroElIndexMap[perEdgeEl1.elIndex];
-	bound.neighObj.elIndex = perEdgeEl1.elIndex;
-	bound.neighObj.elType = macroElIndexTypeMap[perEdgeEl1.elIndex];
+	bound.neighObj.el = macroElIndexMap[perFaceEl1.elIndex];
+	bound.neighObj.elIndex = perFaceEl1.elIndex;
+	bound.neighObj.elType = macroElIndexTypeMap[perFaceEl1.elIndex];
 	bound.neighObj.subObj = FACE;
-	bound.neighObj.ithObj = perEdgeEl1.ithObject;
+	bound.neighObj.ithObj = perFaceEl1.ithObject;
 	
 	bound.type = it->second;
 	
@@ -1435,9 +1447,9 @@ namespace AMDiS {
 	b = bound;
      
 	if (mpiRank > otherElementRank)
-	  b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+	  b.neighObj.reverseMode = elObjects.getFaceReverseMode(perFaceEl0, perFaceEl1);
 	else
-	  b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+	  b.rankObj.reverseMode = elObjects.getFaceReverseMode(perFaceEl0, perFaceEl1);
       }
     }
     
@@ -1952,7 +1964,6 @@ namespace AMDiS {
 
     SerUtil::serialize(out, elemWeights);
     SerUtil::serialize(out, partitionMap);
-    SerUtil::serialize(out, oldPartitionMap);
     SerUtil::serialize(out, nRankDofs);
     SerUtil::serialize(out, nOverallDofs);
 
@@ -1993,7 +2004,6 @@ namespace AMDiS {
 
     SerUtil::deserialize(in, elemWeights);
     SerUtil::deserialize(in, partitionMap);
-    SerUtil::deserialize(in, oldPartitionMap);
     SerUtil::deserialize(in, nRankDofs);    
     SerUtil::deserialize(in, nOverallDofs);
 

AMDiS/src/parallel/MeshDistributor.h

@@ -31,6 +31,7 @@
 
 #include "parallel/MeshPartitioner.h"
 #include "parallel/InteriorBoundary.h"
+#include "AMDiS_fwd.h"
 #include "Global.h"
 #include "ProblemTimeInterface.h"
 #include "ProblemIterationInterface.h"
@@ -38,7 +39,7 @@
 #include "Serializer.h"
 #include "BoundaryManager.h"
 #include "ElementObjectData.h"
-#include "AMDiS_fwd.h"
+#include "SystemVector.h"
 
 namespace AMDiS {
 
@@ -94,6 +95,13 @@ namespace AMDiS {
       interchangeVectors.push_back(vec);
     }
 
+    /// Adds all DOFVectors of a SystemVector to \ref interchangeVecs.
+    void addInterchangeVector(SystemVector *vec)
+    {
+      for (int i = 0; i < vec->getSize(); i++)
+	interchangeVectors.push_back(vec->getDOFVector(i));
+    }
+    
     /** \brief
      * This function checks if the mesh has changed on at least on rank. In this case,
      * the interior boundaries are adapted on all ranks such that they fit together on
@@ -448,12 +456,6 @@ namespace AMDiS {
      */
     map<int, int> partitionMap;
 
-    /** \brief
-     * Stores an old partitioning of elements. To every macro element index the
-     * number of the rank it corresponds to is stored.
-     */
-    map<int, int> oldPartitionMap; 
-   
     /// Number of DOFs in the rank mesh.
     int nRankDofs;
 

AMDiS/src/parallel/MeshPartitioner.h

@@ -58,6 +58,8 @@ namespace AMDiS {
     virtual bool partition(map<int, double> &elemWeights,
 			   PartitionMode mode = INITIAL) = 0;
 
+    virtual void createPartitionMap(map<int, int>& partitionMap) = 0;
+
     /// Write partitioner state to disk.
     void serialize(ostream &out);
 
@@ -84,11 +86,6 @@ namespace AMDiS {
       return elementInRank;
     }
 
-    map<int, int>& getPartitionMap()
-    {
-      return partitionMap;
-    }
-
     map<int, vector<int> >& getRecvElements()
     {
       return recvElements;

AMDiS/src/parallel/ParMetisPartitioner.cc

@@ -421,16 +421,11 @@ namespace AMDiS {
 
     // === Distribute new partition data. ===
 
-    bool b = distributePartitioning(&(part[0]));
-
-    if (b)
-      createPartitionMap();
-
-    return b;
+    return distributePartitioning(&(part[0]));
   }
 
 
-  void ParMetisPartitioner::createPartitionMap()
+  void ParMetisPartitioner::createPartitionMap(map<int, int>& pMap)
   {
     FUNCNAME("ParMetisPartitioner::createPartitionMap()");
 
@@ -468,6 +463,8 @@ namespace AMDiS {
     for (int i = 0; i < mpiSize; i++)
       for (int j = 0; j < nPartitionElements[i]; j++)
 	partitionMap[partitionElements[elmdist[i] + j]] = i;
+
+    pMap = partitionMap;