From a8523376c0f998bcc1c3654b7cb246ecc29a1d2c Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 18 Apr 2012 19:18:26 +0000
Subject: [PATCH] MeshDistributer makes use of parallel DOF mapping.

---
 AMDiS/src/parallel/MeshDistributor.cc         | 133 ++++++------------
 AMDiS/src/parallel/MeshDistributor.h          |  94 +++----------
 AMDiS/src/parallel/MeshPartitioner.h          |   7 +-
 AMDiS/src/parallel/ParMetisPartitioner.cc     |   5 +-
 AMDiS/src/parallel/ParMetisPartitioner.h      |   3 +-
 AMDiS/src/parallel/ParallelDebug.cc           |  31 ++--
 .../parallel/PetscSolverGlobalBlockMatrix.cc  |   5 +-
 AMDiS/src/parallel/PetscSolverGlobalMatrix.cc |   5 +-
 8 files changed, 94 insertions(+), 189 deletions(-)

diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index f88f3860..b4cd1ff9 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -1065,7 +1065,7 @@ namespace AMDiS {
 
   
   void MeshDistributor::deserialize(istream &in, DofContainer &data,
-				    map<int, const DegreeOfFreedom*> &dofMap)
+				    map<int, const DegreeOfFreedom*> &dofIndexMap)
   {
     FUNCNAME("MeshDistributor::deserialize()");
 
@@ -1077,17 +1077,17 @@ namespace AMDiS {
       int dofIndex = 0;
       SerUtil::deserialize(in, dofIndex);
 
-      TEST_EXIT_DBG(dofMap.count(dofIndex) != 0)
+      TEST_EXIT_DBG(dofIndexMap.count(dofIndex) != 0)
 	("Dof index could not be deserialized correctly!\n");
 
-      data[i] = dofMap[dofIndex];
+      data[i] = dofIndexMap[dofIndex];
     }
   }
 
 
   void MeshDistributor::deserialize(istream &in, 
 				    map<int, map<const FiniteElemSpace*, DofContainer> > &data,
-				    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > &dofMap)
+				    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > &dofIndexMap)
   {
     data.clear();
 
@@ -1098,7 +1098,7 @@ namespace AMDiS {
       SerUtil::deserialize(in, rank);
 
       for (unsigned int j = 0; j < feSpaces.size(); j++)
-	deserialize(in, data[rank][feSpaces[j]], dofMap[feSpaces[j]]);
+	deserialize(in, data[rank][feSpaces[j]], dofIndexMap[feSpaces[j]]);
     }
   }
 
@@ -1200,7 +1200,8 @@ namespace AMDiS {
 
     // === Run mesh partitioner to calculate a new mesh partitioning.  ===
 
-    partitioner->setLocalGlobalDofMap(&(dofFeData[feSpaces[0]].mapDofToGlobal));
+    partitioner->setLocalGlobalDofMap(&(dofMap[feSpaces[0]].getMap()));
+
     bool partitioningSucceed = 
       partitioner->partition(elemWeights, ADAPTIVE_REPART);
     if (!partitioningSucceed) {
@@ -1966,29 +1967,30 @@ namespace AMDiS {
     recvDofs.init(nLevels);
     boundaryDofInfo.resize(nLevels);
 
+    dofMap.init(mpiComm, feSpaces, feSpaces, true, true);
+    dofMap.setDofComm(sendDofs, recvDofs);
+
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       updateLocalGlobalNumbering(feSpaces[i]);
 
+    dofMap.update();
 
 #if (DEBUG != 0)
     MSG("------------- Debug information -------------\n");
     for (unsigned int i = 0; i < feSpaces.size(); i++) {
       MSG("FE space %d:\n", i);
-      MSG("   nRankDofs = %d\n", dofFeData[feSpaces[i]].nRankDofs);
-      MSG("   nOverallDofs = %d\n", dofFeData[feSpaces[i]].nOverallDofs);
-      MSG("   rStartDofs = %d\n", dofFeData[feSpaces[i]].rStartDofs);
+      MSG("   nRankDofs = %d\n", dofMap[feSpaces[i]].nRankDofs);
+      MSG("   nOverallDofs = %d\n", dofMap[feSpaces[i]].nOverallDofs);
+      MSG("   rStartDofs = %d\n", dofMap[feSpaces[i]].rStartDofs);
     }
 
     stringstream oss;
     oss << debugOutputDir << "elementIndex-" << mpiRank << ".vtu";
     debug::writeElementIndexMesh(mesh, oss.str());
-
     ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
-
     debug::testSortedDofs(mesh, elMap);
     ParallelDebug::testCommonDofs(*this, true);
-    ParallelDebug::testGlobalIndexByCoords(*this);
-    
+    ParallelDebug::testGlobalIndexByCoords(*this);   
 #else
     int tmp = 0;
     Parameters::get(name + "->write parallel debug file", tmp);
@@ -2004,6 +2006,8 @@ namespace AMDiS {
 
     mesh->dofCompress();
 
+    dofMap.clear();
+
     // === Get all DOFs in ranks partition. ===
 
     std::set<const DegreeOfFreedom*> rankDofSet;
@@ -2031,20 +2035,13 @@ namespace AMDiS {
       }
     }
 	
-    // Get displacment for global rank DOF ordering and global DOF number.
-    dofFeData[feSpace].nRankDofs = rankDofs.size();
-    mpi::getDofNumbering(mpiComm, 
-			 dofFeData[feSpace].nRankDofs, 
-			 dofFeData[feSpace].rStartDofs, 
-			 dofFeData[feSpace].nOverallDofs);
-
-
-    // Stores for all rank owned DOFs a new global index.
-    DofIndexMap rankDofsNewGlobalIndex;
     for (unsigned int i = 0; i < rankDofs.size(); i++)
-      rankDofsNewGlobalIndex[rankDofs[i]] = i + dofFeData[feSpace].rStartDofs;
-
+      dofMap[feSpace].insertRankDof(*(rankDofs[i]));
 
+    for (DofComm::Iterator it(recvDofs, feSpace); !it.end(); it.nextRank())
+      for (; !it.endDofIter(); it.nextDof())
+	dofMap[feSpace].insert(it.getDofIndex());
+    
     // === Send and receive new DOF indices. ===
 
 #if (DEBUG != 0)
@@ -2053,52 +2050,15 @@ namespace AMDiS {
 						 recvDofs.getData());
 #endif
     
-    StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
-    for (DofComm::Iterator it(sendDofs, feSpace); !it.end(); it.nextRank()) {
-      stdMpi.getSendData(it.getRank()).resize(0);
-      stdMpi.getSendData(it.getRank()).reserve(it.getDofs().size());
-
-      for (; !it.endDofIter(); it.nextDof())
-	stdMpi.getSendData(it.getRank()).
-	  push_back(rankDofsNewGlobalIndex[it.getDof()]);
-    }
-
-    stdMpi.updateSendDataSize();
-
-    for (DofComm::Iterator it(recvDofs); !it.end(); it.nextRank())
-      stdMpi.recv(it.getRank());
-
-    stdMpi.startCommunication();
-
-
     // First, we set all DOFs in ranks partition to be owend by the rank. Than,
     // the DOFs in ranks partition that are owned by other rank are set to false.
     dofFeData[feSpace].isRankDof.clear();
     for (int i = 0; i < nRankAllDofs; i++)
       dofFeData[feSpace].isRankDof[i] = true;    
 
-    for (DofComm::Iterator it(recvDofs, feSpace); !it.end(); it.nextRank()) {
-      for (; !it.endDofIter(); it.nextDof()) {
-	rankDofsNewGlobalIndex[it.getDof()] = 
-	  stdMpi.getRecvData(it.getRank())[it.getDofCounter()];
+    for (DofComm::Iterator it(recvDofs, feSpace); !it.end(); it.nextRank())
+      for (; !it.endDofIter(); it.nextDof())
 	dofFeData[feSpace].isRankDof[it.getDofIndex()] = false;
-      }
-    }
-
-
-    // === Create now the local to global index and local to DOF    ===
-    // === index mappings.                                          ===
-
-    dofFeData[feSpace].mapDofToGlobal.clear();
-    dofFeData[feSpace].mapLocalToDof.clear();
-
-    for (DofIndexMap::iterator dofIt = rankDofsNewGlobalIndex.begin();
-	 dofIt != rankDofsNewGlobalIndex.end(); ++dofIt)
-      dofFeData[feSpace].mapDofToGlobal[*(dofIt->first)] = dofIt->second;   
-
-    for (unsigned int i = 0; i < rankDofs.size(); i++)
-      dofFeData[feSpace].mapLocalToDof[i] = *(rankDofs[i]); 
-
 
     // === Update DOF admins due to new number of DOFs. ===
   
@@ -2163,10 +2123,8 @@ namespace AMDiS {
 	  BoundaryType type = bound.type;
 
 	  for (unsigned int j = 0; j < dofs0.size(); j++) {
-	    DegreeOfFreedom globalDof0 = 
-	      dofFeData[feSpace].mapDofToGlobal[*(dofs0[j])];
-	    DegreeOfFreedom globalDof1 = 
-	      dofFeData[feSpace].mapDofToGlobal[*(dofs1[j])];
+	    DegreeOfFreedom globalDof0 = dofMap[feSpace][*(dofs0[j])].global;
+	    DegreeOfFreedom globalDof1 = dofMap[feSpace][*(dofs1[j])].global;
 
 	    if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0))
 	      periodicMap.add(feSpace, type, globalDof0, globalDof1);
@@ -2191,7 +2149,7 @@ namespace AMDiS {
 	// Send the global indices to the rank on the other side.
 	stdMpi.getSendData(it->first).reserve(dofs.size());
 	for (unsigned int i = 0; i < dofs.size(); i++)
-	  stdMpi.getSendData(it->first).push_back(dofFeData[feSpace].mapDofToGlobal[*(dofs[i])]);
+	  stdMpi.getSendData(it->first).push_back(dofMap[feSpace][*(dofs[i])].global);
 	
 	// Receive from this rank the same number of dofs.
 	stdMpi.recv(it->first, dofs.size());
@@ -2217,7 +2175,7 @@ namespace AMDiS {
 
       // Added the received DOFs to the mapping.
       for (unsigned int i = 0; i < dofs.size(); i++) {
-	int globalDofIndex = dofFeData[feSpace].mapDofToGlobal[*(dofs[i])];
+	int globalDofIndex = dofMap[feSpace][*(dofs[i])].global;
 	int mapGlobalDofIndex = stdMpi.getRecvData(it->first)[i];
 	BoundaryType type = types[i];
 
@@ -2250,8 +2208,7 @@ namespace AMDiS {
 	boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs);
 
 	for (unsigned int i = 0; i < dofs.size(); i++) {
-	  DegreeOfFreedom globalDof = 
-	    dofFeData[feSpace].mapDofToGlobal[*dofs[i]];
+	  DegreeOfFreedom globalDof = dofMap[feSpace][*dofs[i]].global;
 
 	  std::set<BoundaryType>& assoc = 
 	    periodicMap.getAssociations(feSpace, globalDof);
@@ -2323,11 +2280,13 @@ namespace AMDiS {
 
     TEST_EXIT_DBG(dofFeData.count(feSpace))("Should not happen!\n");
 
+    ERROR_EXIT("REIMPLEMENT THAT SHIT!\n");
+    /*
     for (DofMap::iterator it = dofFeData[feSpace].mapDofToGlobal.begin();
 	 it != dofFeData[feSpace].mapDofToGlobal.end(); ++it) 
       if (it->second == dof)
 	return it->first;
-
+    */
     return -1;
   }
 
@@ -2358,13 +2317,13 @@ namespace AMDiS {
     SerUtil::serialize(out, nFeSpace);
 
     for (unsigned int i = 0; i < nFeSpace; i++) {
-      SerUtil::serialize(out, dofFeData[feSpaces[i]].nRankDofs);
-      SerUtil::serialize(out, dofFeData[feSpaces[i]].nOverallDofs);
-      SerUtil::serialize(out, dofFeData[feSpaces[i]].rStartDofs);
+      //      SerUtil::serialize(out, dofFeData[feSpaces[i]].nRankDofs);
+      //      SerUtil::serialize(out, dofFeData[feSpaces[i]].nOverallDofs);
+      //      SerUtil::serialize(out, dofFeData[feSpaces[i]].rStartDofs);
       
       SerUtil::serialize(out, dofFeData[feSpaces[i]].isRankDof);
-      SerUtil::serialize(out, dofFeData[feSpaces[i]].mapDofToGlobal);
-      SerUtil::serialize(out, dofFeData[feSpaces[i]].mapLocalToDof);
+      //      SerUtil::serialize(out, dofFeData[feSpaces[i]].mapDofToGlobal);
+      //      SerUtil::serialize(out, dofFeData[feSpaces[i]].mapLocalToDof);
     }
 
     periodicMap.serialize(out, feSpaces);
@@ -2395,7 +2354,7 @@ namespace AMDiS {
     // Create two maps: one from from element indices to the corresponding element 
     // pointers, and one map from Dof indices to the corresponding dof pointers.
     map<int, Element*> elIndexMap;
-    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofMap;
+    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap;
     for (unsigned int i = 0; i < feSpaces.size(); i++) {
       ElementDofIterator elDofIter(feSpaces[i]);
       TraverseStack stack;
@@ -2407,7 +2366,7 @@ namespace AMDiS {
 	if (el->isLeaf()) {
 	  elDofIter.reset(el);
 	  do {
-	    dofMap[feSpaces[i]][elDofIter.getDof()] = elDofIter.getDofPtr();
+	    dofIndexMap[feSpaces[i]][elDofIter.getDof()] = elDofIter.getDofPtr();
 	  } while (elDofIter.next());      
 	}
 	
@@ -2421,8 +2380,8 @@ namespace AMDiS {
     otherIntBoundary.deserialize(in, elIndexMap);
     periodicBoundary.deserialize(in, elIndexMap);
 
-    deserialize(in, sendDofs.getData(), dofMap);
-    deserialize(in, recvDofs.getData(), dofMap);
+    deserialize(in, sendDofs.getData(), dofIndexMap);
+    deserialize(in, recvDofs.getData(), dofIndexMap);
 
     // === Deerialieze FE space dependent data ===
     
@@ -2433,13 +2392,13 @@ namespace AMDiS {
        feSpaces.size(), nFeSpace);
 
     for (unsigned int i = 0; i < nFeSpace; i++) {
-      SerUtil::deserialize(in, dofFeData[feSpaces[i]].nRankDofs);    
-      SerUtil::deserialize(in, dofFeData[feSpaces[i]].nOverallDofs);
-      SerUtil::deserialize(in, dofFeData[feSpaces[i]].rStartDofs);
+//       SerUtil::deserialize(in, dofFeData[feSpaces[i]].nRankDofs);    
+//       SerUtil::deserialize(in, dofFeData[feSpaces[i]].nOverallDofs);
+//       SerUtil::deserialize(in, dofFeData[feSpaces[i]].rStartDofs);
       
       SerUtil::deserialize(in, dofFeData[feSpaces[i]].isRankDof);
-      SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapDofToGlobal);
-      SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapLocalToDof);
+      //      SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapDofToGlobal);
+      //      SerUtil::deserialize(in, dofFeData[feSpaces[i]].mapLocalToDof);
     }
 
     periodicMap.deserialize(in, feSpaces);
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index 27be4794..f1aef6c0 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -31,6 +31,7 @@
 #include "parallel/MeshLevelData.h"
 #include "parallel/MeshPartitioner.h"
 #include "parallel/InteriorBoundary.h"
+#include "parallel/ParallelDofMapping.h"
 #include "parallel/PeriodicMap.h"
 #include "parallel/StdMpi.h"
 #include "AMDiS_fwd.h"
@@ -56,25 +57,10 @@ namespace AMDiS {
 
   struct DofData
   {
-    /// Number of DOFs in the rank mesh.
-    int nRankDofs;
-
-    /// Is the index of the first global DOF index, which is owned by the rank.
-    int rStartDofs;
-
-    /// Number of DOFs in the whole domain.
-    int nOverallDofs;
-
     /// Maps all DOFs in ranks partition to a bool value. If it is true, the DOF 
     /// is owned by the rank. Otherwise, its an interior boundary DOF that is 
     /// owned by another rank.
     DofIndexToBool isRankDof;
-
-    /// Maps local to global dof indices.
-    DofMap mapDofToGlobal;
-
-    /// Maps local dof indices to real dof indices.
-    DofMap mapLocalToDof;
   };
 
 
@@ -171,11 +157,7 @@ namespace AMDiS {
     /// Returns the number of DOFs in rank's domain for a given FE space.
     inline int getNumberRankDofs(const FiniteElemSpace *feSpace) 
     {
-      FUNCNAME("MeshDistributor::getNumberRankDofs()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))("Should not happen!\n");
-
-      return dofFeData[feSpace].nRankDofs;
+      return dofMap[feSpace].nRankDofs;
     }
 
     /// Returns the number of DOFs in rank's domain for a set of FE spaces.
@@ -184,10 +166,8 @@ namespace AMDiS {
       FUNCNAME("MeshDistributor::getNumberRankDofs()");
 
       int result = 0;
-      for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	TEST_EXIT_DBG(dofFeData.count(feSpaces[i]))("Should not happen!\n");
-	result += dofFeData[feSpaces[i]].nRankDofs;
-      }
+      for (unsigned int i = 0; i < feSpaces.size(); i++)
+	result += dofMap[feSpaces[i]].nRankDofs;
 
       return result;
     }
@@ -195,11 +175,7 @@ namespace AMDiS {
     /// Returns the first global DOF index of an FE space, owned by rank.
     inline int getStartDofs(const FiniteElemSpace *feSpace)
     {
-      FUNCNAME("MeshDistributor::getStartDofs()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))("Should not happen!\n");
-
-      return dofFeData[feSpace].rStartDofs;
+      return dofMap[feSpace].rStartDofs;
     }
 
     /// Returns the first global DOF index for a set of FE spaces, owned by rank.
@@ -208,11 +184,8 @@ namespace AMDiS {
       FUNCNAME("MeshDistributor::getStartDofs()");
 
       int result = 0;
-      for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	TEST_EXIT_DBG(dofFeData.count(feSpaces[i]))("Should not happen!\n");
-
-	result += dofFeData[feSpaces[i]].rStartDofs;
-      }
+      for (unsigned int i = 0; i < feSpaces.size(); i++)
+	result += dofMap[feSpaces[i]].rStartDofs;
 
       return result;
     }
@@ -220,11 +193,7 @@ namespace AMDiS {
     ///  Returns the global number of DOFs for a given FE space.
     inline int getNumberOverallDofs(const FiniteElemSpace *feSpace)
     {
-      FUNCNAME("MeshDistributor::getNumberOverallDofs()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))("Should not happen!\n");
-
-      return dofFeData[feSpace].nOverallDofs;
+      return dofMap[feSpace].nOverallDofs;
     }
 
     /// Returns the global number of DOFs for a set of FE spaces.
@@ -233,34 +202,22 @@ namespace AMDiS {
       FUNCNAME("MeshDistributor::getNumberOverallDofs()");
 
       int result = 0;
-      for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	TEST_EXIT_DBG(dofFeData.count(feSpaces[i]))("Should not happen!\n");
-
-	result += dofFeData[feSpaces[i]].nOverallDofs;
-      }
+      for (unsigned int i = 0; i < feSpaces.size(); i++)
+	result += dofMap[feSpaces[i]].nOverallDofs;
 
       return result;
     }
 
-    inline DofMap& getMapDofToGlobal(const FiniteElemSpace *feSpace)
+    inline map<DegreeOfFreedom, MultiIndex>& getMapDofToGlobal(const FiniteElemSpace *feSpace)
     {
-      FUNCNAME("MeshDistributor::getMapDofToGlobal()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))("Should not happen!\n");
-
-      return dofFeData[feSpace].mapDofToGlobal;
+      return dofMap[feSpace].getMap();
     }
 
     /// Maps a local DOF to its global index.
     inline DegreeOfFreedom mapDofToGlobal(const FiniteElemSpace *feSpace,
 					  DegreeOfFreedom dof)
     {
-      FUNCNAME("MeshDistributor::mapDofToGlobal()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))
-	("No DOF data for FE space at addr %p!\n", feSpace);
-
-      return dofFeData[feSpace].mapDofToGlobal[dof];
+      return dofMap[feSpace][dof].global;
     }
 
     /// Returns for a global index the DOF index in rank's subdomain. As there
@@ -270,18 +227,6 @@ namespace AMDiS {
     DegreeOfFreedom mapGlobalToLocal(const FiniteElemSpace *feSpace,
 				     DegreeOfFreedom dof);
 
-    /// Maps a local DOF to its local index.
-    inline DegreeOfFreedom mapLocalToDof(const FiniteElemSpace *feSpace,
-					 DegreeOfFreedom dof)
-    {
-      FUNCNAME("MeshDistributor::mapLocalToDof()");
-
-      TEST_EXIT_DBG(dofFeData.count(feSpace))
-	("No DOF data for FE space at addr %p!\n", feSpace);
-
-      return dofFeData[feSpace].mapLocalToDof[dof];
-    }
-
     /// Returns the periodic mapping handler, \ref periodicMap.
     inline PeriodicMap& getPeriodicMap()
     {
@@ -527,12 +472,12 @@ namespace AMDiS {
 
     /// Reads a vector of dof pointers from an input stream.
     void deserialize(istream &in, DofContainer &data,
-		     map<int, const DegreeOfFreedom*> &dofMap);
+		     map<int, const DegreeOfFreedom*> &dofIndexMap);
 
     /// Reads a \ref RankToDofContainer from an input stream.
     void deserialize(istream &in, 
 		     map<int, map<const FiniteElemSpace*, DofContainer> > &data,
-		     map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > &dofMap);
+		     map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > &dofIndexMap);
 
     /// Writes a mapping from dof pointers to some values to an output stream.
     template<typename T>
@@ -554,7 +499,7 @@ namespace AMDiS {
     /// Reads a mapping from dof pointer to some values from an input stream.
     template<typename T>
     void deserialize(istream &in, map<const DegreeOfFreedom*, T> &data,
-		     map<int, const DegreeOfFreedom*> &dofMap)
+		     map<int, const DegreeOfFreedom*> &dofIndexMap)
     {
       FUNCNAME("ParallelDomainBase::deserialize()");
 
@@ -566,9 +511,10 @@ namespace AMDiS {
 	SerUtil::deserialize(in, v1);
 	SerUtil::deserialize(in, v2);
 
-	TEST_EXIT_DBG(dofMap.count(v1) != 0)("Cannot find DOF %d in map!\n", v1);
+	TEST_EXIT_DBG(dofIndexMap.count(v1) != 0)
+	  ("Cannot find DOF %d in map!\n", v1);
 
-	data[dofMap[v1]] = v2;
+	data[dofIndexMap[v1]] = v2;
       }
     }
 
@@ -621,6 +567,8 @@ namespace AMDiS {
 
     map<const FiniteElemSpace*, DofData> dofFeData;
 
+    ParallelDofMapping dofMap;
+
     /// Database to store and query all sub-objects of all elements of the 
     /// macro mesh.
     ElementObjectDatabase elObjDb;
diff --git a/AMDiS/src/parallel/MeshPartitioner.h b/AMDiS/src/parallel/MeshPartitioner.h
index 1131bed2..9e20ddef 100644
--- a/AMDiS/src/parallel/MeshPartitioner.h
+++ b/AMDiS/src/parallel/MeshPartitioner.h
@@ -30,6 +30,7 @@
 #include "AMDiS_fwd.h"
 #include "Mesh.h"
 #include "parallel/MpiHelper.h"
+#include "parallel/ParallelDofMapping.h"
 
 
 namespace AMDiS {
@@ -106,7 +107,11 @@ namespace AMDiS {
       boxPartitioning = b;
     }
 
+#if 0
     void setLocalGlobalDofMap(map<DegreeOfFreedom, DegreeOfFreedom> *m)
+#else
+      void setLocalGlobalDofMap(map<DegreeOfFreedom, MultiIndex> *m)
+#endif
     {
       mapLocalGlobal = m;
     }
@@ -168,7 +173,7 @@ namespace AMDiS {
     /// macro element index the box number it belongs to.
     map<int, int> elInBox;
 
-    map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal;
+    map<DegreeOfFreedom, MultiIndex> *mapLocalGlobal;
 
     map<int, vector<int> > elNeighbours;
 
diff --git a/AMDiS/src/parallel/ParMetisPartitioner.cc b/AMDiS/src/parallel/ParMetisPartitioner.cc
index 2116c755..4c7e2b8c 100644
--- a/AMDiS/src/parallel/ParMetisPartitioner.cc
+++ b/AMDiS/src/parallel/ParMetisPartitioner.cc
@@ -14,6 +14,7 @@
 #include <mpi.h>
 
 #include "parallel/ParMetisPartitioner.h"
+#include "parallel/ParallelDofMapping.h"
 #include "parallel/MpiHelper.h"
 #include "Serializer.h"
 #include "Mesh.h"
@@ -28,7 +29,7 @@ namespace AMDiS {
 
   ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm, 
 			     std::map<int, bool>& elementInRank,
-			     map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal)
+			     map<DegreeOfFreedom, MultiIndex> *mapLocalGlobal)
     : dim(mesh->getDim()),
       nElements(0),
       mpiComm(comm)
@@ -102,7 +103,7 @@ namespace AMDiS {
 	// write eind entries (element nodes)
 	for (int i = 0; i < dim + 1; i++) {
 	  if (mapLocalGlobal)
-	    *ptr_eind = (*mapLocalGlobal)[element->getDof(i, 0)];
+	    *ptr_eind = (*mapLocalGlobal)[element->getDof(i, 0)].global;
 	  else
 	    *ptr_eind = element->getDof(i, 0);	  
 
diff --git a/AMDiS/src/parallel/ParMetisPartitioner.h b/AMDiS/src/parallel/ParMetisPartitioner.h
index cb1ef465..c95798a5 100644
--- a/AMDiS/src/parallel/ParMetisPartitioner.h
+++ b/AMDiS/src/parallel/ParMetisPartitioner.h
@@ -31,6 +31,7 @@
 #include "AMDiS_fwd.h"
 #include "Global.h"
 #include "parallel/MeshPartitioner.h"
+#include "parallel/ParallelDofMapping.h"
 
 namespace AMDiS {
 
@@ -43,7 +44,7 @@ namespace AMDiS {
   public:
     ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm, 
 		 map<int, bool>& elementInRank,
-		 map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal);
+		 map<DegreeOfFreedom, MultiIndex> *mapLocalGlobal);
 
     ~ParMetisMesh();
 
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index 427599af..ece764f4 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -165,11 +165,6 @@ namespace AMDiS {
       WorldVector<double> c;
       pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
       int nAssoc = it->second.size();
-#if 0
-      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7))
-	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
-	 it->first, c[0], c[1], (pdb.mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
-#endif
     }    
 
 
@@ -484,9 +479,7 @@ namespace AMDiS {
 
     DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
     for (it.reset(); !it.end(); ++it) {
-      coordsToIndex[(*it)] = 
-	pdb.dofFeData[feSpace].mapDofToGlobal[it.getDOFIndex()];    
-
+      coordsToIndex[(*it)] = pdb.dofMap[feSpace][it.getDOFIndex()].global;
 //       MSG("   CHECK FOR DOF %d AT COORDS %f %f %f\n",
 // 	  coordsToIndex[(*it)], (*it)[0], (*it)[1], (*it)[2]);
     }
@@ -646,16 +639,11 @@ namespace AMDiS {
       const FiniteElemSpace *feSpace = pdb.feSpaces[0];
 
       cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
-      
-      for (DofMap::iterator it = pdb.dofFeData[feSpace].mapDofToGlobal.begin();
-	   it != pdb.dofFeData[feSpace].mapDofToGlobal.end(); it++) {
-	DegreeOfFreedom localdof = -1;
-	if (pdb.dofFeData[feSpace].mapLocalToDof.count(it->first) > 0)
-	  localdof = pdb.dofFeData[feSpace].mapLocalToDof[it->first];
-	
-	cout << "DOF " << it->first << " " 
-		  << it->second << " " 
-		  << localdof << endl;
+            
+      map<DegreeOfFreedom, MultiIndex> &dofMap = pdb.dofMap[feSpace].getMap();
+      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
+	   it != dofMap.end(); ++it) {
+	cout << "DOF " << it->first << " " << it->second.global << "\n";
 	WorldVector<double> coords;
 	pdb.mesh->getDofIndexCoords(it->first, feSpace, coords);
 	coords.print();
@@ -702,8 +690,9 @@ namespace AMDiS {
 	cout << endl;
 
 	DegreeOfFreedom localdof = -1;
-	for (DofMap::iterator dofIt = pdb.dofFeData[feSpace].mapDofToGlobal.begin();
-	     dofIt != pdb.dofFeData[feSpace].mapDofToGlobal.end(); ++dofIt)
+	map<DegreeOfFreedom, MultiIndex> &dofMap = pdb.dofMap[feSpace].getMap();
+	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
+	     it != dofMap.end(); ++it)
 	  if (dofIt->second == it->first)
 	    localdof = dofIt->first;
 
@@ -819,7 +808,7 @@ namespace AMDiS {
     DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
     for (it.reset(); !it.end(); ++it) {
       file << it.getDOFIndex() << " " 
-	   << pdb.dofFeData[feSpace].mapDofToGlobal[it.getDOFIndex()] << " "
+	   << pdb.dofMap[feSpace][it.getDOFIndex()].global << " "
 	   << pdb.getIsRankDof(feSpace, it.getDOFIndex());
       for (int i = 0; i < pdb.mesh->getDim(); i++)
 	file << " " << (*it)[i];
diff --git a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
index 26771da7..410c420c 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
@@ -154,8 +154,9 @@ namespace AMDiS {
       PetscScalar *vecPointer;
       VecGetArray(tmp, &vecPointer);
 
-      for (int j = 0; j < nRankDofs; j++)
-	dofvec[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[j]; 
+      ERROR_EXIT("REWRITE CODE!\n");
+      //      for (int j = 0; j < nRankDofs; j++)
+      //	dofvec[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[j]; 
 
       VecRestoreArray(tmp, &vecPointer);
     }
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index f9949953..9a963d28 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -197,8 +197,9 @@ namespace AMDiS {
       const FiniteElemSpace *feSpace = dv.getFeSpace();
       int nRankDofs = meshDistributor->getNumberRankDofs(feSpace);
 
-      for (int j = 0; j < nRankDofs; j++)
-	dv[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[c++]; 
+      ERROR_EXIT("REWRITE CODE!\n");
+      //      for (int j = 0; j < nRankDofs; j++)
+      //	dv[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[c++]; 
     }
 
     VecRestoreArray(petscSolVec, &vecPointer);
-- 
GitLab