Work on bla.

AMDiS/src/parallel/ElementObjectDatabase.cc

@@ -41,9 +41,12 @@ namespace AMDiS {
 
       elInfo = stack.traverseNext(elInfo);
     }
-
-    // Create periodic data, if there are periodic boundary conditions.
-    createPeriodicData();
+    
+    // Handle periodic boudaries
+    if (removePeriodicBoundary == false) 
+      createPeriodicData();
+    else
+      removePeriodicData();
 
     // Create data about the reverse modes of neighbouring elements.
     createReverseModeData();
@@ -99,40 +102,21 @@ namespace AMDiS {
 	  DofEdge edge1 = neigh->getEdge(elInfo->getOppVertex(i));
 	  BoundaryType boundaryType = elInfo->getBoundary(EDGE, i);
 
-
-	  if (removePeriodicBoundary) {
-	    ElementObjectData elObjEdge(el->getIndex(), i);
-	    edgeElements[edge1].push_back(elObjEdge);
-	    edgeLocalMap[elObjEdge] = edge1;
-
-#if 0
-	    ElementObjectData elObjVertex0(el->getIndex(), 
-					   el->getVertexOfEdge(i, 0));
-	    vertexElements[edge1.first].push_back(elObjVertex0);
-	    vertexLocalMap[elObjVertex0] = edge1.first;
-
-	    ElementObjectData elObjVertex1(el->getIndex(), 
-					   el->getVertexOfEdge(i, 1));
-	    vertexElements[edge1.second].push_back(elObjVertex1);
-	    vertexLocalMap[elObjVertex1] = edge1.first;
-#endif
-	  } else {  
-	    // Add the periodic edge.
-	    periodicEdges[make_pair(edge0, edge1)] = boundaryType;
-	    periodicEdgeAssoc[edge0].insert(edge1);
-	    
-	    // Add both vertices of the edge to be periodic.
-	    periodicVertices[make_pair(edge0.first, edge1.first)] = boundaryType;
-	    periodicVertices[make_pair(edge0.second, edge1.second)] = boundaryType;
-	    periodicDofAssoc[edge0.first].insert(boundaryType);
-	    periodicDofAssoc[edge0.second].insert(boundaryType);
-	    
-	    TEST_EXIT_DBG(edge0.first == 
-			  mesh->getPeriodicAssociations(boundaryType)[edge1.first] &&
-			  edge0.second == 
-			  mesh->getPeriodicAssociations(boundaryType)[edge1.second])
-	      ("Should not happen!\n");
-	  }
+	  // Add the periodic edge.
+	  periodicEdges[make_pair(edge0, edge1)] = boundaryType;
+	  periodicEdgeAssoc[edge0].insert(edge1);
+	  
+	  // Add both vertices of the edge to be periodic.
+	  periodicVertices[make_pair(edge0.first, edge1.first)] = boundaryType;
+	  periodicVertices[make_pair(edge0.second, edge1.second)] = boundaryType;
+	  periodicDofAssoc[edge0.first].insert(boundaryType);
+	  periodicDofAssoc[edge0.second].insert(boundaryType);
+	  
+	  TEST_EXIT_DBG(edge0.first == 
+			mesh->getPeriodicAssociations(boundaryType)[edge1.first] &&
+			edge0.second == 
+			mesh->getPeriodicAssociations(boundaryType)[edge1.second])
+	    ("Should not happen!\n");
 	}
       }
       break;
@@ -148,9 +132,6 @@ namespace AMDiS {
 	  DofFace face1 = neigh->getFace(elInfo->getOppVertex(i));
 	  BoundaryType boundaryType = elInfo->getBoundary(FACE, i);
 
-	  TEST_EXIT(removePeriodicBoundary == false)
-	    ("Not yet implemented for 3D!\n");
-	  
 	  // Add the periodic face.
 	  periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i);
 	  
@@ -197,6 +178,46 @@ namespace AMDiS {
     }       
   }
 
+  void ElementObjectDatabase::addVertex(Element *el, int ith)
+  {
+    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;
+  }
+
+
+  void ElementObjectDatabase::addEdge(Element *el, int ith)
+  {
+    FUNCNAME("ElementObjectDatabase::addEdge()");
+
+    DofEdge edge = el->getEdge(ith);
+    int elIndex = el->getIndex();
+    ElementObjectData elObj(elIndex, ith);
+    
+    edgeElements[edge].push_back(elObj);
+    edgeLocalMap[elObj] = edge;
+  }
+
+  
+  void ElementObjectDatabase::addFace(Element *el, int ith)
+  {
+    DofFace face = el->getFace(ith);
+    int elIndex = el->getIndex();
+    ElementObjectData elObj(elIndex, ith);
+    
+    faceElements[face].push_back(elObj);
+    faceLocalMap[elObj] = face;
+  }
+
 
   void ElementObjectDatabase::createPeriodicData()
   {
@@ -391,6 +412,10 @@ namespace AMDiS {
   }
 
 
+void ElementObjectDatabase::removePeriodicData()
+{
+}
+
   BoundaryType ElementObjectDatabase::getNewBoundaryType()
   {
     FUNCNAME("ElementObjectDatabase::getNewBoundaryType()");

AMDiS/src/parallel/ElementObjectDatabase.h

@@ -51,7 +51,8 @@ namespace AMDiS {
   struct ElementObjectData {
     ElementObjectData(int a = -1, int b = 0)
       : elIndex(a),
-	ithObject(b)
+	ithObject(b),
+	mappedOnPeriodicBoundary(false)
     {}
 
     /// Index of the element this object is part of.
@@ -59,12 +60,17 @@ 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.
@@ -72,6 +78,7 @@ namespace AMDiS {
     {
       SerUtil::deserialize(in, elIndex);
       SerUtil::deserialize(in, ithObject);
+      SerUtil::deserialize(in, mappedOnPeriodicBoundary);
     }
 
     /// Compare this element object with another one.
@@ -426,37 +433,13 @@ namespace AMDiS {
     void addElement(ElInfo *elInfo);
 
     /// Adds the i-th DOF vertex of an element to the object database.
-    void addVertex(Element *el, int ith)
-    {
-      DegreeOfFreedom vertex = el->getDof(ith, 0);
-      int elIndex = el->getIndex();
-      ElementObjectData elObj(elIndex, ith);
-
-      vertexElements[vertex].push_back(elObj);
-      vertexLocalMap[elObj] = vertex;
-    }
+    void addVertex(Element *el, int ith);
 
     /// Adds the i-th edge of an element to the object database.
-    void addEdge(Element *el, int ith)
-    {
-      DofEdge edge = el->getEdge(ith);
-      int elIndex = el->getIndex();
-      ElementObjectData elObj(elIndex, ith);
-
-      edgeElements[edge].push_back(elObj);
-      edgeLocalMap[elObj] = edge;
-    }
+    void addEdge(Element *el, int ith);
 
     /// Adds the i-th face of an element to the object database.
-    void addFace(Element *el, int ith)
-    {
-      DofFace face = el->getFace(ith);
-      int elIndex = el->getIndex();
-      ElementObjectData elObj(elIndex, ith);
-
-      faceElements[face].push_back(elObj);
-      faceLocalMap[elObj] = face;
-    }
+    void addFace(Element *el, int ith);
 
     /** \brief
      * Creates final data of the periodic boundaries. Must be called after all
@@ -468,7 +451,8 @@ namespace AMDiS {
      * indirectly connected.
      */
     void createPeriodicData();
-
+    
+    void removePeriodicData();
 
     /// Creates on all boundaries the reverse mode flag.
     void createReverseModeData();

AMDiS/src/parallel/InteriorBoundary.cc

@@ -48,6 +48,8 @@ namespace AMDiS {
     for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
       GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
 
+      MSG("---> start with geo index = %d\n", geoIndex);
+
       while (elObjDb.iterate(geoIndex)) {
 	map<int, ElementObjectData>& objData = elObjDb.getIterateData();
 
@@ -75,6 +77,9 @@ 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);
@@ -82,7 +87,8 @@ namespace AMDiS {
 	bound.rankObj.elType = elObjDb.getElementType(rankBoundEl.elIndex);
 	bound.rankObj.subObj = geoIndex;
 	bound.rankObj.ithObj = rankBoundEl.ithObject;
-	
+
+
 	if (geoIndex == FACE) {
 	  for (int edgeNo = 0; edgeNo < 3; edgeNo++) {
 	    int edgeOfFace = 
@@ -92,7 +98,7 @@ namespace AMDiS {
 	  }
 	}
 	
-	
+
 	if (owner == globalMpiRank) {
 	  for (map<int, ElementObjectData>::iterator it2 = objData.begin();
 	       it2 != objData.end(); ++it2) {
@@ -125,6 +131,12 @@ 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;

AMDiS/src/parallel/MeshDistributor.cc

@@ -1558,15 +1558,11 @@ namespace AMDiS {
     elObjDb.updateRankData();
 
     intBoundary.create(levelData, 0, elObjDb);
-#if (DEBUG != 0)
     ParallelDebug::printBoundaryInfo(intBoundary);
-#endif
 
     if (levelData.getLevelNumber() > 1) {
       intBoundarySd.create(levelData, 1, elObjDb);
-#if (DEBUG != 0)
       ParallelDebug::printBoundaryInfo(intBoundarySd, 0, true);
-#endif
     }
   }
 
@@ -1692,6 +1688,10 @@ namespace AMDiS {
   
     lastMeshChangeIndex = mesh->getChangeIndex();
 
+    MSG("--------> REMOVE THIS LINE <------------\n");
+    debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex-" + 
+				 lexical_cast<string>(mpiRank) + ".vtu");
+
 
 #if (DEBUG != 0)
     ParallelDebug::testDofContainerCommunication(*this);

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -436,21 +436,23 @@ namespace AMDiS {
     // === Define all vertices on the interior boundaries of the macro mesh ===
     // === to be primal variables.                                          ===
 
-    /// Set of DOF indices that are considered to be primal variables.
-    
+    // Set of DOF indices that are considered to be primal variables.
     DofContainerSet& vertices = 
       meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
 
     DofIndexSet primals;
     for (DofContainerSet::iterator it = vertices.begin(); 
 	 it != vertices.end(); ++it) {
-      WorldVector<double> c;
-      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
-
       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;
+	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
+
 	if (fabs(c[0]) < e || fabs(c[1]) < e ||
 	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
 	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {