New approach for deleting double DOFs in mesh repartitioning.

AMDiS/src/Debug.cc

@@ -768,6 +768,30 @@ namespace AMDiS {
     }
 
 
+    void testDofsByCoords(FiniteElemSpace *feSpace,
+			  DofContainer &dofs0, DofContainer &dofs1)
+    {
+      FUNCNAME("debug::testDofsByCoords()");
+
+      TEST_EXIT(dofs0.size() == dofs1.size())
+	("The dof containers have different sizes %d %d!\n",
+	 dofs0.size(), dofs1.size());
+
+      DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
+      feSpace->getMesh()->getDofIndexCoords(feSpace, coords);
+
+      for (unsigned int i = 0; i < dofs0.size(); i++) {
+	WorldVector<double> tmp = coords[*(dofs0[i])];
+	tmp -= coords[*(dofs1[i])];
+
+	TEST_EXIT(norm(tmp) < 1e-13)
+	  ("DOFs %d and %d (i = %d) have different coords!\n", 
+	   *(dofs0[i]), *(dofs1[i]), i);
+      }
+      
+    }
+
+
   } // namespace debug
   
 } // namespace AMDiS

AMDiS/src/Debug.h

@@ -193,6 +193,9 @@ namespace AMDiS {
 		  const DegreeOfFreedom* dof3,
 		  DofContainer &vec);    
 
+    void testDofsByCoords(FiniteElemSpace *feSpace,
+			  DofContainer &dofs0, DofContainer &dofs1);
+
   }
 }
 

AMDiS/src/parallel/ElementObjectData.cc

@@ -97,6 +97,32 @@ namespace AMDiS {
 
 
 
+    nSize = vertexLocalMap.size();
+    SerUtil::serialize(out, nSize);
+    for (std::map<ElementObjectData, DegreeOfFreedom>::iterator it = vertexLocalMap.begin();
+	 it != vertexLocalMap.end(); ++it) {
+      it->first.serialize(out);
+      SerUtil::serialize(out, it->second);
+    }
+
+    nSize = edgeLocalMap.size();
+    SerUtil::serialize(out, nSize);
+    for (std::map<ElementObjectData, DofEdge>::iterator it = edgeLocalMap.begin();
+	 it != edgeLocalMap.end(); ++it) {
+      it->first.serialize(out);
+      SerUtil::serialize(out, it->second);
+    }
+
+    nSize = faceLocalMap.size();
+    SerUtil::serialize(out, nSize);
+    for (std::map<ElementObjectData, DofFace>::iterator it = faceLocalMap.begin();
+	 it != faceLocalMap.end(); ++it) {
+      it->first.serialize(out);
+      SerUtil::serialize(out, it->second);
+    }
+
+
+
     SerUtil::serialize(out, vertexOwner);
     SerUtil::serialize(out, edgeOwner);
     SerUtil::serialize(out, faceOwner);
@@ -165,6 +191,38 @@ namespace AMDiS {
     }
    
 
+
+    SerUtil::deserialize(in, nSize);
+    vertexLocalMap.clear();
+    for (int i = 0; i < nSize; i++) {
+      ElementObjectData data;
+      DegreeOfFreedom dof;
+      data.deserialize(in);
+      SerUtil::deserialize(in, dof);
+      vertexLocalMap[data] = dof;
+    }
+
+    SerUtil::deserialize(in, nSize);
+    edgeLocalMap.clear();
+    for (int i = 0; i < nSize; i++) {
+      ElementObjectData data;
+      DofEdge edge;
+      data.deserialize(in);
+      SerUtil::deserialize(in, edge);
+      edgeLocalMap[data] = edge;
+    }
+
+    SerUtil::deserialize(in, nSize);
+    faceLocalMap.clear();
+    for (int i = 0; i < nSize; i++) {
+      ElementObjectData data;
+      DofFace face;
+      data.deserialize(in);
+      SerUtil::deserialize(in, face);
+      faceLocalMap[data] = face;
+    }
+
+
     SerUtil::deserialize(in, vertexOwner);
     SerUtil::deserialize(in, edgeOwner);
     SerUtil::deserialize(in, faceOwner);

AMDiS/src/parallel/ElementObjectData.h

@@ -48,7 +48,7 @@ namespace AMDiS {
     BoundaryType boundaryType;
 
 
-    void serialize(std::ostream &out)
+    void serialize(std::ostream &out) const
     {
       SerUtil::serialize(out, elIndex);
       SerUtil::serialize(out, ithObject);
@@ -63,6 +63,20 @@ namespace AMDiS {
       SerUtil::deserialize(in, boundaryType);
     }
 
+
+    bool operator==(ElementObjectData& cmp) const
+    {
+      return (elIndex == cmp.elIndex &&
+	      ithObject == cmp.ithObject &&
+	      boundaryType == cmp.boundaryType);
+    }
+
+    bool operator<(const ElementObjectData& rhs) const
+    {
+      return (elIndex < rhs.elIndex ||
+	      (elIndex == rhs.elIndex &&
+	       ithObject < rhs.ithObject));
+    }
   };
 
 
@@ -73,26 +87,56 @@ namespace AMDiS {
       : iterGeoPos(CENTER)
     {}
 
-    void addVertex(DegreeOfFreedom vertex, 
-		   int elIndex, int ith, BoundaryType bound = INTERIOR)
+
+    void addVertex(Element *el, int ith, BoundaryType bound = INTERIOR)
     {
-      vertexElements[vertex].push_back(ElementObjectData(elIndex, ith, bound));
+      DegreeOfFreedom vertex = el->getDof(ith, 0);
+      int elIndex = el->getIndex();
+      ElementObjectData elObj(elIndex, ith, bound);
+
+      vertexElements[vertex].push_back(elObj);
+      vertexLocalMap[elObj] = vertex;
     }
 
-    void addEdge(DofEdge edge, 
-		 int elIndex, int ith, BoundaryType bound = INTERIOR)
+
+    void addEdge(Element *el, int ith, BoundaryType bound = INTERIOR)
     {
-      edgeElements[edge].push_back(ElementObjectData(elIndex, ith, bound));
+      DofEdge edge = el->getEdge(ith);
+      int elIndex = el->getIndex();
+      ElementObjectData elObj(elIndex, ith, bound);
+
+      edgeElements[edge].push_back(elObj);
+      edgeLocalMap[elObj] = edge;
+    }
+
+
+    void addFace(Element *el, int ith, BoundaryType bound = INTERIOR)
+    {
+      DofFace face = el->getFace(ith);
+      int elIndex = el->getIndex();
+      ElementObjectData elObj(elIndex, ith, bound);
+
+      faceElements[face].push_back(elObj);
+      faceLocalMap[elObj] = face;
     }
 
-    void addFace(DofFace face, 
-		 int elIndex, int ith, BoundaryType bound = INTERIOR)
+
+    void addElement(Element *el, BoundaryType bound = INTERIOR)
     {
-      faceElements[face].push_back(ElementObjectData(elIndex, ith, bound));
+      for (int i = 0; i < el->getGeo(VERTEX); i++)
+	addVertex(el, i);
+
+      for (int i = 0; i < el->getGeo(EDGE); i++)
+	addEdge(el, i);      
+
+      for (int i = 0; i < el->getGeo(FACE); i++)
+	addFace(el, i);
     }
 
+
     void createRankData(std::map<int, int>& macroElementRankMap);
 
+
     bool iterate(GeoIndex pos)
     {
       if (iterGeoPos == CENTER) {
@@ -258,6 +302,21 @@ namespace AMDiS {
       return faceInRank[face];
     }
 
+    DegreeOfFreedom getVertexLocalMap(ElementObjectData &data)
+    {
+      return vertexLocalMap[data];
+    }
+
+    DofEdge getEdgeLocalMap(ElementObjectData &data)
+    {
+      return edgeLocalMap[data];
+    }
+
+    DofFace getFaceLocalMap(ElementObjectData &data)
+    {
+      return faceLocalMap[data];
+    }
+
     void serialize(std::ostream &out);
 
     void deserialize(std::istream &in);
@@ -276,14 +335,22 @@ namespace AMDiS {
     std::map<DofEdge, std::vector<ElementObjectData> > edgeElements;
     std::map<DofFace, std::vector<ElementObjectData> > faceElements;
 
+
+    std::map<ElementObjectData, DegreeOfFreedom> vertexLocalMap;
+    std::map<ElementObjectData, DofEdge> edgeLocalMap;
+    std::map<ElementObjectData, DofFace> faceLocalMap;
+
+
     std::map<DegreeOfFreedom, int> vertexOwner;
     std::map<DofEdge, int> edgeOwner;
     std::map<DofFace, int> faceOwner;
 
+
     std::map<DegreeOfFreedom, std::map<int, ElementObjectData> > vertexInRank;
     std::map<DofEdge, std::map<int, ElementObjectData> > edgeInRank;
     std::map<DofFace, std::map<int, ElementObjectData> > faceInRank;
 
+
     std::map<DegreeOfFreedom, std::map<int, ElementObjectData> >::iterator vertexIter;
     std::map<DofEdge, std::map<int, ElementObjectData> >::iterator edgeIter;
     std::map<DofFace, std::map<int, ElementObjectData> >::iterator faceIter;

AMDiS/src/parallel/InteriorBoundary.h

@@ -44,13 +44,13 @@ namespace AMDiS {
 	excludedSubstructures(0)
     {}
 
-    BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj)
+    BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj, bool rMode = false)
       : el(e),
 	elIndex(e->getIndex()),
 	elType(eType),
 	subObj(sObj),
 	ithObj(iObj),
-	reverseMode(false),
+	reverseMode(rMode),
 	excludedSubstructures(0)
     {}
 

AMDiS/src/parallel/MeshDistributor.cc

-
+//
 // Software License for AMDiS
 //
 // Copyright (c) 2010 Dresden University of Technology 
@@ -68,6 +68,7 @@ namespace AMDiS {
       deserialized(false),
       writeSerializationFile(false),
       repartitioningAllowed(false),
+      repartitionIthChange(20),
       nTimestepsAfterLastRepartitioning(0),
       repartCounter(0),
       debugOutputDir(""),
@@ -83,7 +84,9 @@ namespace AMDiS {
     GET_PARAMETER(0, name + "->repartitioning", "%d", &tmp);
     repartitioningAllowed = (tmp > 0);
 
-    GET_PARAMETER(0, name + "->debug output dir", &debugOutputDir);    
+    GET_PARAMETER(0, name + "->debug output dir", &debugOutputDir);
+
+    GET_PARAMETER(0, name + "->repartition ith change", "%d", &repartitionIthChange);
   }
 
 
@@ -581,7 +584,7 @@ namespace AMDiS {
     nTimestepsAfterLastRepartitioning++;
 
     if (repartitioningAllowed) {
-      if (nTimestepsAfterLastRepartitioning >= 20) {
+      if (nTimestepsAfterLastRepartitioning >= repartitionIthChange) {
 	repartitionMesh();
 	nTimestepsAfterLastRepartitioning = 0;
       }
@@ -1001,7 +1004,9 @@ namespace AMDiS {
       return;
 
 #if (DEBUG != 0)
-    if (repartCounter == 0) {
+    ParallelDebug::testDoubleDofs(mesh);
+
+    if (repartCounter == 0) {    
       std::stringstream oss;
       oss << debugOutputDir << "partitioning-" << repartCounter << ".vtu";
       DOFVector<double> tmpa(feSpace, "tmp");
@@ -1191,8 +1196,8 @@ namespace AMDiS {
 
     // === Remove double DOFs. ===
 
-    MeshManipulation meshManipulation(mesh);
-    meshManipulation.deleteDoubleDofs(newMacroEl);
+    MeshManipulation meshManipulation(feSpace);
+    meshManipulation.deleteDoubleDofs(newMacroEl, elObjects);
 
     mesh->dofCompress();
 
@@ -1207,21 +1212,24 @@ namespace AMDiS {
     MSG("USED-SIZE B: %d\n", mesh->getDofAdmin(0).getUsedDofs());
 
     ParallelDebug::testAllElements(*this);
+    ParallelDebug::testDoubleDofs(mesh);
 #endif
 
     partitioner->fillCoarsePartitionVec(&partitionVec);
 
     updateInteriorBoundaryInfo();
 
+#if (DEBUG != 0)
+    ParallelDebug::printBoundaryInfo(*this);
+#endif
+
     updateLocalGlobalNumbering();
 
-   #if (DEBUG != 0)
+#if (DEBUG != 0)
     MSG("AMDiS runs in debug mode, so make some test ...\n");
 
     ParallelDebug::testAllElements(*this);
     ParallelDebug::testInteriorBoundary(*this);
-    ParallelDebug::testCommonDofs(*this, true);
-    ParallelDebug::testGlobalIndexByCoords(*this);
 
     debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
 
@@ -1273,14 +1281,7 @@ namespace AMDiS {
       
       // === Add all sub object of the element to the variable elObjects. ===
 
-      for (int i = 0; i < el->getGeo(VERTEX); i++)
-	elObjects.addVertex(el->getDof(i, 0), el->getIndex(), i);
-
-      for (int i = 0; i < el->getGeo(EDGE); i++)
-	elObjects.addEdge(el->getEdge(i), el->getIndex(), i);      
-
-      for (int i = 0; i < el->getGeo(FACE); i++)
-	elObjects.addFace(el->getFace(i), el->getIndex(), i);
+      elObjects.addElement(el);
 
 
       // === Get periodic boundary information. === 
@@ -1801,7 +1802,8 @@ namespace AMDiS {
     // === Update dof admins due to new number of dofs. ===
   
     lastMeshChangeIndex = mesh->getChangeIndex();
-    
+
+
 #if (DEBUG != 0)
     std::stringstream oss;
     oss << debugOutputDir << "elementIndex-" << mpiRank << ".vtu";

AMDiS/src/parallel/MeshDistributor.h

@@ -568,6 +568,8 @@ namespace AMDiS {
     /// If true, it is possible to repartition the mesh during computations.
     bool repartitioningAllowed;
 
+    int repartitionIthChange;
+
     int nTimestepsAfterLastRepartitioning;
 
     int repartCounter;

AMDiS/src/parallel/MeshManipulation.cc

@@ -12,118 +12,118 @@
 
 #include "parallel/MeshManipulation.h"
 #include "Mesh.h"
+#include "BasisFunction.h"
 #include "Traverse.h"
 #include "Debug.h"
 
 namespace AMDiS {
 
-  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl)
+  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, 
+					  ElementObjects &objects)
   {
-    std::map<int, MacroElement*> leafInMacroEl;
-    
+    FUNCNAME("MeshManipulation::deleteDoubleDofs()");
+
+    TEST_EXIT(mesh->getDim() == 2)("Not yet supported for dim != 2!\n");
+
+    std::map<int, MacroElement*> macroIndexMap;
+    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
+	 it != newMacroEl.end(); ++it)
+      macroIndexMap[(*it)->getIndex()] = *it;
+
+    std::set<int> macrosProcessed;
+    std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
+
+
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
     while (elInfo) {
-      leafInMacroEl[elInfo->getElement()->getIndex()] = elInfo->getMacroElement();
+      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
+	int index = elInfo->getMacroElement()->getIndex();
+
+	macrosProcessed.insert(index);
+	macroIndexMap[index] = elInfo->getMacroElement();
+      }
+
       elInfo = stack.traverseNext(elInfo);
     }
 
-    deleteDoubleDofs(leafInMacroEl, newMacroEl, 0);
-    deleteDoubleDofs(leafInMacroEl, newMacroEl, 1);
-  }
 
+    for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); 
+	 it != newMacroEl.end(); ++it) {
 
-  void MeshManipulation::deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl,
-					  std::set<MacroElement*>& newMacroEl,
-					  int mode)
-  {
-    FUNCNAME("MeshManipulation::deleteDoubleDofs()");
+      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
+	ElementObjectData elObj((*it)->getIndex(), i);
+	DegreeOfFreedom vertex = objects.getVertexLocalMap(elObj);
+	std::vector<ElementObjectData> &vertexEl = objects.getElements(vertex);
 
-    std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
-    std::map<const DegreeOfFreedom*, int> mapDofsMacro;
+	for (std::vector<ElementObjectData>::iterator elIt = vertexEl.begin();
+	     elIt != vertexEl.end(); ++elIt) {
+	  if (elIt->elIndex == (*it)->getIndex())
+	    continue;
 
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
-    while (elInfo) {
-      Element *el = elInfo->getElement();
-
-      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
-
-	Element *neigh = elInfo->getNeighbour(i);
+	  if (macrosProcessed.count(elIt->elIndex) == 1) {
+	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1)
+	      ("Should not happen!\n");
+
+	    Element *el0 = (*it)->getElement();
+	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();
+