Fixed periodic boundaries for FETI-DP

AMDiS/src/parallel/ElementObjectDatabase.cc

@@ -43,10 +43,7 @@ namespace AMDiS {
     }
     
     // Handle periodic boudaries
-    if (removePeriodicBoundary == false) 
-      createPeriodicData();
-    else
-      removePeriodicData();
+    createPeriodicData();
 
     // Create data about the reverse modes of neighbouring elements.
     createReverseModeData();
@@ -183,13 +180,7 @@ namespace AMDiS {
     DegreeOfFreedom vertex = el->getDof(ith, 0);
     int elIndex = el->getIndex();
     ElementObjectData elObj(elIndex, ith);
-
-    if (elIndex == 53 && ith == 0)
-      MSG("A: 53/0 ON DOF %d\n", vertex);
-
-    if (elIndex == 229 && ith == 0)
-      MSG("A: 229/0 ON DOF %d\n", vertex);
-    
+   
     vertexElements[vertex].push_back(elObj);
     vertexLocalMap[elObj] = vertex;
   }
@@ -412,10 +403,6 @@ namespace AMDiS {
   }
 
 
-void ElementObjectDatabase::removePeriodicData()
-{
-}
-
   BoundaryType ElementObjectDatabase::getNewBoundaryType()
   {
     FUNCNAME("ElementObjectDatabase::getNewBoundaryType()");
@@ -626,26 +613,51 @@ void ElementObjectDatabase::removePeriodicData()
 
     TEST_EXIT_DBG(macroElementRankMap)("Should not happen!\n");
 
-    int owner = -1;
-    vector<ElementObjectData> *objData;
-    
     switch (iterGeoPos) {
     case VERTEX:
-      objData = &(vertexElements[vertexIter->first]);
+      return getOwner(vertexElements[vertexIter->first], level);
       break;
     case EDGE:
-      objData = &(edgeElements[edgeIter->first]);
+      return getOwner(edgeElements[edgeIter->first], level);
       break;
     case FACE:
-      objData = &(faceElements[faceIter->first]);
+      return getOwner(faceElements[faceIter->first], level);
       break;
     }
     
+    ERROR_EXIT("There is something reallllly wrong!\n");
+    return -1;
+  }
+
+
+  int ElementObjectDatabase::getOwner(DegreeOfFreedom vertex, int level)
+  {
+    return getOwner(vertexElements[vertex], level);
+  }
+
+
+  int ElementObjectDatabase::getOwner(DofEdge edge, int level)
+  {
+    return getOwner(edgeElements[edge], level);
+  }
+
+
+  int ElementObjectDatabase::getOwner(DofFace face, int level)
+  {
+    return getOwner(faceElements[face], level);
+  }
+
+
+  int ElementObjectDatabase::getOwner(vector<ElementObjectData>& objData, 
+				      int level)
+  {
+    int owner = -1;
+
     std::set<int> &levelRanks = levelData->getLevelRanks(level);
     bool allRanks = (levelRanks.size() == 1 && *(levelRanks.begin()) == -1);
     
-    for (vector<ElementObjectData>::iterator it = objData->begin();
-	 it != objData->end(); ++it) {
+    for (vector<ElementObjectData>::iterator it = objData.begin();
+	 it != objData.end(); ++it) {
       int elRank = (*macroElementRankMap)[it->elIndex];
       if (allRanks || levelRanks.count(elRank))
 	owner = std::max(owner, elRank);

AMDiS/src/parallel/ElementObjectDatabase.h

@@ -51,8 +51,7 @@ namespace AMDiS {
   struct ElementObjectData {
     ElementObjectData(int a = -1, int b = 0)
       : elIndex(a),
-	ithObject(b),
-	mappedOnPeriodicBoundary(false)
+	ithObject(b)
     {}
 
     /// Index of the element this object is part of.
@@ -61,16 +60,11 @@ namespace AMDiS {
     /// Index of the object within the element.
     int ithObject;
 
-    /// If true, the element does not exists in mesh but is due to a mapping
-    /// on a periodic boundary.
-    bool mappedOnPeriodicBoundary;
-    
     /// Write this element object to disk.
     void serialize(ostream &out) const
     {
       SerUtil::serialize(out, elIndex);
       SerUtil::serialize(out, ithObject);
-      SerUtil::serialize(out, mappedOnPeriodicBoundary);
     }
 
     /// Read this element object from disk.
@@ -78,7 +72,6 @@ namespace AMDiS {
     {
       SerUtil::deserialize(in, elIndex);
       SerUtil::deserialize(in, ithObject);
-      SerUtil::deserialize(in, mappedOnPeriodicBoundary);
     }
 
     /// Compare this element object with another one.
@@ -114,11 +107,8 @@ namespace AMDiS {
 	mesh(NULL),
 	iterGeoPos(CENTER),
 	macroElementRankMap(NULL),
-	levelData(NULL),
-	removePeriodicBoundary(false)
-    {
-      Parameters::get("parallel->remove periodic boundary", removePeriodicBoundary);
-    }
+	levelData(NULL)
+    {}
 
     void setFeSpace(const FiniteElemSpace *fe)
     {
@@ -242,6 +232,15 @@ namespace AMDiS {
     /// Returns the rank owner of the current iterator position.
     int getIterateOwner(int level);
 
+    /// Returns the owner of a macro element vertex.
+    int getOwner(DegreeOfFreedom vertex, int level);
+
+    /// Returns the owner of a macro element edge.
+    int getOwner(DofEdge edge, int level);
+	
+    /// Returns the owner of a macro element face.
+    int getOwner(DofFace face, int level);
+
     /// Returns the rank owner of the current iterator position.
     int getIterateMaxLevel();
 
@@ -452,8 +451,6 @@ namespace AMDiS {
      */
     void createPeriodicData();
     
-    void removePeriodicData();
-
     /// Creates on all boundaries the reverse mode flag.
     void createReverseModeData();
 
@@ -474,6 +471,7 @@ namespace AMDiS {
     /// Some auxiliary function to read the element object database from disk.
     void deserialize(istream &in, map<int, ElementObjectData>& data);
 
+    int getOwner(vector<ElementObjectData>& objData, int level);
   private:
     const FiniteElemSpace* feSpace;
 
@@ -562,12 +560,6 @@ namespace AMDiS {
     map<int, int> macroElIndexTypeMap;
 
     MeshLevelData* levelData;
-
-    /// If this variable is set to true, the mesh distributor removes all 
-    /// periodic boundary conditions. The element neighbourhood relation is
-    /// not changed. Thus, when using some domain decomposition method, this is
-    /// a natural way to deal with periodic boundary conditions.
-    bool removePeriodicBoundary;
   };
 
 }

AMDiS/src/parallel/InteriorBoundary.cc

@@ -77,9 +77,6 @@ namespace AMDiS {
 	int owner = elObjDb.getIterateOwner(level);
 	ElementObjectData& rankBoundEl = objData[globalMpiRank];
 	
-	if (rankBoundEl.mappedOnPeriodicBoundary)
-	  continue;
-
 	AtomicBoundary bound;
 	bound.maxLevel = elObjDb.getIterateMaxLevel();
 	bound.rankObj.el = elObjDb.getElementPtr(rankBoundEl.elIndex);
@@ -131,13 +128,7 @@ namespace AMDiS {
 	    ("Should not happen!\n");
 	  
 	  ElementObjectData& ownerBoundEl = objData[owner];
-
-	  if (rankBoundEl.elIndex == 53 &&
-	      geoIndex == 1 &&
-	      rankBoundEl.ithObject == 0) {
-	    MSG("OWNER: %d %d\n", owner, ownerBoundEl);
-	  }
-	  
+  
 	  bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
 	  bound.neighObj.elIndex = ownerBoundEl.elIndex;
 	  bound.neighObj.elType = -1;
@@ -161,6 +152,10 @@ namespace AMDiS {
 
     // === Create periodic boundary data structure. ===
 
+    int removePeriodicBoundary = 0;
+    Parameters::get("parallel->remove periodic boundary", removePeriodicBoundary);
+
+
     for (PerBoundMap<DegreeOfFreedom>::iterator it = elObjDb.getPeriodicVertices().begin();
 	 it != elObjDb.getPeriodicVertices().end(); ++it) {
       if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
@@ -188,19 +183,57 @@ namespace AMDiS {
 	bound.neighObj.subObj = VERTEX;
 	bound.neighObj.ithObj = perDofEl1.ithObject;
 
-	bound.type = it->second;
-
-	if (MeshDistributor::sebastianMode) {
-	  if (bound.rankObj.elIndex == 3 &&
-	      bound.rankObj.ithObj == 1 ||
-	      bound.rankObj.elIndex == 78 &&
-	      bound.rankObj.ithObj == 2) {
-	    AtomicBoundary& b = getNewPeriodic(otherElementRank);
-	    b = bound;	    	    
+	if (removePeriodicBoundary) {
+	  bound.type = INTERIOR;
+	  
+	  map<int, ElementObjectData> objData = 
+	    elObjDb.getElementsInRank(it->first.first);
+	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
+			 elObjDb.getElementsInRank(it->first.second).end());
+	  
+	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
+			       elObjDb.getOwner(it->first.second, level));
+	  if (owner == globalMpiRank) {
+	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
+	     	 it2 != objData.end(); ++it2) {
+	      if (it2->first == globalMpiRank)
+	     	continue;
+	      
+	      if (!allRanks && levelRanks.count(it2->first) == 0)
+	     	continue;
+	      
+	      AtomicBoundary& b = getNewOwn(it2->first);	  
+	      b = bound;	  
+	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
+	      b.neighObj.elIndex = it2->second.elIndex;
+	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
+	      b.neighObj.ithObj = it2->second.ithObject;
+	    }
+	  } else {
+	    AtomicBoundary& b = getNewOther(owner);
+	    b = bound;
+	    
+	    ElementObjectData& ownerBoundEl = objData[owner];    
+	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
+	    b.neighObj.elIndex = ownerBoundEl.elIndex;
+	    b.neighObj.elType = -1;
+	    b.neighObj.ithObj = ownerBoundEl.ithObject;
 	  }
 	} else {
-	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
-	  b = bound;	    	  
+	  bound.type = it->second;
+	  
+	  if (MeshDistributor::sebastianMode) {
+	    if (bound.rankObj.elIndex == 3 &&
+		bound.rankObj.ithObj == 1 ||
+		bound.rankObj.elIndex == 78 &&
+		bound.rankObj.ithObj == 2) {
+	      AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	      b = bound;	    	    
+	    }
+	  } else {
+	    AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	    b = bound;	    	  
+	  }
 	}
       }
     }
@@ -232,18 +265,62 @@ namespace AMDiS {
 	bound.neighObj.elType = elObjDb.getElementType(perEdgeEl1.elIndex);
 	bound.neighObj.subObj = EDGE;
 	bound.neighObj.ithObj = perEdgeEl1.ithObject;
-	
-	bound.type = it->second;
 
-	AtomicBoundary& b = getNewPeriodic(otherElementRank);
-	b = bound;
-     
-	if (globalMpiRank > otherElementRank)
-	  b.neighObj.reverseMode = 
-	    elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
-	else
-	  b.rankObj.reverseMode = 
-	    elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
+	if (removePeriodicBoundary) {
+	  bound.type = INTERIOR;
+
+	  map<int, ElementObjectData> objData = 
+	    elObjDb.getElementsInRank(it->first.first);
+	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
+			 elObjDb.getElementsInRank(it->first.second).end());
+	  
+	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
+			       elObjDb.getOwner(it->first.second, level));
+	  ElementObjectData& rankBoundEl = objData[globalMpiRank];
+
+	  if (owner == globalMpiRank) {
+	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
+	     	 it2 != objData.end(); ++it2) {
+	      if (it2->first == globalMpiRank)
+	     	continue;
+	      
+	      if (!allRanks && levelRanks.count(it2->first) == 0)
+	     	continue;
+	      
+	      AtomicBoundary& b = getNewOwn(it2->first);	  
+	      b = bound;	  
+	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
+	      b.neighObj.elIndex = it2->second.elIndex;
+	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
+	      b.neighObj.ithObj = it2->second.ithObject;
+	      b.neighObj.reverseMode = 
+		elObjDb.getEdgeReverseMode(rankBoundEl, it2->second);
+	    }
+	  } else {
+	    AtomicBoundary& b = getNewOther(owner);
+	    b = bound;
+	    
+	    ElementObjectData& ownerBoundEl = objData[owner];    
+	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
+	    b.neighObj.elIndex = ownerBoundEl.elIndex;
+	    b.neighObj.elType = -1;
+	    b.neighObj.ithObj = ownerBoundEl.ithObject;
+	    b.rankObj.reverseMode = 
+	      elObjDb.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
+	  }
+	} else {
+	  bound.type = it->second;
+	  
+	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	  b = bound;
+	  
+	  if (globalMpiRank > otherElementRank)
+	    b.neighObj.reverseMode = 
+	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
+	  else
+	    b.rankObj.reverseMode = 
+	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
+	}
       }
     }
 
@@ -280,17 +357,61 @@ namespace AMDiS {
 	bound.neighObj.subObj = FACE;
 	bound.neighObj.ithObj = perFaceEl1.ithObject;
 	
-	bound.type = it->second;
-	
-	AtomicBoundary& b = getNewPeriodic(otherElementRank);
-	b = bound;
-     
-	if (globalMpiRank > otherElementRank)
-	  b.neighObj.reverseMode = 
-	    elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
-	else
-	  b.rankObj.reverseMode = 
-	    elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
+	if (removePeriodicBoundary) {
+	  bound.type = INTERIOR;
+
+	  map<int, ElementObjectData> objData = 
+	    elObjDb.getElementsInRank(it->first.first);
+	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
+			 elObjDb.getElementsInRank(it->first.second).end());
+	  
+	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
+			       elObjDb.getOwner(it->first.second, level));
+	  ElementObjectData& rankBoundEl = objData[globalMpiRank];
+
+	  if (owner == globalMpiRank) {
+	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
+	     	 it2 != objData.end(); ++it2) {
+	      if (it2->first == globalMpiRank)
+	     	continue;
+	      
+	      if (!allRanks && levelRanks.count(it2->first) == 0)
+	     	continue;
+	      
+	      AtomicBoundary& b = getNewOwn(it2->first);	  
+	      b = bound;	  
+	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
+	      b.neighObj.elIndex = it2->second.elIndex;
+	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
+	      b.neighObj.ithObj = it2->second.ithObject;
+	      b.neighObj.reverseMode = 
+		elObjDb.getFaceReverseMode(rankBoundEl, it2->second);
+	    }
+	  } else {
+	    AtomicBoundary& b = getNewOther(owner);
+	    b = bound;
+	    
+	    ElementObjectData& ownerBoundEl = objData[owner];    
+	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
+	    b.neighObj.elIndex = ownerBoundEl.elIndex;
+	    b.neighObj.elType = -1;
+	    b.neighObj.ithObj = ownerBoundEl.ithObject;
+	    b.rankObj.reverseMode = 
+	      elObjDb.getFaceReverseMode(rankBoundEl, ownerBoundEl);
+	  }
+	} else {
+	  bound.type = it->second;
+	  
+	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	  b = bound;
+	  
+	  if (globalMpiRank > otherElementRank)
+	    b.neighObj.reverseMode = 
+	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
+	  else
+	    b.rankObj.reverseMode = 
+	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
+	}
       }
     }
 
@@ -358,6 +479,8 @@ namespace AMDiS {
     }
 
 
+#if 0
+
     // === Do the same for the periodic boundaries. ===
 
     if (periodic.size() > 0) {
@@ -412,6 +535,8 @@ namespace AMDiS {
 	}
       }     
     } // periodicBoundary.boundary.size() > 0
+
+#endif
   }
 
 

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -445,9 +445,6 @@ namespace AMDiS {
 	 it != vertices.end(); ++it) {
       double e = 1e-8;
       if (meshLevel == 0) {
-	WorldVector<double> c;
-	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
-	MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
 	primals.insert(**it);
       } else {
 	WorldVector<double> c;