FETI-DP

AMDiS/src/parallel/FeSpaceMapping.cc

@@ -11,6 +11,7 @@
 
 
 #include "parallel/FeSpaceMapping.h"
+#include "parallel/StdMpi.h"
 
 namespace AMDiS {
 
@@ -28,7 +29,7 @@ namespace AMDiS {
   void GlobalDofMap::update(bool add)
   {
     for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
-      if (it->second == -1)
+      if (it->second == -1 && nonRankDofs.count(it->first) == 0)
 	it->second = nRankDofs++;
     
     nOverallDofs = 0;
@@ -37,6 +38,9 @@ namespace AMDiS {
     
     if (add)
       addOffset(rStartDofs);
+
+    if (overlap)
+      computeNonLocalIndices();
   }
 
 
@@ -46,4 +50,55 @@ namespace AMDiS {
       it->second += offset;
   }
 
+  void GlobalDofMap::computeNonLocalIndices()
+  {
+    FUNCNAME("GlobalDofMap::computeNonLocalIndices()");
+
+    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
+
+    StdMpi<vector<int> > stdMpi(*mpiComm);
+    for (DofComm::Iterator it(*sendDofs, feSpace); !it.end(); it.nextRank())
+      for (; !it.endDofIter(); it.nextDof())
+	if (dofMap.count(it.getDofIndex()) && !nonRankDofs.count(it.getDofIndex()))
+	  stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()]);
+
+    stdMpi.updateSendDataSize();
+
+    for (DofComm::Iterator it(*recvDofs, feSpace); !it.end(); it.nextRank()) {
+      bool recvFromRank = false;
+      for (; !it.endDofIter(); it.nextDof()) {
+	if (nonRankDofs.count(it.getDofIndex())) {
+	  recvFromRank = true;
+	  break;
+	}
+      }
+
+      if (recvFromRank) 
+	stdMpi.recv(it.getRank());
+    }
+
+    stdMpi.startCommunication();
+
+    for (DofComm::Iterator it(*recvDofs, feSpace);
+	 !it.end(); it.nextRank()) {
+      int i = 0;
+      for (; !it.endDofIter(); it.nextDof())
+	if (nonRankDofs.count(it.getDofIndex()))
+	  dofMap[it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[i++];
+    }
+  }
+
+
+  void GlobalDofMap::print()
+  {
+    FUNCNAME("GlobalDofMap::print()");
+    
+    MSG("Local to global mapping on this rank: \n");
+
+    for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
+      if (nonRankDofs.count(it->first) == 0)	
+	MSG("  %d -> %d  (rank-dof)\n", it->first, it->second);
+      else
+	MSG("  %d -> %d  \n", it->first, it->second);
+  }
 }

AMDiS/src/parallel/FeSpaceMapping.h

@@ -22,6 +22,8 @@
 
 #include <vector>
 #include <map>
+#include <set>
+#include "parallel/DofComm.h"
 #include "parallel/MpiHelper.h"
 #include "parallel/ParallelTypes.h"
 
@@ -43,9 +45,11 @@ namespace AMDiS {
       
     GlobalDofMap(MPI::Intracomm* m)
       : mpiComm(m),
+	feSpace(NULL),
 	nRankDofs(0),
 	nOverallDofs(0),
-	rStartDofs(0)
+	rStartDofs(0),
+	overlap(false)
     {}
     
     void clear();
@@ -57,22 +61,24 @@ namespace AMDiS {
       return dofMap[d];
     }
     
-    void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
+    void insertRankDof(DegreeOfFreedom dof0)
     {
       FUNCNAME("GlobalDofMap::insertRankDof()");
       
       TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
       
-      dofMap[dof0] = dof1;
+      dofMap[dof0] = -1;
     }
     
-    void insert(DegreeOfFreedom dof0, DegreeOfFreedom dof1)
+    void insert(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
     {
       FUNCNAME("GlobalDofMap::insert()");
       
       TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
       
       dofMap[dof0] = dof1;
+
+      nonRankDofs.insert(dof0);
     }
     
     bool isSet(DegreeOfFreedom dof)
@@ -94,15 +100,44 @@ namespace AMDiS {
 
     void addOffset(int offset);
 
+    void computeNonLocalIndices();
+
+    void print();
+
+    void setFeSpace(const FiniteElemSpace *fe)
+    {
+      feSpace = fe;
+    }
+
+    void setOverlap(bool b)
+    {
+      overlap = b;
+    }
+
+    void setDofComm(DofComm &pSend, DofComm &pRecv)
+    {
+      sendDofs = &pSend;
+      recvDofs = &pRecv;
+    }
+
   private:
     MPI::Intracomm* mpiComm;
-    
+
+    const FiniteElemSpace *feSpace;
+
     /// 
-    DofMapping dofMap;    
+    DofMapping dofMap;
 
+    std::set<DegreeOfFreedom> nonRankDofs;
+
+    DofComm *sendDofs;
+    
+    DofComm *recvDofs;
   public:
     /// 
     int nRankDofs, nOverallDofs, rStartDofs;   
+
+    bool overlap;
   };
   
 
@@ -140,6 +175,8 @@ namespace AMDiS {
 	data.find(feSpace)->second.clear();
       else
 	data.insert(make_pair(feSpace, T(mpiComm)));
+
+      data.find(feSpace)->second.setFeSpace(feSpace);
     }    
 
     int getRankDofs(vector<const FiniteElemSpace*> &fe)
@@ -203,7 +240,6 @@ namespace AMDiS {
       feSpaces = fe;
       for (unsigned int i = 0; i < feSpaces.size(); i++)
 	addFeSpace(feSpaces[i]);
-
     }
 
     void update()

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -251,60 +251,18 @@ namespace AMDiS {
     for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
       if (meshDistributor->getIsRankDof(feSpace, *it))
 	primalDofMap[feSpace].insertRankDof(*it);
+      else
+  	primalDofMap[feSpace].insert(*it);
 
+    primalDofMap[feSpace].setOverlap(true);
+    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
+  				     meshDistributor->getRecvDofs());
     primalDofMap[feSpace].update();
 
     MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
 	primalDofMap[feSpace].nRankDofs, 
 	primalDofMap[feSpace].nOverallDofs);
 
-
-    // === Communicate primal's global index from ranks that own the     ===
-    // === primals to ranks that contain this primals but are not owning ===
-    // === them.                                                         ===
-    
-    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
-
-    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
-    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
-	 !it.end(); it.nextRank())
-      for (; !it.endDofIter(); it.nextDof())
-	if (primalDofMap[feSpace].isSet(it.getDofIndex())) {
-	  DegreeOfFreedom d = primalDofMap[feSpace][it.getDofIndex()];
-	  stdMpi.getSendData(it.getRank()).push_back(d);
-	}
-
-    stdMpi.updateSendDataSize();
-
-    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
-	 !it.end(); it.nextRank()) {
-      bool recvFromRank = false;
-      for (; !it.endDofIter(); it.nextDof()) {
-	if (primals.count(it.getDofIndex()) && 
-	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
-	  recvFromRank = true;
-	  break;
-	}
-      }
-
-      if (recvFromRank) 
-	stdMpi.recv(it.getRank());
-    }
-
-    stdMpi.startCommunication();
-
-    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
-	 !it.end(); it.nextRank()) {
-      int i = 0;
-      for (; !it.endDofIter(); it.nextDof()) {
-	if (primals.count(it.getDofIndex()) && 
-	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
-	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[i++];
-	  primalDofMap[feSpace].insert(it.getDofIndex(), d);
-	}
-      }
-    }
-
     TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
       ("Number of primals %d, but number of global primals on this rank is %d!\n",
        primals.size(), primalDofMap[feSpace].size());
@@ -472,7 +430,7 @@ namespace AMDiS {
 	nLocalInterior++;
       }
     }
-
+    
     localDofMap[feSpace].update(false);
 
     TEST_EXIT_DBG(nLocalInterior + 
@@ -919,13 +877,7 @@ namespace AMDiS {
 
       for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
 	   it != localDofMap[feSpace].getMap().end(); ++it) {
-
-#if 1
 	int petscIndex = localDofMap.mapLocal(it->first, i);
-#else	  
-	int petscIndex = it->second * nComponents + i;
-#endif
-
 	dofVec[it->first] = localSolB[petscIndex];
       }
 
@@ -1106,12 +1058,7 @@ namespace AMDiS {
 		
 	      if (rowIndex < nLocalInterior) {
 		if (colIndex < nLocalInterior) {
-#if 1
 		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
-#else
-		  int colIndex2 = localDofMap[feSpace][col(*icursor)] * nComponents + j;
-#endif
-
 		  colsLocal.push_back(colIndex2);
 		  valuesLocal.push_back(value(*icursor));
 		} else {
@@ -1123,12 +1070,7 @@ namespace AMDiS {
 		}
 	      } else {
 		if (colIndex < nLocalInterior) {
-#if 1
 		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
-#else
-		  int colIndex2 = localDofMap[feSpace][col(*icursor)] * nComponents + j;
-#endif
-
 		  colsLocalOther.push_back(colIndex2);
 		  valuesLocalOther.push_back(value(*icursor));
 		} else {
@@ -1159,41 +1101,22 @@ namespace AMDiS {
 
 	    if (colsOther.size()) {
 	      for (unsigned int k = 0; k < colsOther.size(); k++)
-#if 1
 		colsOther[k] = localDofMap.mapGlobal(colsOther[k], j);
-#else
-	      colsOther[k] = (localDofMap[feSpace][colsOther[k]] + 
- 				localDofMap[feSpace].rStartDofs) * nComponents + j;
-#endif
-
-
- 	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
+ 	      
+	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	    }
 	  } else {
-#if 1 
 	    int localRowIndex = localDofMap.mapLocal(*cursor, i);
-#else
-	    int localRowIndex = localDofMap[feSpace][*cursor] * nComponents + i; 
-#endif
 
 	    for (unsigned int k = 0; k < cols.size(); k++)
-#if 1
 	      cols[k] = localDofMap.mapLocal(cols[k], j);
-#else
-	      cols[k] = localDofMap[feSpace][cols[k]] * nComponents + j; 
-#endif
-
+	    
   	    MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
   			 &(cols[0]), &(values[0]), ADD_VALUES);
 
 	    if (colsOther.size()) {
-#if 1
 	      int globalRowIndex = localDofMap.mapGlobal(*cursor, i);
-#else
-	      int globalRowIndex = 
-		(localDofMap[feSpace][*cursor] + localDofMap[feSpace].rStartDofs) * nComponents + i;
-#endif
 
 	      for (unsigned int k = 0; k < colsOther.size(); k++)
 		colsOther[k] = 
@@ -1338,12 +1261,7 @@ namespace AMDiS {
 	  double value = *dofIt;
 	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
 	} else {
-#if 1
 	  index = localDofMap.mapGlobal(index, i);
-#else
-	  index = (localDofMap[feSpace][index] + localDofMap[feSpace].rStartDofs) * nComponents + i;
-#endif
-
 	  VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
 	}      
       }

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -712,9 +712,6 @@ namespace AMDiS {
  	  int globalMatIndex =  
  	    globalIndex - meshDistributor->getStartDofs(feSpaces[i]) + offset;
 
-// 	  int globalMatIndex = 
-// 	    globalIndex * feSpaces.size() + i;
-
 	  dofToMatIndex.add(i, globalIndex, globalMatIndex);
 	}