Have a nice weekend.

AMDiS/src/ProblemStat.cc

@@ -870,17 +870,6 @@ namespace AMDiS {
     INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
 		  TIME_USED(first, clock()));
 #endif
-
-#if 0
-    VtkWriter::writeFile(rhs, "rhs0.vtu");
-    for (int i = 0; i < nComponents; i++) {
-      double s = rhs->getDOFVector(i)->sum();
-      s /= -static_cast<double>(rhs->getDOFVector(i)->getUsedSize());
-      MSG("MOVE RHS: %f\n", s);
-      add(*(rhs->getDOFVector(i)), s, *(rhs->getDOFVector(i)));
-    }
-    VtkWriter::writeFile(rhs, "rhs1.vtu");
-#endif
   }
 
 

AMDiS/src/parallel/FeSpaceMapping.cc

@@ -26,17 +26,17 @@ namespace AMDiS {
   }
 
 
-  void GlobalDofMap::update(bool add)
+  void GlobalDofMap::update()
   {
-    for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
-      if (it->second == -1 && nonRankDofs.count(it->first) == 0)
-	it->second = nRankDofs++;
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
+      if (it->second.local == -1 && nonRankDofs.count(it->first) == 0)
+	it->second.local = nRankDofs++;
     
     nOverallDofs = 0;
     rStartDofs = 0;
     mpi::getDofNumbering(*mpiComm, nRankDofs, rStartDofs, nOverallDofs);
     
-    if (add)
+    if (needGlobalMapping)
       addOffset(rStartDofs);
 
     if (overlap)
@@ -48,8 +48,8 @@ namespace AMDiS {
 
   void GlobalDofMap::addOffset(int offset)
   {
-    for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
-      it->second += offset;
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
+      it->second.local += offset;
   }
 
   void GlobalDofMap::computeNonLocalIndices()
@@ -62,7 +62,7 @@ namespace AMDiS {
     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.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].local);
 
     stdMpi.updateSendDataSize();
 
@@ -86,7 +86,7 @@ namespace AMDiS {
       int i = 0;
       for (; !it.endDofIter(); it.nextDof())
 	if (nonRankDofs.count(it.getDofIndex()))
-	  dofMap[it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[i++];
+	  dofMap[it.getDofIndex()].local = stdMpi.getRecvData(it.getRank())[i++];
     }
   }
 
@@ -97,9 +97,9 @@ namespace AMDiS {
 
     map<DegreeOfFreedom, int> dofToMatIndex;
 
-    for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
       if (!nonRankDofs.count(it->first)) {
-	int globalMatIndex = it->second - rStartDofs + offset;
+	int globalMatIndex = it->second.local - rStartDofs + offset;
 	dofToMatIndex[it->first] = globalMatIndex;
       }
     }
@@ -143,10 +143,10 @@ namespace AMDiS {
     
     MSG("Local to global mapping on this rank: \n");
 
-    for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
       if (nonRankDofs.count(it->first) == 0)	
-	MSG("  %d -> %d  (rank-dof)\n", it->first, it->second);
+	MSG("  %d -> %d  (rank-dof)\n", it->first, it->second.local);
       else
-	MSG("  %d -> %d  \n", it->first, it->second);
+	MSG("  %d -> %d  \n", it->first, it->second.local);
   }
 }

AMDiS/src/parallel/FeSpaceMapping.h

@@ -33,6 +33,11 @@
 namespace AMDiS {
   using namespace std;
 
+  struct MultiIndex
+  {
+    int local, global;
+  };
+
   class GlobalDofMap
   {
   public:
@@ -48,6 +53,7 @@ namespace AMDiS {
 	feSpace(NULL),
 	sendDofs(NULL),
 	recvDofs(NULL),
+	needGlobalMapping(false),
 	nRankDofs(0),
 	nOverallDofs(0),
 	rStartDofs(0),
@@ -56,20 +62,23 @@ namespace AMDiS {
     
     void clear();
     
-    DegreeOfFreedom operator[](DegreeOfFreedom d)
+    MultiIndex& operator[](DegreeOfFreedom d)
     {
       TEST_EXIT_DBG(dofMap.count(d))("Should not happen!\n");
 
       return dofMap[d];
     }
     
-    void insertRankDof(DegreeOfFreedom dof0)
+    void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
     {
       FUNCNAME("GlobalDofMap::insertRankDof()");
       
       TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
       
-      dofMap[dof0] = -1;
+      dofMap[dof0].local = dof1;
+      
+      if (dof1 != -1)
+	nRankDofs++;
     }
     
     void insert(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
@@ -78,7 +87,7 @@ namespace AMDiS {
       
       TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
       
-      dofMap[dof0] = dof1;
+      dofMap[dof0].local = dof1;
 
       nonRankDofs.insert(dof0);
     }
@@ -93,12 +102,12 @@ namespace AMDiS {
       return dofMap.size();
     }
     
-    DofMapping& getMap()
+    map<DegreeOfFreedom, MultiIndex>& getMap()
     {
       return dofMap;
     }
     
-    void update(bool add = true);
+    void update();
 
     void addOffset(int offset);
 
@@ -124,19 +133,26 @@ namespace AMDiS {
       recvDofs = &pRecv;
     }
 
+    void setNeedGlobalMapping(bool b)
+    {
+      needGlobalMapping = b;
+    }
+
   private:
     MPI::Intracomm* mpiComm;
 
     const FiniteElemSpace *feSpace;
 
     /// 
-    DofMapping dofMap;
+    map<DegreeOfFreedom, MultiIndex> dofMap;
 
     std::set<DegreeOfFreedom> nonRankDofs;
 
     DofComm *sendDofs;
     
     DofComm *recvDofs;
+
+    bool needGlobalMapping;
   public:
     /// 
     int nRankDofs, nOverallDofs, rStartDofs;   
@@ -157,11 +173,6 @@ namespace AMDiS {
 	rStartDofs(-1)
     {} 
 
-    void setMpiComm(MPI::Intracomm *m)
-    {
-      mpiComm = m;
-    }
-
     T& operator[](const FiniteElemSpace* feSpace)
     {
       FUNCNAME("FeSpaceData::operator[]()");
@@ -239,15 +250,23 @@ namespace AMDiS {
       return rStartDofs;
     }
 
-    void setFeSpaces(vector<const FiniteElemSpace*> &fe)
+    void init(MPI::Intracomm *m,
+	      vector<const FiniteElemSpace*> &fe,
+	      bool needGlobalMapping)
     {
+      mpiComm = m;
       feSpaces = fe;
-      for (unsigned int i = 0; i < feSpaces.size(); i++)
+      for (unsigned int i = 0; i < feSpaces.size(); i++) {
 	addFeSpace(feSpaces[i]);
+	data[feSpaces[i]].setNeedGlobalMapping(needGlobalMapping);
+      }
     }
 
     void update()
     {
+//       for (unsigned int i = 0; i < feSpaces.size(); i++)
+// 	data[feSpaces[i]].update();
+
       nRankDofs = getRankDofs(feSpaces);
       nOverallDofs = getOverallDofs(feSpaces);
       rStartDofs = getStartDofs(feSpaces);
@@ -258,7 +277,7 @@ namespace AMDiS {
       int result = 0;
       for (int i = 0; i < ithFeSpace; i++)
 	result += data[feSpaces[i]].nRankDofs;
-      result += data[feSpaces[ithFeSpace]][index];
+      result += data[feSpaces[ithFeSpace]][index].local;
 
       return result;
     }
@@ -277,7 +296,7 @@ namespace AMDiS {
       int result = rStartDofs;
       for (int i = 0; i < ithFeSpace; i++)
 	result += data[feSpaces[i]].nRankDofs;
-      result += data[feSpaces[ithFeSpace]][index];
+      result += data[feSpaces[ithFeSpace]][index].local;
 
       return result;
     }

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -224,11 +224,76 @@ namespace AMDiS {
       createLagrange(feSpace);      
       createIndexB(feSpace);
     }
+
+    primalDofMap.update();
+    dualDofMap.update();
+    lagrangeMap.update();
+    localDofMap.update();
+
+    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
+      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
+
+      MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
+	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
+      
+      MSG("nRankDuals = %d   nOverallDuals = %d\n",
+	  dualDofMap[feSpace].nRankDofs, 
+	  dualDofMap[feSpace].nOverallDofs);
+
+      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
+	  lagrangeMap[feSpace].nRankDofs, 
+	  lagrangeMap[feSpace].nOverallDofs);
+    }
+
+    // === Communicate lagrangeMap to all other ranks. ===
+
+    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
+      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
+
+      StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
+      
+      for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
+	   !it.end(); it.nextRank()) {
+	for (; !it.endDofIter(); it.nextDof())
+	  if (!isPrimal(feSpace, it.getDofIndex())) {
+	    DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()].local;
+	    stdMpi.getSendData(it.getRank()).push_back(d);
+	  }
+      }
+      
+      stdMpi.updateSendDataSize();
+      
+      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
+	   !it.end(); it.nextRank()) {
+	bool recvData = false;
+	for (; !it.endDofIter(); it.nextDof())
+	  if (!isPrimal(feSpace, it.getDofIndex())) {
+	    recvData = true;
+	    break;
+	  }
+	
+	if (recvData)
+	  stdMpi.recv(it.getRank());
+      }
+      
+      stdMpi.startCommunication();
+
+      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
+	   !it.end(); it.nextRank()) {
+	int counter = 0;
+	for (; !it.endDofIter(); it.nextDof()) {
+	  if (!isPrimal(feSpace, it.getDofIndex())) {
+	    DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
+	    lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
+	  }
+	}
+      }
+    }
   }
 
 
   void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
-  {
+ {
     FUNCNAME("PetscSolverFeti::createPrimals()");  
 
     // === Define all vertices on the interior boundaries of the macro mesh ===
@@ -257,14 +322,8 @@ namespace AMDiS {
     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);
-
-    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());
+    primalDofMap[feSpace].update();
   }
 
 
@@ -335,11 +394,7 @@ namespace AMDiS {
       if (!isPrimal(feSpace, **it))
 	dualDofMap[feSpace].insertRankDof(**it);
 
-    dualDofMap[feSpace].update(false);
-
-    MSG("nRankDuals = %d   nOverallDuals = %d\n",
-	dualDofMap[feSpace].nRankDofs, 
-	dualDofMap[feSpace].nOverallDofs);
+    dualDofMap[feSpace].update();
   }
 
   
@@ -351,8 +406,8 @@ namespace AMDiS {
     // === appropriate number of Lagrange constraints.                      ===
 
     int nRankLagrange = 0;
-    DofMapping& dualMap = dualDofMap[feSpace].getMap();
-    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
+    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
       if (meshDistributor->getIsRankDof(feSpace, it->first)) {
 	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
 	int degree = boundaryDofRanks[it->first].size();
@@ -361,52 +416,6 @@ namespace AMDiS {
     }
     lagrangeMap[feSpace].nRankDofs = nRankLagrange;
     lagrangeMap[feSpace].update();
-    
-    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
-	lagrangeMap[feSpace].nRankDofs, 
-	lagrangeMap[feSpace].nOverallDofs);
-
-
-    // === Communicate lagrangeMap to all other ranks. ===
-
-    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
-
-    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
-	 !it.end(); it.nextRank()) {
-      for (; !it.endDofIter(); it.nextDof())
-	if (!isPrimal(feSpace, it.getDofIndex())) {
-	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
-	  stdMpi.getSendData(it.getRank()).push_back(d);
-	}
-    }
-
-    stdMpi.updateSendDataSize();
-
-    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
-	 !it.end(); it.nextRank()) {
-      bool recvData = false;
-      for (; !it.endDofIter(); it.nextDof())
-	if (!isPrimal(feSpace, it.getDofIndex())) {
-	  recvData = true;
-	  break;
-	}
-	  
-      if (recvData)
-	stdMpi.recv(it.getRank());
-    }
-
-    stdMpi.startCommunication();
-
-    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
-	 !it.end(); it.nextRank()) {
-      int counter = 0;
-      for (; !it.endDofIter(); it.nextDof()) {
-	if (!isPrimal(feSpace, it.getDofIndex())) {
-	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
-	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
-	}
-      }
-    }
   }
 
 
@@ -425,13 +434,11 @@ namespace AMDiS {
       if (admin->isDofFree(i) == false && 
 	  isPrimal(feSpace, i) == false &&
 	  dualDofMap[feSpace].isSet(i) == false) {
-	localDofMap[feSpace].insertRankDof(i);
+	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
 	nLocalInterior++;
       }
     }
     
-    localDofMap[feSpace].update(false);
-
     TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
 		  dualDofMap[feSpace].size() == 
 		  static_cast<unsigned int>(admin->getUsedDofs()))
@@ -439,14 +446,11 @@ namespace AMDiS {
 
     // === And finally, add the global indicies of all dual nodes. ===
 
-    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
+    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
 	 it != dualDofMap[feSpace].getMap().end(); ++it)
       localDofMap[feSpace].insertRankDof(it->first);
 
-    localDofMap[feSpace].update(false);
-
-    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
-				  nLocalInterior);
+    localDofMap[feSpace].update();
   }
 
 
@@ -471,13 +475,13 @@ namespace AMDiS {
     // === m == r, than the rank sets -1.0 for the corresponding           ===
     // === constraint.                                                     ===
 
-    DofMapping &dualMap = dualDofMap[feSpace].getMap();
-    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
+    map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();