Parallel AMDiS

AMDiS/src/InteriorBoundary.cc

 #include "InteriorBoundary.h"
 
 namespace AMDiS {
+
+  AtomicBoundary& InteriorBoundary::getNewAtomicBoundary(int rank)
+  {
+    boundary[rank].resize(boundary[rank].size() + 1);
+    return boundary[rank][boundary[rank].size() - 1];
+  }
+
 }

AMDiS/src/InteriorBoundary.h

@@ -23,6 +23,7 @@
 #define AMDIS_INTERIORBOUNDARY_H
 
 #include <vector>
+#include <map>
 
 #include "MacroElement.h"
 
@@ -32,7 +33,7 @@ namespace AMDiS {
   struct BoundaryObject {
 
     /// The macro element to which the boundary element corresponds to.
-    MacroElement& el;
+    Element* el;
 
     /** \brief
      * Defines the geometrical object at the boundary. It must be "a part" of the
@@ -58,13 +59,10 @@ namespace AMDiS {
    */
   struct AtomicBoundary {
     /// The rank's part of the boundary.
-    BoundaryObject& rankObject;
+    BoundaryObject rankObject;
 
     /// The object on the other side of the boundary.
-    BoundaryObject& neighbourObject;
-
-    /// The rank number of the process at the other side of the boundary.
-    int neighbourRank;
+    BoundaryObject neighbourObject;
   };
 
   /** \brief
@@ -73,8 +71,12 @@ namespace AMDiS {
    */
   class InteriorBoundary {
   public:
-  protected:
-    std::vector<AtomicBoundary> boundary;
+    InteriorBoundary() {}
+
+    AtomicBoundary& getNewAtomicBoundary(int rank);
+
+  public:
+    std::map<int, std::vector<AtomicBoundary> > boundary;
   };
 
 }

AMDiS/src/ParallelDomainProblem.cc

@@ -45,31 +45,26 @@ namespace AMDiS {
     // and now partition the mesh
     partitionMesh(adaptInfo);   
 
-    std::map<int, std::set<int> > partitionDofs;
 
+    /// === Determine to each dof the set of partitions the dof belongs to. ===
+
+    std::map<const DegreeOfFreedom*, std::set<int> > partitionDofs;
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
-    int nLeaves = 0;
     while (elInfo) {
       Element *element = elInfo->getElement();
 
-      // Hidde elements which are not part of ranks partition.
-      PartitionElementData *partitionData = 
-	dynamic_cast<PartitionElementData*>
-	(element->getElementData(PARTITION_ED));
-      if (partitionData->getPartitionStatus() != IN) {
-      } else {
-      }
-
       // Determine to each dof the partition(s) it corresponds to.
       for (int i = 0; i < 3; i++) 
-	partitionDofs[element->getDOF(i, 0)].insert(partitionVec[element->getIndex()]);
+	partitionDofs[element->getDOF(i)].insert(partitionVec[element->getIndex()]);
           
       elInfo = stack.traverseNext(elInfo);
     }
 
-    std::vector<int> rankDofs;
-    for (std::map<int, std::set<int> >::iterator it = partitionDofs.begin();
+    /// === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
+
+    std::vector<const DegreeOfFreedom*> rankDofs;
+    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDofs.begin();
 	 it != partitionDofs.end();
 	 ++it) {
       for (std::set<int>::iterator itpart1 = it->second.begin();
@@ -79,30 +74,64 @@ namespace AMDiS {
 	  if (it->second.size() == 1) {
 	    rankDofs.push_back(it->first);
 	  } else {	    
+	    // This dof is at the ranks boundary. It is owned by the rank only if
+	    // the rank number is the highest of all ranks containing this dof.
+
 	    bool insert = true;
+	    int highestRank = mpiRank;
 	    for (std::set<int>::iterator itpart2 = it->second.begin();
 		 itpart2 != it->second.end();
 		 ++itpart2) {
-	      if (*itpart2 > mpiRank) {
+	      if (*itpart2 > mpiRank)
 		insert = false;
-		break;
-	      }
+
+	      if (*itpart2 > highestRank)
+		highestRank = *itpart2;
 	    }
-	    if (insert) {
+
+	    if (insert)
 	      rankDofs.push_back(it->first);
-	    }
+
+	    boundaryDofs[it->first] = highestRank;
 	  }
 	}
       }
     }
 
-    if (mpiRank == 1) {
-      std::cout << "RANKS dofs = ";
-      for (int i = 0; i < rankDofs.size(); i++) 
-	std::cout << rankDofs[i] << " ";
-      std::cout << std::endl;
+
+    // === Create interior boundary information ===
+
+    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+    while (elInfo) {
+      Element *element = elInfo->getElement();
+
+      // Hidde elements which are not part of ranks partition.
+      PartitionElementData *partitionData = 
+	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
+      if (partitionData->getPartitionStatus() == IN) {
+	for (int i = 0; i < 3; i++) {
+	  if (!elInfo->getNeighbour(i))
+	    continue;
+
+	  PartitionElementData *neighbourPartitionData =
+	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
+ 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
+ 	    AtomicBoundary& bound = interiorBoundary.
+	      getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]);
+ 	    bound.rankObject.el = element;
+ 	    bound.rankObject.subObjAtBoundary = EDGE;
+ 	    bound.rankObject.ithObjAtBoundary = i;
+ 	    bound.neighbourObject.el = elInfo->getNeighbour(i);
+ 	    bound.neighbourObject.subObjAtBoundary = EDGE;
+ 	    bound.neighbourObject.ithObjAtBoundary = -1;
+ 	  }
+	}
+      }
+
+      elInfo = stack.traverseNext(elInfo);
     }
 
+
     // === Remove all macro elements that are not part of the rank partition. ===
 
     std::vector<MacroElement*> macrosToRemove;
@@ -124,15 +153,15 @@ namespace AMDiS {
     int *gOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
     int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
 
-    for (std::vector<int>::iterator it = rankDofs.begin();
+    for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin();
 	 it != rankDofs.end(); ++it) {
-      gOrder[nRankDOFs++] = *it;
+      gOrder[nRankDOFs++] = (*it)[0];
     }
 
     int rstart = 0;
     MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
     rstart -= nRankDOFs;
-    
+   
     for (int i = 0; i < nRankDOFs; i++) {
       lOrder[i] = rstart + i;
     }
@@ -141,6 +170,103 @@ namespace AMDiS {
     
     free(gOrder);
     free(lOrder);
+
+    /// === Create information which dof indices must be send and which must be received. ===
+
+    std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> > sendNewDofs;
+    std::map<int, std::vector<DegreeOfFreedom> > recvNewDofs;
+
+    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDofs.begin();
+	 it != boundaryDofs.end();
+	 ++it) {
+      if (it->second == mpiRank) {
+	int oldDofIndex = (it->first)[0];
+	int newDofIndex = 0;
+	for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) {
+	  if (rankDofs[i] == it->first) {
+	    newDofIndex = rstart + i;
+	    break;
+	  }
+	}
+
+	for (std::set<int>::iterator itRanks = partitionDofs[it->first].begin();
+	     itRanks != partitionDofs[it->first].end();
+	     ++itRanks) {
+	  if (*itRanks != mpiRank) {
+	    sendNewDofs[*itRanks][oldDofIndex] = newDofIndex;
+	  }
+	}
+      } else {
+	recvNewDofs[it->second].push_back((it->first)[0]);
+      }
+    }
+
+    /// === Send and receive the dof indices at boundary. ===
+
+    std::vector<int*> sendBuffers(sendNewDofs.size());
+    std::vector<int*> recvBuffers(recvNewDofs.size());
+    
+    int i = 0;
+    for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
+	 sendIt != sendNewDofs.end();
+	 ++sendIt, i++) {
+      sendBuffers[i] = new int[sendIt->second.size() * 2];
+      int c = 0;
+      for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = sendIt->second.begin();
+	   dofIt != sendIt->second.end();
+	   ++dofIt, c += 2) {
+	sendBuffers[i][c] = dofIt->first;
+	sendBuffers[i][c + 1] = dofIt->second;
+      }
+
+      mpiComm.Isend(sendBuffers[i], sendIt->second.size() * 2, MPI_INT, sendIt->first, 0);
+    }
+
+    i = 0;
+    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end();
+	 ++recvIt, i++) {
+      recvBuffers[i] = new int[recvIt->second.size() * 2];
+      
+      mpiComm.Irecv(recvBuffers[i], recvIt->second.size() * 2, MPI_INT, recvIt->first, 0);
+    }
+
+
+    mpiComm.Barrier();
+
+    /// === Change dof indices at boundary from other ranks. ===
+
+    i = 0;
+    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end();
+	 ++recvIt, i++) {
+
+      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
+	for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDofs.begin();
+	     dofIt != boundaryDofs.end();
+	     ++dofIt) {
+	  if ((dofIt->first)[0] == recvBuffers[i][j * 2]) {
+	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = recvBuffers[i][j * 2 + 1];
+	    break;
+	  }
+	}
+      }
+
+      delete [] recvBuffers[i];
+    }
+
+    i = 0;
+    for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
+	 sendIt != sendNewDofs.end();
+	 ++sendIt, i++) {
+      delete [] sendBuffers[i];
+    }
+
+    /// === Change dof indices for rank partition. ===
+
+    for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) {
+      const_cast<DegreeOfFreedom*>(rankDofs[i])[0] = rstart + i;
+    }
   }
 
   void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)

AMDiS/src/ParallelDomainProblem.h

@@ -29,6 +29,7 @@
 #include "ProblemIterationInterface.h"
 #include "FiniteElemSpace.h"
 #include "AdaptInfo.h"
+#include "InteriorBoundary.h"
 #include "petscao.h"
 #include "mpi.h"
 
@@ -156,6 +157,14 @@ namespace AMDiS {
 
     /// Number of DOFs in the rank mesh.
     int nRankDOFs;
+
+    /** \brief 
+     * Defines the interioir boundaries of the domain that result from partitioning
+     * the whole mesh.
+     */
+    InteriorBoundary interiorBoundary;
+
+    std::map<const DegreeOfFreedom*, int> boundaryDofs;
   };
 
   class ParallelDomainProblemScal : public ParallelDomainProblemBase