From 3056a9f6c6ee133d7a82c925e5ee622a7f9796cd Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 9 Feb 2012 09:36:04 +0000
Subject: [PATCH] Fixed problem for repartitioning with periodic boundary
 conditions.

---
 AMDiS/src/parallel/ElementObjectData.h |  55 +++++------
 AMDiS/src/parallel/MeshDistributor.cc  |  43 +++++----
 AMDiS/src/parallel/MeshPartitioner.h   |  16 +++-
 AMDiS/src/parallel/PetscSolverFeti.cc  | 124 +++++++++++--------------
 AMDiS/src/parallel/PetscSolverFeti.h   |  12 +--
 5 files changed, 128 insertions(+), 122 deletions(-)

diff --git a/AMDiS/src/parallel/ElementObjectData.h b/AMDiS/src/parallel/ElementObjectData.h
index 42123d2b..cb0d444b 100644
--- a/AMDiS/src/parallel/ElementObjectData.h
+++ b/AMDiS/src/parallel/ElementObjectData.h
@@ -90,13 +90,13 @@ namespace AMDiS {
 
   /** \brief
    * This class is a database of element objects. An element object is either a
-   * vertex, edge or the face of a specific element. This database is used to store
-   * all objects of all elements of a mesh. The information is stored in a way that
-   * makes it possible to identify all elements, which have a given vertex, edge or
-   * face in common. If is is known which element is owned by which rank in parallel
-   * computations, it is thus possible to get all interior boundaries on object 
-   * level. This is required, because two elements may share a common vertex without
-   * beging neighbours in the definition of AMDiS.
+   * vertex, edge or the face of a specific element. This database is used to
+   * store all objects of all elements of a mesh. The information is stored in a
+   * way that makes it possible to identify all elements, which have a given
+   * vertex, edge or face in common. If is is known which element is owned by 
+   * which rank in parallel computations, it is thus possible to get all interior
+   * boundaries on object level. This is required, because two elements may share
+   * a common vertex without beging neighbours in the definition of AMDiS.
    */
   class ElementObjects {
   public:
@@ -114,9 +114,9 @@ namespace AMDiS {
 
 
     /** \brief
-     * Adds an element to the object database. If the element is part of a periodic
-     * boundary, all information about subobjects of the element on this boundary
-     * are collected.
+     * Adds an element to the object database. If the element is part of a
+     * periodic boundary, all information about subobjects of the element on
+     * this boundary are collected.
      *
      * \param[in]  elInfo    ElInfo object of the element. 
      */
@@ -136,15 +136,16 @@ namespace AMDiS {
 
 
     /** \brief
-     * Create for a filled object database the membership information for all element
-     * objects. An object is owned by a rank, if the rank has the heighest rank
-     * number of all ranks where the object is part of.
+     * Create for a filled object database the membership information for all
+     * element objects. An object is owned by a rank, if the rank has the
+     * heighest rank number of all ranks where the object is part of.
      *
-     * \param[in]  macroElementRankMap   Maps to each macro element of the mesh the
-     *                                   rank that owns this macro element.
+     * \param[in]  macroElementRankMap   Maps to each macro element of the mesh
+     *                                   the rank that owns this macro element.
      */
     void createRankData(map<int, int>& macroElementRankMap);
 
+
     /** \brief
      * Creates on all boundaries the reverse mode flag.
      *
@@ -162,12 +163,12 @@ namespace AMDiS {
 
 
     /** \brief
-     * Iterates over all elements for one geometrical index, i.e., over all vertices,
-     * edges or faces in the mesh. The function returns true, if the result is valid.
-     * Otherwise the iterator is at the end position.
+     * Iterates over all elements for one geometrical index, i.e., over all
+     * vertices, edges or faces in the mesh. The function returns true, if the
+     * result is valid. Otherwise the iterator is at the end position.
      *
-     * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the elements
-     *                   that should be traversed.
+     * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the
+     *                   elements that should be traversed.
      */
     bool iterate(GeoIndex pos)
     {
@@ -419,7 +420,8 @@ namespace AMDiS {
       return periodicFaces;
     }
 
-    inline bool getEdgeReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    inline bool getEdgeReverseMode(ElementObjectData &obj0, 
+				   ElementObjectData &obj1)
     {
       if (mesh->getDim() == 2)
 	return true;
@@ -430,7 +432,8 @@ namespace AMDiS {
       return edgeReverseMode[make_pair(obj0, obj1)];
     }
 
-    inline bool getFaceReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    inline bool getFaceReverseMode(ElementObjectData &obj0, 
+				   ElementObjectData &obj1)
     {
       TEST_EXIT_DBG(faceReverseMode.count(make_pair(obj0, obj1)))
 	("Should not happen!\n");
@@ -540,12 +543,12 @@ namespace AMDiS {
     /// objects that have this vertex DOF in common.
     map<DegreeOfFreedom, map<int, ElementObjectData> > vertexInRank;
 
-    /// Defines to each edge a map that maps to each rank number the element objects
-    /// that have this edge in common.
+    /// Defines to each edge a map that maps to each rank number the element 
+    /// objects that have this edge in common.
     map<DofEdge, map<int, ElementObjectData> > edgeInRank;
 
-    /// Defines to each face a map that maps to each rank number the element objects
-    /// that have this face in common.
+    /// Defines to each face a map that maps to each rank number the element 
+    /// objects that have this face in common.
     map<DofFace, map<int, ElementObjectData> > faceInRank;
 
 
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index d36930c1..e38fe776 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -94,16 +94,10 @@ namespace AMDiS {
     mpiSize = MPI::COMM_WORLD.Get_size();
     mpiComm = MPI::COMM_WORLD;
 
-    int tmp = 0;
-    Parameters::get(name + "->repartitioning", tmp);
-    repartitioningAllowed = (tmp > 0);
-
+    Parameters::get(name + "->repartitioning", repartitioningAllowed);
     Parameters::get(name + "->debug output dir", debugOutputDir);
     Parameters::get(name + "->repartition ith change", repartitionIthChange);
-
-    tmp = 0;
-    Parameters::get(name + "->log main rank", tmp);
-    Msg::outputMainRank = (tmp > 0);
+    Parameters::get(name + "->log main rank", Msg::outputMainRank);
 
     string partStr = "parmetis";
     Parameters::get(name + "->partitioner", partStr);
@@ -125,7 +119,7 @@ namespace AMDiS {
     if (partStr == "simple")
       partitioner = new SimplePartitioner(&mpiComm);
 
-    tmp = 0;
+    int tmp = 0;
     Parameters::get(name + "->box partitioning", tmp);
     partitioner->setBoxPartitioning(static_cast<bool>(tmp));
 
@@ -874,7 +868,7 @@ namespace AMDiS {
       repartitionMesh();
       nMeshChangesAfterLastRepartitioning = 0;
     } else {
-      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   reparttitioningAllowed = %d   nMeshChange =%d    repartitionIthChange = %d\n",
+      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   repartitioningAllowed = %d   nMeshChange = %d   repartitionIthChange = %d\n",
 	      tryRepartition, repartitioningAllowed, 
 	      nMeshChangesAfterLastRepartitioning, repartitionIthChange);
     }
@@ -896,15 +890,18 @@ namespace AMDiS {
     if (mpiRank == 0) {
       int nOverallDofs = 0;
       int maxDofs = numeric_limits<int>::min();
+      int minDofs = numeric_limits<int>::max();
       for (int i = 0; i < mpiSize; i++) {
 	nOverallDofs += nDofsInRank[i];
 	maxDofs = std::max(maxDofs, nDofsInRank[i]);
+	minDofs = std::min(minDofs, nDofsInRank[i]);
       }      
-      int avrgDofs = nOverallDofs / mpiSize;
-      double imbalance = 
-	(static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;
+//       int avrgDofs = nOverallDofs / mpiSize;
+//       double imbalance0 = 
+// 	(static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;
+      double imbalance1 = (static_cast<double>(maxDofs) / minDofs - 1.0) * 100.0;
 
-      MSG("Imbalancing factor: %.1f\%\n", imbalance);
+      MSG("Imbalancing factor: %.1f\%\n", imbalance1);
     }
   }
 
@@ -1157,12 +1154,18 @@ namespace AMDiS {
     // === Run mesh partitioner to calculate a new mesh partitioning.  ===
 
     partitioner->setLocalGlobalDofMap(&(dofFeData[feSpaces[0]].mapDofToGlobal));
-    bool partitioningSucceed = partitioner->partition(elemWeights, ADAPTIVE_REPART);
+    bool partitioningSucceed = 
+      partitioner->partition(elemWeights, ADAPTIVE_REPART);
     if (!partitioningSucceed) {
       MSG("Mesh partitioner created empty partition!\n");
       return;
     }
 
+    // In the case the partitioner does not create a new mesh partition, return
+    // without and changes.
+    if (!partitioner->meshChanged())
+      return;
+
 
     // === Create map that maps macro element indices to pointers to the ===
     // === macro elements.                                               ===
@@ -2023,13 +2026,17 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createPeriodicMap()");
 
-    if (periodicBoundary.boundary.size() == 0)
-      return;
-
     // Clear all periodic DOF mappings calculated before. We do it from scratch.
     periodicDofs.clear();
     periodicMap.clear();
 
+    // If there are no periodic boundaries, return. Note that periodicDofs and
+    // periodicMap must be still cleared before: if we do repartitioning and
+    // there were periodic boundaries in subdomain before and after repartitioning
+    // there are no more periodic boundaries.
+    if (periodicBoundary.boundary.size() == 0)
+      return;
+
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       createPeriodicMap(feSpaces[i]);
   }
diff --git a/AMDiS/src/parallel/MeshPartitioner.h b/AMDiS/src/parallel/MeshPartitioner.h
index e89ed1a6..1131bed2 100644
--- a/AMDiS/src/parallel/MeshPartitioner.h
+++ b/AMDiS/src/parallel/MeshPartitioner.h
@@ -29,13 +29,13 @@
 
 #include "AMDiS_fwd.h"
 #include "Mesh.h"
+#include "parallel/MpiHelper.h"
 
-namespace AMDiS {
 
+namespace AMDiS {
 
   using namespace std;
 
-
   enum PartitionMode {
     INITIAL = 0,          // initial partitioning of a unpartitioned mesh 
     ADAPTIVE_REPART = 1,  // repartitioning of a adaptively refined mesh
@@ -130,6 +130,16 @@ namespace AMDiS {
     {
       return sendElements;
     }
+    
+    /// After mesh repartition this function returns true if the mesh must be 
+    /// redistributed on at least one rank.
+    bool meshChanged()
+    {
+      int nChanges = recvElements.size() + sendElements.size();
+      mpi::globalAdd(nChanges);
+
+      return static_cast<bool>(nChanges);
+    }
 
   protected:
     /// Pointer to the MPI communicator the mesh partitioner should make use of.
@@ -170,6 +180,8 @@ namespace AMDiS {
     /// partitiong mode) the rank number the element belongs to.
     map<int, int> partitionMap;
 
+    /// After mesh repartitioning these maps stores which elements are communicated
+    /// from this rank to other ranks.
     map<int, vector<int> > recvElements, sendElements;
   };
 }
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index cd451918..853a96d1 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -398,6 +398,8 @@ namespace AMDiS {
     MSG("nRankDuals = %d   nOverallDuals = %d\n",
 	dualDofMap[feSpace].nRankDofs, 
 	dualDofMap[feSpace].nOverallDofs);
+
+    nLocalDuals = dualDofMap[feSpace].size();
   }
 
   
@@ -409,26 +411,26 @@ namespace AMDiS {
     // === appropriate number of Lagrange constraints.                      ===
 
     const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-    dofFirstLagrange.addFeSpace(feSpace);
+    lagrangeMap.addFeSpace(feSpace);
 
     int nRankLagrange = 0;
     DofMapping& dualMap = dualDofMap[feSpace].getMap();
     for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
       if (meshDistributor->getIsRankDof(feSpace, it->first)) {
-	dofFirstLagrange[feSpace].insert(it->first, nRankLagrange);
+	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
 	int degree = boundaryDofRanks[it->first].size();
 	nRankLagrange += (degree * (degree - 1)) / 2;
       }
     }
-    dofFirstLagrange[feSpace].nRankDofs = nRankLagrange;
-    dofFirstLagrange[feSpace].update();
+    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
+    lagrangeMap[feSpace].update();
     
     MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
-	dofFirstLagrange[feSpace].nRankDofs, 
-	dofFirstLagrange[feSpace].nOverallDofs);
+	lagrangeMap[feSpace].nRankDofs, 
+	lagrangeMap[feSpace].nOverallDofs);
 
 
-    // === Communicate dofFirstLagrange to all other ranks. ===
+    // === Communicate lagrangeMap to all other ranks. ===
 
     StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
 
@@ -436,9 +438,7 @@ namespace AMDiS {
 	 !it.end(); it.nextRank()) {
       for (; !it.endDofIter(); it.nextDof())
 	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
-	  TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it.getDofIndex()))
-	    ("Should not happen!\n");
-	  DegreeOfFreedom d = dofFirstLagrange[feSpace][it.getDofIndex()];
+	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
 	  stdMpi.getSendData(it.getRank()).push_back(d);
 	}
     }
@@ -466,7 +466,7 @@ namespace AMDiS {
       for (; !it.endDofIter(); it.nextDof()) {
 	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
 	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
-	  dofFirstLagrange[feSpace].insert(it.getDofIndex(), d); 
+	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
 	}
       }
     }
@@ -475,36 +475,29 @@ namespace AMDiS {
 
   void PetscSolverFeti::createIndexB()
   {
-    FUNCNAME("PetscSolverFeti::createIndeB()");
+    FUNCNAME("PetscSolverFeti::createIndexB()");
 
     const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-    localIndexB.clear();
+    localDofMap.addFeSpace(feSpace);
     DOFAdmin* admin = feSpace->getAdmin();
 
     // === To ensure that all interior node on each rank are listen first in ===
     // === the global index of all B nodes, insert all interior nodes first, ===
     // === without defining a correct index.                                 ===
-    
-    for (int i = 0; i < admin->getUsedSize(); i++)
+
+    nLocalInterior = 0;    
+    for (int i = 0; i < admin->getUsedSize(); i++) {
       if (admin->isDofFree(i) == false && 
 	  primalDofMap[feSpace].isSet(i) == false && 
-	  dualDofMap[feSpace].isSet(i) == false)
-	localIndexB[i] = -1;
-
-
-    // === Get correct index for all interior nodes. ===
-
-    nLocalInterior = 0;
-    for (DofMapping::iterator it = localIndexB.begin(); 
-	 it != localIndexB.end(); ++it) {
-      it->second = nLocalInterior;
-      nLocalInterior++;
+	  dualDofMap[feSpace].isSet(i) == false) {
+	localDofMap[feSpace].insertRankDof(i);
+	nLocalInterior++;
+      }
     }
-    nLocalDuals = dualDofMap[feSpace].size();
 
     TEST_EXIT_DBG(nLocalInterior + 
 		  primalDofMap[feSpace].size() + 
-		  dualDofMap[feSpace].size() == 
+		  nLocalDuals == 
 		  static_cast<unsigned int>(admin->getUsedDofs()))
       ("Should not happen!\n");
 
@@ -513,7 +506,8 @@ namespace AMDiS {
 
     for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
 	 it != dualDofMap[feSpace].getMap().end(); ++it)
-      localIndexB[it->first] = it->second - rStartB;
+      localDofMap[feSpace].insert(it->first, it->second - rStartB);
+      //      localDofMap[feSpace].insertRankDof(it->first);
   }
 
 
@@ -526,9 +520,9 @@ namespace AMDiS {
     // === Create distributed matrix for Lagrange constraints. ===
 
     MatCreateMPIAIJ(PETSC_COMM_WORLD,
-		    dofFirstLagrange[feSpace].nRankDofs * nComponents, 
+		    lagrangeMap[feSpace].nRankDofs * nComponents, 
 		    nRankB * nComponents,
-		    dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
+		    lagrangeMap[feSpace].nOverallDofs * nComponents, 
 		    nOverallB * nComponents,
 		    2, PETSC_NULL, 2, PETSC_NULL,
 		    &mat_lagrange);
@@ -542,13 +536,11 @@ namespace AMDiS {
 
     DofMapping &dualMap = dualDofMap[feSpace].getMap();
     for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
-      TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it->first))
-	("Should not happen!\n");
       TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
 	("Should not happen!\n");
 
       // Global index of the first Lagrange constriant for this node.
-      int index = dofFirstLagrange[feSpace][it->first];
+      int index = lagrangeMap[feSpace][it->first];
       // Copy set of all ranks that contain this dual node.
       vector<int> W(boundaryDofRanks[it->first].begin(), 
 		    boundaryDofRanks[it->first].end());
@@ -742,18 +734,18 @@ namespace AMDiS {
 		 nRankB * nComponents, nOverallB * nComponents,
 		 &(fetiData.tmp_vec_b));
     VecCreateMPI(PETSC_COMM_WORLD,
-		 dofFirstLagrange[feSpace].nRankDofs * nComponents, 
-		 dofFirstLagrange[feSpace].nOverallDofs * nComponents,
+		 lagrangeMap[feSpace].nRankDofs * nComponents, 
+		 lagrangeMap[feSpace].nOverallDofs * nComponents,
 		 &(fetiData.tmp_vec_lagrange));
     VecCreateMPI(PETSC_COMM_WORLD, 
 		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
 		 &(fetiData.tmp_vec_primal));
 
     MatCreateShell(PETSC_COMM_WORLD,
-		   dofFirstLagrange[feSpace].nRankDofs * nComponents, 
-		   dofFirstLagrange[feSpace].nRankDofs * nComponents,
-		   dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
-		   dofFirstLagrange[feSpace].nOverallDofs * nComponents,
+		   lagrangeMap[feSpace].nRankDofs * nComponents, 
+		   lagrangeMap[feSpace].nRankDofs * nComponents,
+		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
+		   lagrangeMap[feSpace].nOverallDofs * nComponents,
 		   &fetiData, &mat_feti);
     MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
 
@@ -942,8 +934,8 @@ namespace AMDiS {
     for (int i = 0; i < nComponents; i++) {
       DOFVector<double>& dofVec = *(vec.getDOFVector(i));
 
-      for (DofMapping::iterator it = localIndexB.begin();
-	   it != localIndexB.end(); ++it) {
+      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
+	   it != localDofMap[feSpace].getMap().end(); ++it) {
 	int petscIndex = it->second * nComponents + i;
 	dofVec[it->first] = localSolB[petscIndex];
       }
@@ -973,7 +965,8 @@ namespace AMDiS {
 
     primalDofMap.setMpiComm(mpiComm);
     dualDofMap.setMpiComm(mpiComm);
-    dofFirstLagrange.setMpiComm(mpiComm);
+    lagrangeMap.setMpiComm(mpiComm);
+    localDofMap.setMpiComm(mpiComm);
     updateDofData();
 
     // === Create matrices for the FETI-DP method. ===
@@ -1100,11 +1093,8 @@ namespace AMDiS {
 	    } else {
 	      // Column is not a primal variable.
 
-	      TEST_EXIT_DBG(localIndexB.count(col(*icursor)))
-		("No global B index for DOF %d!\n", col(*icursor));
-	      
 	      int colIndex = 
-		(localIndexB[col(*icursor)] + rStartB) * nComponents + j;
+		(localDofMap[feSpace][col(*icursor)] + rStartB) * nComponents + j;
 
 	      if (rowPrimal) {
 		colsOther.push_back(colIndex);
@@ -1120,17 +1110,17 @@ namespace AMDiS {
 	    // === For preconditioner ===
 
 	    if (!rowPrimal && !colPrimal) {
-	      int rowIndex = localIndexB[*cursor];
-	      int colIndex = localIndexB[col(*icursor)];
+	      int rowIndex = localDofMap[feSpace][*cursor];
+	      int colIndex = localDofMap[feSpace][col(*icursor)];
 		
 	      if (rowIndex < nLocalInterior) {
 		if (colIndex < nLocalInterior) {
-		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
+		  int colIndex2 = localDofMap[feSpace][col(*icursor)] * nComponents + j;
 
 		  colsLocal.push_back(colIndex2);
 		  valuesLocal.push_back(value(*icursor));
 		} else {
-		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
+		  int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) *
 		    nComponents + j;
 
 		  colsLocalOther.push_back(colIndex2);
@@ -1138,12 +1128,12 @@ namespace AMDiS {
 		}
 	      } else {
 		if (colIndex < nLocalInterior) {
-		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
+		  int colIndex2 = localDofMap[feSpace][col(*icursor)] * nComponents + j;
 
 		  colsLocalOther.push_back(colIndex2);
 		  valuesLocalOther.push_back(value(*icursor));
 		} else {
-		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
+		  int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) *
 		    nComponents + j;
 
 		  colsLocal.push_back(colIndex2);
@@ -1165,10 +1155,7 @@ namespace AMDiS {
 	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	  } else {
-	    TEST_EXIT_DBG(localIndexB.count(*cursor))
-	      ("Should not happen!\n");
-
-	    int localRowIndex = localIndexB[*cursor] * nComponents + i; 
+	    int localRowIndex = localDofMap[feSpace][*cursor] * nComponents + i; 
 	    for (unsigned int k = 0; k < cols.size(); k++)
 	      cols[k] -= rStartB * nComponents;
 	    MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
@@ -1176,7 +1163,7 @@ namespace AMDiS {
 
 	    if (colsOther.size()) {
 	      int globalRowIndex = 
-		(localIndexB[*cursor] + rStartB) * nComponents + i; 
+		(localDofMap[feSpace][*cursor] + rStartB) * nComponents + i; 
 	      MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	    }
@@ -1188,10 +1175,10 @@ namespace AMDiS {
 	    switch (fetiPreconditioner) {
 	    case FETI_DIRICHLET:
 	      {
-		int rowIndex = localIndexB[*cursor];
+		int rowIndex = localDofMap[feSpace][*cursor];
 		
 		if (rowIndex < nLocalInterior) {
-		  int rowIndex2 = localIndexB[*cursor] * nComponents + i;
+		  int rowIndex2 = localDofMap[feSpace][*cursor] * nComponents + i;
 		  
 		  MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
 			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
@@ -1201,7 +1188,7 @@ namespace AMDiS {
 				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
 		} else {
 		  int rowIndex2 = 
-		    (localIndexB[*cursor] - nLocalInterior) * nComponents + i;
+		    (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i;
 		  
 		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
 			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
@@ -1216,11 +1203,11 @@ namespace AMDiS {
 
 	    case FETI_LUMPED:
 	      {
-		int rowIndex = localIndexB[*cursor];
+		int rowIndex = localDofMap[feSpace][*cursor];
 		
 		if (rowIndex >= nLocalInterior) {
 		  int rowIndex2 = 
-		    (localIndexB[*cursor] - nLocalInterior) * nComponents + i;
+		    (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i;
 		  
 		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
 			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
@@ -1316,10 +1303,7 @@ namespace AMDiS {
 	  double value = *dofIt;
 	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
 	} else {
-	  TEST_EXIT_DBG(localIndexB.count(index))
-	    ("Should not happen!\n");
-
-	  index = (localIndexB[index] + rStartB) * nComponents + i;
+	  index = (localDofMap[feSpace][index] + rStartB) * nComponents + i;
 	  VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
 	}      
       }
@@ -1376,7 +1360,7 @@ namespace AMDiS {
     int nOverallNest = 
       (nOverallB + 
        primalDofMap[feSpace].nOverallDofs + 
-       dofFirstLagrange[feSpace].nOverallDofs) * nComponents;
+       lagrangeMap[feSpace].nOverallDofs) * nComponents;
 
     Mat mat_lagrange_transpose;
     MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
@@ -1460,8 +1444,8 @@ namespace AMDiS {
 
     // === mat_lagrange ===
 
-    for (int i = 0; i < dofFirstLagrange[feSpace].nRankDofs * nComponents; i++) {
-      int rowIndex = dofFirstLagrange[feSpace].rStartDofs * nComponents + i;
+    for (int i = 0; i < lagrangeMap[feSpace].nRankDofs * nComponents; i++) {
+      int rowIndex = lagrangeMap[feSpace].rStartDofs * nComponents + i;
       MatGetRow(mat_lagrange, rowIndex, &nCols, &cols, &vals);
 
       int rowIndexA = (nOverallB + primalDofMap[feSpace].nOverallDofs) * nComponents + rowIndex;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 2c949ff0..15425922 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -156,19 +156,19 @@ namespace AMDiS {
 
     /// Stores to each dual DOF index the index of the first Lagrange
     /// constraint that is assigned to this DOF.
-    FeSpaceData<GlobalDofMap> dofFirstLagrange;
-
-    /// Stores to each dual boundary DOF the set of ranks in which the DOF
-    /// is contained in.
-    DofIndexToPartitions boundaryDofRanks;
+    FeSpaceData<GlobalDofMap> lagrangeMap;
     
     /// Index for each non primal variables to the global index of
     /// B variables.
-    DofMapping globalIndexB, localIndexB;
+    FeSpaceData<GlobalDofMap> localDofMap;
 
     /// Number of non primal, thus B, variables on rank and globally.
     int nRankB, nOverallB, rStartB;    
 
+    /// Stores to each dual boundary DOF the set of ranks in which the DOF
+    /// is contained in.
+    DofIndexToPartitions boundaryDofRanks;
+
     /// Global PETSc matrix of non primal variables.
     Mat mat_b_b;
 
-- 
GitLab