Fixed some bugs for periodic BCs in 3D in parrallel computing.

AMDiS/src/ElInfo.cc

@@ -172,6 +172,7 @@ namespace AMDiS {
     }
   }
 
+
   BoundaryType ElInfo::getBoundary(GeoIndex pos, int i)
   {
     static int indexOffset[3][3] = {

AMDiS/src/parallel/InteriorBoundary.cc

@@ -14,11 +14,13 @@
 #include "FiniteElemSpace.h"
 #include "BasisFunction.h"
 #include "Serializer.h"
+#include "VertexVector.h"
 
 namespace AMDiS {
 
   void BoundaryObject::setReverseMode(BoundaryObject &otherBound, 
-				      FiniteElemSpace *feSpace)
+				      FiniteElemSpace *feSpace,
+				      BoundaryType boundary)
   {
     FUNCNAME("BoundaryObject::setReverseMode()");
 
@@ -31,7 +33,8 @@ namespace AMDiS {
 
     case 3:
       TEST_EXIT_DBG(otherBound.elType == 0)
-	("Only 3D macro elements with level 0 are supported!\n");
+	("Only 3D macro elements with level 0 are supported. This element has level %d!\n", otherBound.elType);
+
 
       if (subObj == EDGE) {	
 	int el0_v0 = el->getVertexOfEdge(ithObj, 0);
@@ -45,15 +48,22 @@ namespace AMDiS {
 	basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
 	basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
 
-	TEST_EXIT_DBG(localDofs0[el0_v0] == localDofs1[el1_v0] ||
-		      localDofs0[el0_v0] == localDofs1[el1_v1])
-	  ("This should not happen!\n");
-	TEST_EXIT_DBG(localDofs0[el0_v1] == localDofs1[el1_v0] ||
-		      localDofs0[el0_v1] == localDofs1[el1_v1])
-	  ("This should not happen!\n");
-	
- 	if (localDofs0[el0_v0] != localDofs1[el1_v0])
-	  otherMode = true; 	
+	Mesh *mesh = feSpace->getMesh();
+
+	if (mesh->isPeriodicAssociation(boundary) == false) {
+	  TEST_EXIT_DBG(localDofs0[el0_v0] == localDofs1[el1_v0] ||
+			localDofs0[el0_v0] == localDofs1[el1_v1])
+	    ("This should not happen!\n");
+	  TEST_EXIT_DBG(localDofs0[el0_v1] == localDofs1[el1_v0] ||
+			localDofs0[el0_v1] == localDofs1[el1_v1])
+	    ("This should not happen!\n");
+	  
+	  if (localDofs0[el0_v0] != localDofs1[el1_v0])
+	    otherMode = true; 	
+	} else {
+	  if (mesh->associated(localDofs0[el0_v0], localDofs1[el1_v0]) == false)
+	    otherMode = true;
+	}
       }
 
       if (subObj == FACE && ithObj != 1) {

AMDiS/src/parallel/InteriorBoundary.h

@@ -54,7 +54,9 @@ namespace AMDiS {
 	excludedSubstructures(0)
     {}
 
-    void setReverseMode(BoundaryObject &otherBound, FiniteElemSpace *feSpace);
+    void setReverseMode(BoundaryObject &otherBound, 
+			FiniteElemSpace *feSpace,
+			BoundaryType boundary);
 
     bool operator==(const BoundaryObject& other) const;
     

AMDiS/src/parallel/MeshDistributor.h

@@ -483,10 +483,6 @@ namespace AMDiS {
     std::map<std::pair<DofEdge, DofEdge>, BoundaryType> periodicEdges;
     std::map<std::pair<DofFace, DofFace>, BoundaryType> periodicFaces;
 
-    // Stores to each DOF all its periodic associations.
-    std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
-
-
     /** \brief 
      * Defines the interior boundaries of the domain that result from partitioning
      * the whole mesh. Contains only the boundaries, which are owned by the rank, i.e.,

AMDiS/src/parallel/ParMetisPartitioner.cc

@@ -364,6 +364,8 @@ namespace AMDiS {
 
     for (int i = 0; i < nElements; i++) {
       wgts[i] = min(wgts[i], kpartMax);
+      wgts[i] = max(wgts[i], 1);
+
       smin = std::min(smin, wgts[i]);
       smax = std::max(smax, wgts[i]);
       ssum += wgts[i];

AMDiS/src/parallel/ParallelDebug.cc

@@ -67,8 +67,11 @@ namespace AMDiS {
 
     // === To the last, do the same of periodic boundaries. ===
 
-    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
+    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();       
 	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+      if (rankIt->first == pdb.mpiRank)
+	continue;
+
       int nValues = rankIt->second.size();
       int* sBuffer = new int[nValues];
       for (int i = 0; i < nValues; i++)
@@ -118,11 +121,17 @@ namespace AMDiS {
 
     for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
 	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+      if (rankIt->first == pdb.mpiRank)
+	continue;
+
       for (unsigned int i = 0; i < rankIt->second.size(); i++) {
 	int elIndex1 = recvBuffers[bufCounter][i];
 	int elIndex2 = pdb.periodicBoundary.boundary[rankIt->first][i].neighObj.elIndex;
 
-	TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at periodic boundary!\n");
+	TEST_EXIT(elIndex1 == elIndex2)
+	  ("Wrong element index at periodic boundary el %d with rank %d: %d %d\n", 
+	   pdb.periodicBoundary.boundary[rankIt->first][i].rankObj.elIndex,
+	   rankIt->first, elIndex1, elIndex2);
       }
 
       delete [] recvBuffers[bufCounter++];

AMDiS/src/parallel/ParallelProblemStatBase.cc

@@ -26,14 +26,16 @@ namespace AMDiS {
     ProblemVec::buildAfterCoarsen(adaptInfo, flag, assembleMatrix, assembleVector);
 
     double vm, rss;
-    processMemUsage(vm, rss);
+    processMemUsage(vm, rss);   
+    vm /= 1024.0;
+    rss /= 1024.0;
     
-    MSG("My memory usage is VM = %f    RSS = %f\n", vm, rss);
+    MSG("My memory usage is VM = %f MB    RSS = %f MB\n", vm, rss);
     
     mpi::globalAdd(vm);
     mpi::globalAdd(rss);
-    
-    MSG("Overall memory usage is VM = %f    RSS = %f\n", vm, rss);
+
+    MSG("Overall memory usage is VM = %f MB    RSS = %f MB\n", vm, rss);
   }
   
   

AMDiS/src/parallel/PetscSolver.cc

@@ -128,8 +128,8 @@ namespace AMDiS {
       }
 
 
-      // === Up to now we have assembled on row. Now, the row must be send to the ===
-      // === corresponding rows to the petsc matrix.                              ===
+      // === Up to now we have assembled one row. Now, the row must be send to the ===
+      // === corresponding rows in the petsc matrix.                               ===
      
       if (periodicRow) {
 	// The row dof is periodic, so send dof to all the corresponding rows.
@@ -144,7 +144,8 @@ namespace AMDiS {
 	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
 		     &(cols[0]), &(values[0]), ADD_VALUES);	
  
-	for (std::set<int>::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) {
+	for (std::set<int>::iterator perIt = perAsc.begin(); 
+	     perIt != perAsc.end(); ++perIt) {
 	  std::vector<int> perCols;
 	  perCols.reserve(300);
 	  std::vector<double> perValues;

AMDiS/src/parallel/StdMpi.cc

@@ -60,7 +60,7 @@ namespace AMDiS {
 
   int intSizeOf(std::vector<AtomicBoundary> &data)
   {
-    return data.size() * 3;
+    return data.size() * 6;
   }
 
   int intSizeOf(std::vector<BoundaryObject> &data)
@@ -242,9 +242,12 @@ namespace AMDiS {
   void makeBuf(std::vector<AtomicBoundary> &data, int *buf)
   {
     for (unsigned int i = 0; i < data.size(); i++) {
-      buf[i * 3] = data[i].rankObj.elIndex;
-      buf[i * 3 + 1] = data[i].rankObj.subObj;
-      buf[i * 3 + 2] = data[i].rankObj.ithObj;
+      buf[i * 6] = data[i].rankObj.elIndex;
+      buf[i * 6 + 1] = data[i].rankObj.subObj;
+      buf[i * 6 + 2] = data[i].rankObj.ithObj;
+      buf[i * 6 + 3] = data[i].neighObj.elIndex;
+      buf[i * 6 + 4] = data[i].neighObj.subObj;
+      buf[i * 6 + 5] = data[i].neighObj.ithObj;
     }
   }
 
@@ -253,13 +256,16 @@ namespace AMDiS {
     if (bufSize == 0)
       return;
 
-    TEST_EXIT(bufSize % 3 == 0)("This should not happen!\n");    
+    TEST_EXIT(bufSize % 6 == 0)("This should not happen!\n");    
 
-    data.resize(bufSize / 3);
-    for (int i = 0; i < bufSize / 3; i++) {
-      data[i].rankObj.elIndex = buf[i * 3];
-      data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 3 + 1]);
-      data[i].rankObj.ithObj = buf[i * 3 + 2];
+    data.resize(bufSize / 6);
+    for (int i = 0; i < bufSize / 6; i++) {
+      data[i].rankObj.elIndex = buf[i * 6];
+      data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 6 + 1]);
+      data[i].rankObj.ithObj = buf[i * 6 + 2];
+      data[i].neighObj.elIndex = buf[i * 6 + 3];
+      data[i].neighObj.subObj = static_cast<GeoIndex>(buf[i * 6 + 4]);
+      data[i].neighObj.ithObj = buf[i * 6 + 5];
     }
   }