Work on parallel periodic boundary condition in 3D.

AMDiS/src/parallel/MeshDistributor.cc

@@ -221,6 +221,9 @@ namespace AMDiS {
 
     createPeriodicMap();    
 
+#if (DEBUG != 0)
+    ParallelDebug::testPeriodicBoundary(*this);
+#endif
 
     // === Global refinements. ===
 
@@ -547,7 +550,7 @@ namespace AMDiS {
 
       for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
 	if (it.getRank() == mpiRank) {
-	  //	  ERROR_EXIT("Na, du weisst schon!\n");
+	  WARNING("Na, du weisst schon!\n");
 	} else {
 	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
 	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
@@ -1261,8 +1264,6 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
 
-    MSG("CREATE BOUNDARY INFO!\n");
-
     createMeshElementData();
 
     createBoundaryData();
@@ -1274,8 +1275,6 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::updateInteriorBoundaryInfo()");
 
-    MSG("UPDATE BOUNDARY INFO!\n");
-
     elObjects.createRankData(partitionVec);
 
     createBoundaryData();
@@ -1328,8 +1327,8 @@ namespace AMDiS {
 
 	    periodicEdgeAssoc[edge1].insert(edge2);
 
-	    periodicDofs[std::make_pair(edge1.first, edge2.first)] = boundaryType;
-	    periodicDofs[std::make_pair(edge1.second, edge2.second)] = boundaryType;
+	    periodicVertices[std::make_pair(edge1.first, edge2.first)] = boundaryType;
+	    periodicVertices[std::make_pair(edge1.second, edge2.second)] = boundaryType;
 
 	    periodicDofAssoc[edge1.first].insert(boundaryType);
 	    periodicDofAssoc[edge1.second].insert(boundaryType);
@@ -1351,9 +1350,9 @@ namespace AMDiS {
 
 	    periodicFaces[std::make_pair(face1, face2)] = elInfo->getBoundary(i);
 
-	    periodicDofs[std::make_pair(face1.get<0>(), face2.get<0>())] = boundaryType;
-	    periodicDofs[std::make_pair(face1.get<1>(), face2.get<1>())] = boundaryType;
-	    periodicDofs[std::make_pair(face1.get<2>(), face2.get<2>())] = boundaryType;
+	    periodicVertices[std::make_pair(face1.get<0>(), face2.get<0>())] = boundaryType;
+	    periodicVertices[std::make_pair(face1.get<1>(), face2.get<1>())] = boundaryType;
+	    periodicVertices[std::make_pair(face1.get<2>(), face2.get<2>())] = boundaryType;
 
 	    periodicDofAssoc[face1.get<0>()].insert(boundaryType);
 	    periodicDofAssoc[face1.get<1>()].insert(boundaryType);
@@ -1395,7 +1394,7 @@ namespace AMDiS {
 
     // === Search for interectly connected vertices in periodic boundaries. ===
 
-    if (periodicDofs.size() > 0) {
+    if (periodicVertices.size() > 0) {
 
       
       // === Search for an unsed boundary index. ===
@@ -1413,50 +1412,49 @@ namespace AMDiS {
 
       // === Get all vertex DOFs that have multiple periodic associations. ===
 
-      std::vector<DegreeOfFreedom> multiplePeriodicDof;
+      std::map<int, std::vector<DegreeOfFreedom> > multiplePeriodicDof;
       for (std::map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
 	   it != periodicDofAssoc.end(); ++it) {
-	if (mesh->getDim() == 2) {
-	  TEST_EXIT_DBG(it->second.size() <= 2)("Should not happen!\n");
-	}
-	if (mesh->getDim() == 3) {
-	  TEST_EXIT_DBG(it->second.size() <= 3)("Should not happen!\n");
-	}
-
+	TEST_EXIT_DBG((mesh->getDim() == 2 && it->second.size() <= 2) ||
+		      (mesh->getDim() == 3 && it->second.size() <= 3))
+	  ("Should not happen!\n");
 	
-	if ((mesh->getDim() == 2 && it->second.size() == 2) ||
-	    (mesh->getDim() == 3 && it->second.size() == 3))
-	  multiplePeriodicDof.push_back(it->first);
+	multiplePeriodicDof[it->second.size()].push_back(it->first);
       }
 
 
-      if (mesh->getDim() == 2) {
-	TEST_EXIT_DBG(multiplePeriodicDof.size() == 0 || 
-		      multiplePeriodicDof.size() == 4)
-	  ("Should not happen (%d)!\n", multiplePeriodicDof.size());
-      }
-
-      if (mesh->getDim() == 3) {
-	TEST_EXIT_DBG(multiplePeriodicDof.size() == 0 || 
-		      multiplePeriodicDof.size() == 8)
-	  ("Should not happen (%d)!\n", multiplePeriodicDof.size());
-      }
+      if (mesh->getDim() == 2)
+	TEST_EXIT_DBG(multiplePeriodicDof[2].size() == 0 || 
+		      multiplePeriodicDof[2].size() == 4)
+	  ("Should not happen (%d)!\n", multiplePeriodicDof[2].size());
+      if (mesh->getDim() == 3)
+	TEST_EXIT_DBG(multiplePeriodicDof[3].size() == 0 || 
+		      multiplePeriodicDof[3].size() == 8)
+	  ("Should not happen (%d)!\n", multiplePeriodicDof[3].size());
 
 
-      int nMultiplePeriodicDofs = multiplePeriodicDof.size();
+      for (int k = 2; k <= 3; k++) {
+	int nMultiplePeriodicDofs = multiplePeriodicDof[k].size();
 	for (int i = 0; i < nMultiplePeriodicDofs; i++) {
-	for (int j = 0; j < nMultiplePeriodicDofs; j++) {
-	  if (i == j) 
-	    continue;
-	  
-	  std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs = 
-	    std::make_pair(multiplePeriodicDof[i], multiplePeriodicDof[j]);
+	  for (int j = i + 1; j < nMultiplePeriodicDofs; j++) {
+	    std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs0 = 
+	      std::make_pair(multiplePeriodicDof[k][i], multiplePeriodicDof[k][j]);
+	    std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs1 = 
+	      std::make_pair(multiplePeriodicDof[k][j], multiplePeriodicDof[k][i]);
+	    
+	    if (periodicVertices.count(perDofs0) == 0) {
+	      TEST_EXIT_DBG(periodicVertices.count(perDofs1) == 0)
+		("Should not happen!\n");
 	      
-	  if (periodicDofs.count(perDofs) == 0)
-	    periodicDofs[perDofs] = newPeriodicBoundaryType;	   	    
+	      periodicVertices[perDofs0] = newPeriodicBoundaryType;
+	      periodicVertices[perDofs1] = newPeriodicBoundaryType;
+	      newPeriodicBoundaryType--;
+	      mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
+		new VertexVector(feSpace->getAdmin(), "");
+	    }
+	  }
 	}
       }
-
 
 
       // === Get all edges that have multiple periodic associations (3D only!). ===
@@ -1474,8 +1472,45 @@ namespace AMDiS {
 
 	  periodicEdges[perEdge0] = newPeriodicBoundaryType;
 	  periodicEdges[perEdge1] = newPeriodicBoundaryType;
+	  newPeriodicBoundaryType--;
+	  mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
+	    new VertexVector(feSpace->getAdmin(), "");      
+	}
+      }
+
+
+      // === In debug mode we make some tests, if the periodic structures are set ===
+      // === in a symmetric way, i.e., if A -> B for a specific boundary type,    ===
+      // === there must be a mapping B -> A with the same boundary type.          ===
+
+#if (DEBUG != 0)
+      for (std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType>::iterator it = periodicVertices.begin();
+	   it != periodicVertices.end(); ++it) {
+	std::pair<DegreeOfFreedom, DegreeOfFreedom> testVertex = 
+	  std::make_pair(it->first.second, it->first.first);
+	
+	TEST_EXIT_DBG(periodicVertices.count(testVertex) == 1)("Should not happen!\n");
+	TEST_EXIT_DBG(periodicVertices[testVertex] == it->second)("Should not happen!\n");
+      }
+
+      for (std::map<std::pair<DofEdge, DofEdge>, BoundaryType>::iterator it = periodicEdges.begin();
+	   it != periodicEdges.end(); ++it) {
+	std::pair<DofEdge, DofEdge> testEdge = 
+	  std::make_pair(it->first.second, it->first.first);
+	
+	TEST_EXIT_DBG(periodicEdges.count(testEdge) == 1)("Should not happen!\n");
+	TEST_EXIT_DBG(periodicEdges[testEdge] == it->second)("Should not happen!\n");
       }
+
+      for (std::map<std::pair<DofFace, DofFace>, BoundaryType>::iterator it = periodicFaces.begin();
+	   it != periodicFaces.end(); ++it) {
+	std::pair<DofFace, DofFace> testFace = 
+	  std::make_pair(it->first.second, it->first.first);
+	
+	TEST_EXIT_DBG(periodicFaces.count(testFace) == 1)("Should not happen!\n");
+	TEST_EXIT_DBG(periodicFaces[testFace] == it->second)("Should not happen!\n");
       }
+#endif
 
 
     }
@@ -1488,7 +1523,7 @@ namespace AMDiS {
     FUNCNAME("MeshDistributor::createBoundaryData()");
 
 
-    // === Clear all relevant data structures, ===
+    // === Clear all relevant data structures. ===
 
     myIntBoundary.clear();
     otherIntBoundary.clear();
@@ -1578,12 +1613,19 @@ namespace AMDiS {
     
     // === Create periodic boundary data structure. ===
 
-    for (std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType>::iterator it = periodicDofs.begin();
-	 it != periodicDofs.end(); ++it) {
+    for (std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType>::iterator it = periodicVertices.begin();
+	 it != periodicVertices.end(); ++it) {
       if (elObjects.isInRank(it->first.first, mpiRank) == false)
 	continue;
 
+      WorldVector<double> c0, c1;
+      mesh->getDofIndexCoords(it->first.first, feSpace, c0);
+      mesh->getDofIndexCoords(it->first.second, feSpace, c1);
+
+      MSG("CREATE BOUNDARY FOR DOF MAP: %d (%.3f %.3f %.3f)<-> %d (%.3f %.3f %.3f)\n", it->first.first, 
+	  c0[0], c0[1], c0[2], it->first.second, c1[0], c1[1], c1[2]);
       ElementObjectData& perDofEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
+      MSG("DATA: %d %d %d\n", perDofEl0.elIndex, VERTEX, perDofEl0.ithObject);
 
       for (std::map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
 	   elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
@@ -2033,6 +2075,8 @@ namespace AMDiS {
     for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
 	 it != periodicBoundary.boundary.end(); ++it) {
 
+      MSG("------------- WITH RANK %d ------------------\n", it->first);
+
       if (it->first == mpiRank) {
 	TEST_EXIT_DBG(it->second.size() % 2 == 0)("Should not happen!\n");
 
@@ -2044,6 +2088,10 @@ namespace AMDiS {
 	  bound.neighObj.el->getVertexDofs(feSpace, bound.neighObj, dofs1);
 	  bound.neighObj.el->getNonVertexDofs(feSpace, bound.neighObj, dofs1);
 
+	  MSG("BOUND-I %d-%d-%d WITH %d-%d-%d\n",
+	      bound.rankObj.elIndex, bound.rankObj.subObj, bound.rankObj.ithObj,
+	      bound.neighObj.elIndex, bound.neighObj.subObj, bound.neighObj.ithObj);
+
 	  TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n");
 
 	  BoundaryType type = bound.type;
@@ -2055,6 +2103,8 @@ namespace AMDiS {
 	    if (periodicDofAssociations[globalDof0].count(type) == 0) {
 	      periodicDof[type][globalDof0] = globalDof1;
 	      periodicDofAssociations[globalDof0].insert(type);
+
+	      MSG("SET(A TYPE %d) DOF %d -> %d\n", type, globalDof0, globalDof1);
 	    }
 	  }
 	}
@@ -2067,12 +2117,20 @@ namespace AMDiS {
 
 	  int nDofs = dofs.size();
 
+	  MSG("BOUND-R (T = %d) %d-%d-%d WITH %d-%d-%d\n",
+	      boundIt->type,
+	      boundIt->rankObj.elIndex, boundIt->rankObj.subObj, boundIt->rankObj.ithObj,
+	      boundIt->neighObj.elIndex, boundIt->neighObj.subObj, boundIt->neighObj.ithObj);
+
+	  
 	  boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
 	  boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
 
-	  for (unsigned int i = 0; i < (dofs.size() - nDofs); i++)
+	  for (unsigned int i = 0; i < (dofs.size() - nDofs); i++) {
+	    MSG("   i = %d   DOF = %d\n", nDofs + i, mapLocalGlobalDofs[*(dofs[nDofs + i])]);
 	    rankToDofType[it->first].push_back(boundIt->type);
 	  }
+	}
 	
 	// Send the global indices to the rank on the other side.
 	stdMpi.getSendData(it->first).reserve(dofs.size());
@@ -2088,6 +2146,8 @@ namespace AMDiS {
     stdMpi.startCommunication<int>(MPI_INT);
 
 
+    MSG("===============================\n");
+
     // === The rank has received the dofs from the rank on the other side of ===
     // === the boundary. Now it can use them to create the mapping between   ===
     // === the periodic dofs in this rank and the corresponding periodic     ===
@@ -2096,6 +2156,9 @@ namespace AMDiS {
 
     for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
 	 it != periodicBoundary.boundary.end(); ++it) {
+
+      MSG("------------- WITH RANK %d ------------------\n", it->first);
+
       DofContainer& dofs = rankPeriodicDofs[it->first];
       std::vector<int>& types = rankToDofType[it->first];
 
@@ -2110,23 +2173,28 @@ namespace AMDiS {
 	// Check if this global dof with the corresponding boundary type was
 	// not added before by another periodic boundary from other rank.
 	if (periodicDofAssociations[globalDofIndex].count(type) == 0) {
+
+	  MSG("SET(B-%d TYPE %d) DOF %d -> %d\n", i, type, globalDofIndex, mapGlobalDofIndex);
+
+
 	  periodicDof[type][globalDofIndex] = mapGlobalDofIndex;
 	  periodicDofAssociations[globalDofIndex].insert(type);
+	} else {
+	  MSG("ASSOC ALREADY SET FOR %d TYPE %d\n", i, type);
 	}
       }
     }
 
-
 #if (DEBUG != 0)
-    // In 2D, a periodic DOF can have either 1 or 3 associations. Check this!
-    if (mesh->getDim() == 2) {
     for (std::map<int, std::set<BoundaryType> >::iterator it = periodicDofAssociations.begin();
 	 it != periodicDofAssociations.end(); ++it) {
+      WorldVector<double> c;
+      mesh->getDofIndexCoords(it->first, feSpace, c);
       int nAssoc = it->second.size();
-	TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3)
-	  ("Should not happen! DOF %d has %d periodic associations!\n", 
-	   it->first, nAssoc);
-      }
+      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || 
+		    (mesh->getDim() == 3 && (nAssoc == 7 || nAssoc == 11)))
+	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
+	 it->first, c[0], c[1], (mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
     }    
 #endif
   }
@@ -2179,7 +2247,7 @@ namespace AMDiS {
 
     elObjects.serialize(out);
 
-    SerUtil::serialize(out, periodicDofs);
+    SerUtil::serialize(out, periodicVertices);
     SerUtil::serialize(out, periodicEdges);
     SerUtil::serialize(out, periodicFaces);
 
@@ -2243,7 +2311,7 @@ namespace AMDiS {
 
     elObjects.deserialize(in);
    
-    SerUtil::deserialize(in, periodicDofs);
+    SerUtil::deserialize(in, periodicVertices);
     SerUtil::deserialize(in, periodicEdges);
     SerUtil::deserialize(in, periodicFaces);
    

AMDiS/src/parallel/MeshDistributor.h

@@ -479,7 +479,7 @@ namespace AMDiS {
     std::map<int, int> macroElIndexTypeMap;
 
     // The following three data structures store periodic DOFs, edges and faces.
-    std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicDofs;
+    std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicVertices;
     std::map<std::pair<DofEdge, DofEdge>, BoundaryType> periodicEdges;
     std::map<std::pair<DofFace, DofFace>, BoundaryType> periodicFaces;
 

AMDiS/src/parallel/ParallelDebug.cc

@@ -139,6 +139,89 @@ namespace AMDiS {
   }
 
 
+  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb)
+  {
+    FUNCNAME("ParallelDebug::testPeriodicBoundary()");
+
+    typedef MeshDistributor::PeriodicDofMap PeriodicDofMap;
+
+    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);
+
+    if (pdb.mpiRank == 0) {
+      for (int i = 1; i < pdb.mpiSize; i++)
+	stdMpi.recv(i);
+    } else {
+      stdMpi.send(0, pdb.periodicDof);
+    }
+
+    stdMpi.startCommunication<int>(MPI_INT);
+
+    int foundError = 0;
+
+    // === The boundary DOFs are checked only on the zero rank. === 
+
+    if (pdb.mpiRank == 0) {
+      std::map<int, PeriodicDofMap> rankToMaps;
+      PeriodicDofMap dofMap = pdb.periodicDof;
+      rankToMaps[0] = dofMap;
+
+      for (int i = 1; i < pdb.mpiSize; i++) {
+	PeriodicDofMap &otherMap = stdMpi.getRecvData(i);
+	rankToMaps[i] = otherMap;
+	
+	for (PeriodicDofMap::iterator it = otherMap.begin(); 
+	     it != otherMap.end(); ++it) {
+	  for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
+	       dofIt != it->second.end(); ++dofIt) {
+	    if (dofMap.count(it->first) == 1 &&
+		dofMap[it->first].count(dofIt->first) == 1) {
+	      TEST_EXIT_DBG(dofMap[it->first][dofIt->first] == dofIt->second)
+		("Should not happen!\n");
+	    } else {
+	      dofMap[it->first][dofIt->first] = dofIt->second;
+	    }
+	  }
+	}
+      }
+
+
+      // === Now we test if global DOF A is mapped to B, then B must be mapped ===
+      // === to A for the same boundary type.                                  ===
+
+      for (PeriodicDofMap::iterator it = dofMap.begin(); 
+	   it != dofMap.end(); ++it) {
+	for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
+	     dofIt != it->second.end(); ++dofIt) {
+	  if (it->second[dofIt->second] != dofIt->first) {
+	    MSG("[DBG]  For boundary type %d: DOF %d -> %d, but %d -> %d!\n",
+		it ->first, 
+		dofIt->first, dofIt->second, 
+		dofIt->second, it->second[dofIt->second]);
+
+	    for (int i = 0; i < pdb.mpiSize; i++) {
+	      if (rankToMaps[i][it->first].count(dofIt->first) == 1) {
+		MSG("[DBG]    %d -> %d in rank %d\n", 
+		    dofIt->first, rankToMaps[i][it->first][dofIt->first], i);
+	      }
+
+	      if (rankToMaps[i][it->first].count(dofIt->second) == 1) {
+		MSG("[DBG]    %d -> %d in rank %d\n", 
+		    dofIt->second, rankToMaps[i][it->first][dofIt->second], i);
+	      }
+	    }
+	    
+	    ERROR("Wrong periodic DOFs!\n");
+	    foundError = 1;
+	  }
+	}
+      }
+    }
+
+    mpi::globalAdd(foundError);
+    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
+  }
+
+
   void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
   {
     FUNCNAME("ParallelDebug::testCommonDofs()");
@@ -573,7 +656,8 @@ namespace AMDiS {
     }
 
     for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
-      MSG("Periodic boundary with rank %d: \n", it.getRank());
+      MSG("Periodic boundary (ID %d) with rank %d: \n", 
+	  it->type, it.getRank());
       MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
 	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
       MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",