Work on general mesh partitioner interface.

AMDiS/src/Mesh.cc

@@ -333,9 +333,9 @@ namespace AMDiS {
 
       // Go through all neighbours of the macro element and remove this macro element
       // to be neighbour of some other macro element.
-      for (int i = 0; i <= dim; i++)
+      for (int i = 0; i < getGeo(NEIGH); i++)
 	if ((*macroIt)->getNeighbour(i))
-	  for (int j = 0; j <= dim; j++)
+	  for (int j = 0; j < getGeo(NEIGH); j++)
 	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
 	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
 
@@ -352,13 +352,13 @@ namespace AMDiS {
 
       // We will delete at least some element DOFs in the next but will keep the
       // element object still in memory. Therefore the DOF pointer must be set to be
-      // invalide.
+      // invalid.
       (*macroIt)->getElement()->setDofValid(false);
     }     
 
 
-    // === Check now all the dofs, that have no owner anymore and therefore have ===
-    // === to be removed.                                                        ===
+    // === Check now all the DOFs that have no owner anymore and therefore have ===
+    // === to be removed.                                                       ===
 
     for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
 	 dofsIt != dofsOwner.end(); ++dofsIt) {

AMDiS/src/io/DataCollector.cc

@@ -217,11 +217,11 @@ namespace AMDiS {
     FUNCNAME("DataCollector::addElementData()");
 
     const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
-    DegreeOfFreedom vertexDOF;
-    WorldVector<double> vertexCoords;
 
-    // MSG("ELEMENT: %d\n", elInfo->getElement()->getIndex());
-    // MSG("DOFs: %d %d %d\n", dof[0][0], dof[1][0], dof[2][0]);
+//     MSG("ELEMENT: %d %d\n", 
+// 	elInfo->getElement()->getIndex(),
+// 	elInfo->getMacroElement()->getIndex());    
+//     MSG("DOFs: %d %d %d %d\n", dof[0][0], dof[1][0], dof[2][0], dof[3][0]);
 
     // create ElementInfo
     ElementInfo elementInfo(dim);
@@ -244,28 +244,28 @@ namespace AMDiS {
     }
 
     // for all vertices
-    for (int i = 0; i < dim + 1; i++) {
+    for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
       // get coords of this vertex
-      vertexCoords = elInfo->getCoord(i);
+      WorldVector<double> &vertexCoords = elInfo->getCoord(i);
       
       // get dof index of this vertex
-      vertexDOF = dof[i][nPreDofs];
+      DegreeOfFreedom vertexDof = dof[i][nPreDofs];
 
       // search for coords at this dof
       std::list<VertexInfo>::iterator it =
-	  find((*vertexInfos)[vertexDOF].begin(), (*vertexInfos)[vertexDOF].end(),
+	  find((*vertexInfos)[vertexDof].begin(), (*vertexInfos)[vertexDof].end(),
 	       vertexCoords);
       
       // coords not yet in list?
-      if (it == (*vertexInfos)[vertexDOF].end()) {
+      if (it == (*vertexInfos)[vertexDof].end()) {
 	// create new vertex info and increment number of vertices
 	VertexInfo newVertexInfo = {vertexCoords, nVertices};
 
 	// add new vertex info to list
-	(*vertexInfos)[vertexDOF].push_front(newVertexInfo);
+	(*vertexInfos)[vertexDof].push_front(newVertexInfo);
 
 	// set iterator to new vertex info
-	it = (*vertexInfos)[vertexDOF].begin();	
+	it = (*vertexInfos)[vertexDof].begin();	
 
 	nVertices++;
       }

AMDiS/src/io/VtkWriter.cc

@@ -56,16 +56,6 @@ namespace AMDiS {
     file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc)); 
     writeFileToStream(file);
 
-#if 0
-
-    std::ofstream file;
-    file.open(name.c_str());
-    TEST_EXIT(file.is_open())("Cannot open file %s for writing\n", name.c_str());
-    writeFileToStream(file);    
-    file.close();
-
-#endif
-      
     return 0;
   }
 
@@ -159,7 +149,7 @@ namespace AMDiS {
     DataCollector dc(values->getFeSpace(), values);
     std::vector<DataCollector*> dcList(0);
     dcList.push_back(&dc);
-    writeFile(dcList,filename,writeParallel);
+    writeFile(dcList, filename, writeParallel);
   }
 
 

AMDiS/src/io/VtkWriter.hh

@@ -19,6 +19,8 @@ namespace AMDiS {
   template<typename T>
   void VtkWriter::writeFileToStream(T &file)
   {
+    FUNCNAME("VtkWriter::writeFileToStream()");
+
     int nVertices = (*dataCollector)[0]->getNumberVertices();
     int nElements = (*dataCollector)[0]->getNumberElements();
     int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);

AMDiS/src/parallel/ElementObjectData.cc

@@ -336,7 +336,7 @@ namespace AMDiS {
     FUNCNAME("ElementObjects::provideConnectedPeriodicBoundary()");
 
     std::pair<BoundaryType, BoundaryType> bConn = 
-      (b0 <= b1 ? std::make_pair(b0, b1) : std::make_pair(b1, b0));
+      (b0 <= b1 ? make_pair(b0, b1) : make_pair(b1, b0));
 
     if (bConnMap.count(bConn) == 0) {
       BoundaryType newPeriodicBoundaryType = getNewBoundaryType();
@@ -417,6 +417,65 @@ namespace AMDiS {
   }
 
 
+  void ElementObjects::createReverseModeData(map<int, Element*> &elIndexMap,
+					     map<int, int> &elIndexTypeMap)
+  {
+    FUNCNAME("ElementObjects::createReverseModeData()");
+
+    // In 2D, all reverse modes are always true!
+    if (mesh->getDim() == 2)
+      return;
+
+    for (map<DofEdge, vector<ElementObjectData> >::iterator edgeIt = edgeElements.begin();
+	 edgeIt != edgeElements.end(); ++edgeIt) {
+
+      vector<ElementObjectData>& els = edgeIt->second;
+
+      for (unsigned int i = 0; i < els.size(); i++) {
+	BoundaryObject obj0(elIndexMap[els[i].elIndex], 
+			    elIndexTypeMap[els[i].elIndex], 
+			    EDGE, els[i].ithObject);
+
+	for (unsigned int j = i + 1; j < els.size(); j++) {
+	  BoundaryObject obj1(elIndexMap[els[j].elIndex], 
+			      elIndexTypeMap[els[j].elIndex], 
+			      EDGE, els[j].ithObject);
+
+	  bool reverseMode = 
+	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);
+
+	  edgeReverseMode[make_pair(els[i], els[j])] = reverseMode;
+	  edgeReverseMode[make_pair(els[j], els[i])] = reverseMode;
+	}
+      }
+    }
+
+    for (map<DofFace, vector<ElementObjectData> >::iterator faceIt = faceElements.begin();
+	 faceIt != faceElements.end(); ++faceIt) {
+
+      vector<ElementObjectData>& els = faceIt->second;
+
+      for (unsigned int i = 0; i < els.size(); i++) {
+	BoundaryObject obj0(elIndexMap[els[i].elIndex], 
+			    elIndexTypeMap[els[i].elIndex], 
+			    FACE, els[i].ithObject);
+
+	for (unsigned int j = i + 1; j < els.size(); j++) {
+	  BoundaryObject obj1(elIndexMap[els[j].elIndex], 
+			      elIndexTypeMap[els[j].elIndex], 
+			      FACE, els[j].ithObject);
+
+	  bool reverseMode = 
+	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);
+
+	  faceReverseMode[make_pair(els[i], els[j])] = reverseMode;
+	  faceReverseMode[make_pair(els[j], els[i])] = reverseMode;
+	}
+      }
+    }
+  }
+
+
   void ElementObjects::serialize(ostream &out)
   {
     FUNCNAME("ElementObjects::serialize()");
@@ -505,6 +564,42 @@ namespace AMDiS {
     SerUtil::serialize(out, periodicVertices);
     SerUtil::serialize(out, periodicEdges);
     SerUtil::serialize(out, periodicFaces);
+
+
+    nSize = periodicDofAssoc.size();
+    SerUtil::serialize(out, nSize);
+    for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
+	 it != periodicDofAssoc.end(); ++it) {
+      SerUtil::serialize(out, it->first);
+      SerUtil::serialize(out, it->second);
+    }
+
+    nSize = periodicEdgeAssoc.size();
+    SerUtil::serialize(out, nSize);
+    for (map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
+	 it != periodicEdgeAssoc.end(); ++it) {
+      SerUtil::serialize(out, it->first);
+      SerUtil::serialize(out, it->second);
+    }
+
+
+    nSize = edgeReverseMode.size();
+    SerUtil::serialize(out, nSize);
+    for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = edgeReverseMode.begin();
+	 it != edgeReverseMode.end(); ++it) {
+      it->first.first.serialize(out);
+      it->first.second.serialize(out);
+      SerUtil::serialize(out, it->second);
+    }
+
+    nSize = faceReverseMode.size();
+    SerUtil::serialize(out, nSize);
+    for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = faceReverseMode.begin();
+	 it != faceReverseMode.end(); ++it) {
+      it->first.first.serialize(out);
+      it->first.second.serialize(out);
+      SerUtil::serialize(out, it->second);
+    }
   }
 
 
@@ -613,7 +708,55 @@ namespace AMDiS {
 
     SerUtil::deserialize(in, periodicVertices);
     SerUtil::deserialize(in, periodicEdges);
-    SerUtil::deserialize(in, periodicFaces);  
+    SerUtil::deserialize(in, periodicFaces);
+
+    
+
+    SerUtil::deserialize(in, nSize);
+    periodicDofAssoc.clear();
+    for (int i = 0; i < nSize; i++) {
+      DegreeOfFreedom dof;
+      std::set<DegreeOfFreedom> dofs;
+      SerUtil::deserialize(in, dof);
+      SerUtil::deserialize(in, dofs);
+      periodicDofAssoc[dof] = dofs;
+    }
+
+    SerUtil::deserialize(in, nSize);
+    periodicEdgeAssoc.clear();
+    for (int i = 0; i < nSize; i++) {
+      DofEdge edge;
+      std::set<DofEdge> edges;
+      SerUtil::deserialize(in, edge);
+      SerUtil::deserialize(in, edges);
+      periodicEdgeAssoc[edge] = edges;
+    }
+
+
+
+    SerUtil::deserialize(in, nSize);
+    edgeReverseMode.clear();
+    for (int i = 0; i < nSize; i++) {
+      ElementObjectData obj0, obj1;
+      bool reverseMode;
+      obj0.deserialize(in);
+      obj1.deserialize(in);
+      SerUtil::deserialize(in, reverseMode);
+
+      edgeReverseMode[make_pair(obj0, obj1)] = reverseMode;
+    }
+
+    SerUtil::deserialize(in, nSize);
+    faceReverseMode.clear();
+    for (int i = 0; i < nSize; i++) {
+      ElementObjectData obj0, obj1;
+      bool reverseMode;
+      obj0.deserialize(in);
+      obj1.deserialize(in);
+      SerUtil::deserialize(in, reverseMode);
+
+      faceReverseMode[make_pair(obj0, obj1)] = reverseMode;
+    }
   }
 
 

AMDiS/src/parallel/ElementObjectData.h

@@ -148,6 +148,9 @@ namespace AMDiS {
      */
     void createRankData(map<int, int>& macroElementRankMap);
 
+    void createReverseModeData(map<int, Element*> &elIndexMap,
+			       map<int, int> &elIndexTypeMap);
+
 
     /** \brief
      * Iterates over all elements for one geometrical index, i.e., over all vertices,
@@ -320,6 +323,32 @@ namespace AMDiS {
       return faceElements[face];
     }
 
+
+    /// Returns a vector with all macro elements that have a given vertex DOF in common.
+    vector<ElementObjectData>& getElementsVertex(int elIndex, int ithVertex)
+    {
+      ElementObjectData elObj(elIndex, ithVertex);
+      DegreeOfFreedom vertex = vertexLocalMap[elObj];
+      return vertexElements[vertex];
+    }
+    
+    /// Returns a vector with all macro elements that have a given edge in common.
+    vector<ElementObjectData>& getElementsEdge(int elIndex, int ithEdge)
+    {
+      ElementObjectData elObj(elIndex, ithEdge);
+      DofEdge edge = edgeLocalMap[elObj];
+      return edgeElements[edge];
+    }
+
+    /// Returns a vector with all macro elements that have a given face in common.
+    vector<ElementObjectData>& getElementsFace(int elIndex, int ithFace)
+    {
+      ElementObjectData elObj(elIndex, ithFace);
+      DofFace face = faceLocalMap[elObj];
+      return faceElements[face];
+    }
+
+
     
     /// Returns a map that maps to each rank all macro elements in this rank that
     /// have a given vertex DOF in common.
@@ -381,6 +410,22 @@ namespace AMDiS {
       return periodicFaces;
     }
 
+    bool getEdgeReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    {
+      TEST_EXIT_DBG(edgeReverseMode.count(make_pair(obj0, obj1)))
+	("Should not happen!\n");
+
+      return edgeReverseMode[make_pair(obj0, obj1)];
+    }
+
+    bool getFaceReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    {
+      TEST_EXIT_DBG(faceReverseMode.count(make_pair(obj0, obj1)))
+	("Should not happen!\n");
+
+      return faceReverseMode[make_pair(obj0, obj1)];
+    }
+
     /// Write the element database to disk.
     void serialize(ostream &out);
     
@@ -505,7 +550,7 @@ namespace AMDiS {
     /// iterators, i.e., if no iteration is running.
     GeoIndex iterGeoPos;
 
-    std::map<std::pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
+    map<pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
 
     // The following three data structures store periodic DOFs, edges and faces.
     PerBoundMap<DegreeOfFreedom>::type periodicVertices;
@@ -513,10 +558,14 @@ namespace AMDiS {
     PerBoundMap<DofFace>::type periodicFaces;
 
     // Stores to each vertex all its periodic associations.
-    std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
+    map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
 
     // Stores to each edge all its periodic associations.
-    std::map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
+    map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
+
+    map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode;
+
+    map<pair<ElementObjectData, ElementObjectData>, bool> faceReverseMode;
   };
 
 }

AMDiS/src/parallel/InteriorBoundary.cc

@@ -18,35 +18,38 @@
 
 namespace AMDiS {
 
-  void BoundaryObject::setReverseMode(BoundaryObject &otherBound, 
-				      FiniteElemSpace *feSpace,
-				      BoundaryType boundary)
+
+  bool BoundaryObject::computeReverseMode(BoundaryObject &obj0, 
+					  BoundaryObject &obj1,
+					  FiniteElemSpace *feSpace,
+					  BoundaryType boundary)
   {
-    FUNCNAME("BoundaryObject::setReverseMode()");
+    FUNCNAME("BoundaryObject::computeReverseMode()");
 
-    bool otherMode = false;
+    bool reverseMode = false;
 
     switch (feSpace->getMesh()->getDim()) {
     case 2:
-      otherMode = true;
+      reverseMode = true;
       break;
 
     case 3:
-      TEST_EXIT_DBG(otherBound.elType == 0)
-	("Only 3D macro elements with level 0 are supported. This element has level %d!\n", otherBound.elType);
+      TEST_EXIT_DBG(obj1.elType == 0)
+	("Only 3D macro elements with level 0 are supported. This element has level %d!\n", 
+	 obj1.elType);
 
 
-      if (subObj == EDGE) {	
-	int el0_v0 = el->getVertexOfEdge(ithObj, 0);
-	int el0_v1 = el->getVertexOfEdge(ithObj, 1);
-	int el1_v0 = el->getVertexOfEdge(otherBound.ithObj, 0);
-	int el1_v1 = el->getVertexOfEdge(otherBound.ithObj, 1);
+      if (obj0.subObj == EDGE) {	
+	int el0_v0 = obj0.el->getVertexOfEdge(obj0.ithObj, 0);
+	int el0_v1 = obj0.el->getVertexOfEdge(obj0.ithObj, 1);
+	int el1_v0 = obj0.el->getVertexOfEdge(obj1.ithObj, 0);
+	int el1_v1 = obj0.el->getVertexOfEdge(obj1.ithObj, 1);
 
 	const BasisFunction *basFcts = feSpace->getBasisFcts();
 	int nBasFcts = basFcts->getNumber();
 	std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
-	basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
-	basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
+	basFcts->getLocalIndices(obj0.el, feSpace->getAdmin(), localDofs0);
+	basFcts->getLocalIndices(obj1.el, feSpace->getAdmin(), localDofs1);
 
 	Mesh *mesh = feSpace->getMesh();
 
@@ -59,30 +62,25 @@ namespace AMDiS {
 	    ("This should not happen!\n");
 
 	  if (localDofs0[el0_v0] != localDofs1[el1_v0])
-	    otherMode = true; 	
+	    reverseMode = true; 	
 	} else {
 	  if (mesh->associated(localDofs0[el0_v0], localDofs1[el1_v0]) == false)
-	    otherMode = true;
-
-	  if ((elIndex == 28 && otherBound.elIndex == 41) ||
-	      (elIndex == 41 && otherBound.elIndex == 28)) {
-	    MSG("HERE TEST B (%d %d): %d\n", el0_v0, el1_v0, otherMode);
-	  }
+	    reverseMode = true;
 	}
       }
 
-      if (subObj == FACE && ithObj != 1) {
+      if (obj0.subObj == FACE && obj0.ithObj != 1) {
 	const BasisFunction *basFcts = feSpace->getBasisFcts();
 	int nBasFcts = basFcts->getNumber();
 	std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
-	basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
-	basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
+	basFcts->getLocalIndices(obj0.el, feSpace->getAdmin(), localDofs0);
+	basFcts->getLocalIndices(obj1.el, feSpace->getAdmin(), localDofs1);
 	
-	if (ithObj == 2 || ithObj == 3)
-	  otherMode = (localDofs0[0] != localDofs1[0]);
+	if (obj0.ithObj == 2 || obj0.ithObj == 3)
+	  reverseMode = (localDofs0[0] != localDofs1[0]);
 	  
-	if (ithObj == 0)
-	  otherMode = (localDofs0[1] != localDofs1[1]);
+	if (obj0.ithObj == 0)
+	  reverseMode = (localDofs0[1] != localDofs1[1]);
       }
       break;
 
@@ -90,14 +88,23 @@ namespace AMDiS {
       ERROR_EXIT("This should not happen!\n");
     }