No, does not work, so merge to juropa and go on with debugging.

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -259,31 +259,28 @@ namespace AMDiS {
     }
 
     MSG("FETI-DP data created on mesh level %d\n", meshLevel);
-    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
+    for (unsigned int i = 0; i < meshDistributor->getComponentFeSpaces().size(); i++) {
       const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
-      MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace);
+
+      MSG("FETI-DP data for %d-ith component (FE space %p):\n", i, feSpace);
 
       if (feSpace == pressureFeSpace) {
 	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
-	    interfaceDofMap[feSpace].nRankDofs, 
-	    interfaceDofMap[feSpace].nOverallDofs);
+	    interfaceDofMap[i].nRankDofs, 
+	    interfaceDofMap[i].nOverallDofs);
       } else {
 	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
-	    primalDofMap[feSpace].nRankDofs, 
-	    primalDofMap[feSpace].nLocalDofs,
-	    primalDofMap[feSpace].nOverallDofs);
+	    primalDofMap[i].nRankDofs, 
+	    primalDofMap[i].nLocalDofs,
+	    primalDofMap[i].nOverallDofs);
 	
 	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
-	    dualDofMap[feSpace].nRankDofs, 
-	    dualDofMap[feSpace].nOverallDofs);
+	    dualDofMap[i].nRankDofs, 
+	    dualDofMap[i].nOverallDofs);
 	
 	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
-	    lagrangeMap[feSpace].nRankDofs, 
-	    lagrangeMap[feSpace].nOverallDofs);
-
-// 	TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
-// 		      static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
-// 	  ("Should not happen!\n");	
+	    lagrangeMap[i].nRankDofs, 
+	    lagrangeMap[i].nOverallDofs);
       }
     }
 
@@ -534,7 +531,7 @@ namespace AMDiS {
     // === appropriate number of Lagrange constraints.                      ===
 
     int nRankLagrange = 0;
-    DofMap& dualMap = dualDofMap[feSpace].getMap();
+    DofMap& dualMap = dualDofMap[component].getMap();
     for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
       if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
 	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
@@ -598,8 +595,8 @@ namespace AMDiS {
     
     // === And finally, add the global indicies of all dual nodes. ===
 
-    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
-	 it != dualDofMap[feSpace].getMap().end(); ++it) {
+    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
+	 it != dualDofMap[component].getMap().end(); ++it) {
       if (meshLevel == 0) {
 	localDofMap[component].insertRankDof(it->first);
       } else {
@@ -638,19 +635,19 @@ namespace AMDiS {
     // === m == r, than the rank sets -1.0 for the corresponding           ===
     // === constraint.                                                     ===
 
-    for (unsigned int k = 0; k < feSpaces.size(); k++) {
-      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
+    for (unsigned int component = 0; component < feSpaces.size(); component++) {
+      DofMap &dualMap = dualDofMap[component].getMap();
 
       for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
-	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
+	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[component]].count(it->first))
 	  ("Should not happen!\n");
 	
 	// Global index of the first Lagrange constriant for this node.
-	int index = lagrangeMap.getMatIndex(k, it->first);
+	int index = lagrangeMap.getMatIndex(component, it->first);
 
 	// Copy set of all ranks that contain this dual node.
-	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
-		      boundaryDofRanks[feSpaces[k]][it->first].end());
+	vector<int> W(boundaryDofRanks[feSpaces[component]][it->first].begin(), 
+		      boundaryDofRanks[feSpaces[component]][it->first].end());
 	// Number of ranks that contain this dual node.
 	int degree = W.size();
 
@@ -662,7 +659,7 @@ namespace AMDiS {
 	    if (W[i] == mpiRank || W[j] == mpiRank) {
 	      MatSetValue(mat_lagrange, 
 			  index + counter, 
-			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
+			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
 			  (W[i] == mpiRank ? 1.0 : -1.0),
 			  INSERT_VALUES);
 	    }
@@ -673,7 +670,7 @@ namespace AMDiS {
 	// === Create scaling factors for scaling the lagrange matrix, which ===
 	// === is required for FETI-DP preconditioners.                      ===
 	
-	if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) {
+	if (meshDistributor->getDofMap()[feSpaces[component]].isRankDof(it->first)) {
 	  int nConstraints = (degree * (degree - 1)) / 2;
 	  for (int i = 0; i < nConstraints; i++) {
 	    VecSetValue(vec_scale_lagrange,
@@ -1228,12 +1225,12 @@ namespace AMDiS {
       TEST_EXIT_DBG(meshLevel == 0)
 	("Should not happen, check usage of localDofMap!\n");
 
-      for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
+      for (unsigned int component = 0; component < feSpaces.size(); component++) {
+	DofMap &dualMap = dualDofMap[component].getMap();
 	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
 	  DegreeOfFreedom d = it->first;	  
-	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
-	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
+	  int matIndexLocal = localDofMap.getLocalMatIndex(component, d);
+	  int matIndexDual = dualDofMap.getLocalMatIndex(component, d);
 	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
 	}
       }
@@ -1271,15 +1268,15 @@ namespace AMDiS {
 	MatGetVecs(mat_duals_duals, PETSC_NULL, 
 		   &(lumpedData->tmp_vec_duals1));
 	
-	for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	  if (stokesMode && i == pressureComponent)
+	for (unsigned int component = 0; component < feSpaces.size(); component++) {
+	  if (stokesMode && component == pressureComponent)
 	    continue;
 	  
-	  DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
+	  DofMap &dualMap = dualDofMap[component].getMap();
 	  for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
 	    DegreeOfFreedom d = it->first;
-	    int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
-	    int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
+	    int matIndexLocal = localDofMap.getLocalMatIndex(component, d);
+	    int matIndexDual = dualDofMap.getLocalMatIndex(component, d);
 	    lumpedData->localToDualMap[matIndexLocal] = matIndexDual;
 	  }
 	}
@@ -1588,7 +1585,7 @@ namespace AMDiS {
 	std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
 	for (std::set<DegreeOfFreedom>::iterator it = dirichletRows.begin();
 	     it != dirichletRows.end(); ++it) {
-	  if (localDofMap[feSpace].isSet(*it)) {
+	  if (localDofMap[i].isSet(*it)) {
 	    MSG("Dirichlet dof %d in component %d with interior mat index %d\n",
 		*it, i, localDofMap.getMatIndex(i, *it));
 	  }
@@ -1640,7 +1637,7 @@ namespace AMDiS {
     // === Create scatter to get solutions of all primal nodes that are ===
     // === contained in rank's domain.                                  ===
     
-    vector<const FiniteElemSpace*> feSpaces = vec.getComponentFeSpaces();
+    int nComponents = vec.getSize();
 
     vector<PetscInt> globalIsIndex, localIsIndex;
     globalIsIndex.reserve(primalDofMap.getLocalDofs());
@@ -1648,10 +1645,10 @@ namespace AMDiS {
 
     {
       int cnt = 0;
-      for (unsigned int i = 0; i < feSpaces.size(); i++) {
-	DofMap& dofMap = primalDofMap[feSpaces[i]].getMap();
+      for (unsigned int component = 0; component < nComponents; component++) {
+	DofMap& dofMap = primalDofMap[component].getMap();
 	for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) {
-	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
+	  globalIsIndex.push_back(primalDofMap.getMatIndex(component, it->first));
 	  localIsIndex.push_back(cnt++);
 	}
       }
@@ -1692,26 +1689,27 @@ namespace AMDiS {
 
     // === And copy from PETSc local vectors to the DOF vectors. ===
     
+    vector<const FiniteElemSpace*> feSpaces = vec.getComponentFeSpaces();
     int cnt = 0;
-    for (int i = 0; i < vec.getSize(); i++) {
-      DOFVector<double>& dofVec = *(vec.getDOFVector(i));
+    for (int component = 0; component < nComponents; component++) {
+      DOFVector<double>& dofVec = *(vec.getDOFVector(component));
       
-      for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap().begin();
-	   it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
+      for (DofMap::iterator it = localDofMap[component].getMap().begin();
+	   it != localDofMap[component].getMap().end(); ++it) {
 	if (meshLevel == 0) {
-	  int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
+	  int petscIndex = localDofMap.getLocalMatIndex(component, it->first);
 	  dofVec[it->first] = localSolB[petscIndex];
 	} else {
-	  if (meshDistributor->getDofMapSd()[feSpaces[i]].isRankDof(it->first)) {
-	    int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
+	  if (meshDistributor->getDofMapSd()[feSpaces[component]].isRankDof(it->first)) {
+	    int petscIndex = localDofMap.getLocalMatIndex(component, it->first);
 	    TEST_EXIT(petscIndex < bsize)("Should not happen!\n");
 	    dofVec[it->first] = localSolB[petscIndex];	 
 	  }
 	}
       }
 
-      for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
-	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
+      for (DofMap::iterator it = primalDofMap[component].getMap().begin();
+	   it != primalDofMap[component].getMap().end(); ++it)
 	dofVec[it->first] = localSolPrimal[cnt++];      
     }