From 82a265365ff96bb2616ea1a4cd1d2b978051fd95 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 28 Aug 2012 11:05:53 +0000
Subject: [PATCH] Work on lumped preconditioner for FETI-DP Stokes.

---
 AMDiS/src/parallel/PetscSolverFeti.cc         | 402 +++++++++---------
 AMDiS/src/parallel/PetscSolverFeti.h          |   3 +
 .../src/parallel/PetscSolverFetiOperators.cc  |  21 +-
 AMDiS/src/parallel/PetscSolverFetiStructs.h   |   4 +
 4 files changed, 235 insertions(+), 195 deletions(-)

diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index f1566373..cab18244 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -108,7 +108,7 @@ namespace AMDiS {
     if (fullInterface != NULL)
       interfaceDofMap.init(levelData, feSpaces, uniqueFe);
 
-    if (fetiPreconditioner != FETI_NONE) {
+    if (fetiPreconditioner == FETI_DIRICHLET) {
       TEST_EXIT(meshLevel == 0)
 	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");
 
@@ -133,7 +133,7 @@ namespace AMDiS {
     dualDofMap.clear();
     lagrangeMap.clear();
     localDofMap.clear();
-    if (fetiPreconditioner != FETI_NONE)
+    if (fetiPreconditioner == FETI_DIRICHLET)
       interiorDofMap.clear();
 
     primalDofMap.setDofComm(meshDistributor->getDofComm());
@@ -143,7 +143,7 @@ namespace AMDiS {
     dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
     lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
     localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
-    if (fetiPreconditioner != FETI_NONE)
+    if (fetiPreconditioner == FETI_DIRICHLET)
       interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
 
     if (meshLevel == 0)
@@ -172,7 +172,7 @@ namespace AMDiS {
     primalDofMap.update();
     dualDofMap.update();
     localDofMap.update();
-    if (fetiPreconditioner != FETI_NONE)
+    if (fetiPreconditioner == FETI_DIRICHLET)
       interiorDofMap.update();
 
     if (fullInterface != NULL)
@@ -504,7 +504,7 @@ namespace AMDiS {
 	if (meshLevel == 0) {
 	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
 
-	  if (fetiPreconditioner != FETI_NONE)
+	  if (fetiPreconditioner == FETI_DIRICHLET)
 	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
 	  
 	  nLocalInterior++;	  
@@ -820,8 +820,8 @@ namespace AMDiS {
       
       MatNullSpace matNullspace;
       MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace);
-      // MatSetNullSpace(mat_feti, matNullspace);
-      // KSPSetNullSpace(ksp_feti, matNullspace);
+      MatSetNullSpace(mat_feti, matNullspace);
+      KSPSetNullSpace(ksp_feti, matNullspace);
     }
 
 
@@ -834,10 +834,13 @@ namespace AMDiS {
 
 
     switch (fetiPreconditioner) {
-    case FETI_DIRICHLET:      
-      TEST_EXIT_DBG(meshLevel == 0)
+    case FETI_DIRICHLET:
+      TEST_EXIT(meshLevel == 0)
 	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");
 
+      TEST_EXIT(!enableStokesMode)
+	("Stokes mode does not yet support the Dirichlet precondition!\n");
+
       KSPCreate(PETSC_COMM_SELF, &ksp_interior);
       KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
 		      SAME_NONZERO_PATTERN);
@@ -898,10 +901,14 @@ namespace AMDiS {
       break;
 
     case FETI_LUMPED:
+      fetiLumpedPreconData.enableStokesMode = enableStokesMode;
       fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
       fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;
 
       for (unsigned int i = 0; i < feSpaces.size(); i++) {
+	if (enableStokesMode && feSpaces[i] == fullInterface)
+	  continue;
+
 	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
 	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
 	  DegreeOfFreedom d = it->first;
@@ -1347,182 +1354,9 @@ namespace AMDiS {
       MPI::COMM_WORLD.Barrier();
       MSG("FETI-DP timing 02: %.5f seconds (creation of interior matrices)\n",
 	  MPI::Wtime() - wtime);
-    }
-   
-    
-    // === Create matrices for FETI-DP preconditioner. ===
-    
-    if (fetiPreconditioner != FETI_NONE) {
-      wtime = MPI::Wtime();
-
-      int nRowsDual = dualDofMap.getRankDofs();
-
-      MatCreateSeqAIJ(PETSC_COMM_SELF,
-		      nRowsDual, nRowsDual, 100, PETSC_NULL,
-		      &mat_duals_duals);
-      MatSetOption(mat_duals_duals, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
-
-      if (fetiPreconditioner == FETI_DIRICHLET) {
-	int nRowsInterior = interiorDofMap.getRankDofs();
-
-	MatCreateSeqAIJ(PETSC_COMM_SELF,
-			nRowsInterior, nRowsInterior, 100, PETSC_NULL,
-			&mat_interior_interior);
-	MatSetOption(mat_interior_interior, 
-		     MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
-
-	MatCreateSeqAIJ(PETSC_COMM_SELF,
-			nRowsInterior, nRowsDual, 100, PETSC_NULL,
-			&mat_interior_duals);
-	MatSetOption(mat_interior_duals, 
-		     MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
-	
-	MatCreateSeqAIJ(PETSC_COMM_SELF,
-			nRowsDual, nRowsInterior, 100, PETSC_NULL,
-			&mat_duals_interior);
-	MatSetOption(mat_duals_interior, 
-		     MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
-      }
-    
-      // === Prepare traverse of sequentially created matrices. ===
-      
-      using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
-      namespace traits = mtl::traits;
-      typedef DOFMatrix::base_matrix_type Matrix;
-      
-      typedef traits::range_generator<row, Matrix>::type cursor_type;
-      typedef traits::range_generator<nz, cursor_type>::type icursor_type;
-      
-      vector<int> colsLocal, colsLocalOther;
-      vector<double> valuesLocal, valuesLocalOther;
-      colsLocal.reserve(300);
-      colsLocalOther.reserve(300);
-      valuesLocal.reserve(300);
-      valuesLocalOther.reserve(300);
-      
-      
-      // === Traverse all sequentially created matrices and add the values to ===
-      // === the global PETSc matrices.                                       ===
-
-      int nComponents = mat->getSize();
-      for (int i = 0; i < nComponents; i++) {
-	for (int j = 0; j < nComponents; j++) {
-	  if (!(*mat)[i][j])
-	    continue;
-
-	  traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
-	  traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
-	
-	  // Traverse all rows.
-	  for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
-		 cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
-
-	    bool rowPrimal = isPrimal(feSpaces[i], *cursor);
-  
-	    colsLocal.clear();
-	    colsLocalOther.clear();
-	    valuesLocal.clear();
-	    valuesLocalOther.clear();
-
-	  
-	    // Traverse all columns.
-	    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
-		 icursor != icend; ++icursor) {
-
-	      bool colPrimal = isPrimal(feSpaces[j], col(*icursor));
-
-	      if (!rowPrimal && !colPrimal) {
-		if (!isDual(feSpaces[i], *cursor)) {
-		  if (!isDual(feSpaces[j], col(*icursor))) {
-		    int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
-		    colsLocal.push_back(colIndex);
-		    valuesLocal.push_back(value(*icursor));
-		  } else {
-		    int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
-		    colsLocalOther.push_back(colIndex);
-		    valuesLocalOther.push_back(value(*icursor));
-		  }
-		} else {
-		  if (!isDual(feSpaces[j], col(*icursor))) {
-		    int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
-		    colsLocalOther.push_back(colIndex);
-		    valuesLocalOther.push_back(value(*icursor));
-		  } else {
-		    int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
-		    colsLocal.push_back(colIndex);
-		    valuesLocal.push_back(value(*icursor));
-		  }
-		}		
-	      }
-	    }  // for each nnz in row
-
-
-	    // === Set matrix values for preconditioner ===
-
-	    if (!rowPrimal) {
-	      switch (fetiPreconditioner) {
-	      case FETI_DIRICHLET:
-		if (!isDual(feSpaces[i], *cursor)) {
-		  int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor);
-		  MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
-			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
-		  if (colsLocalOther.size()) 
-		    MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
-				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
-		} else {
-		  int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
-		  MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
-			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);  
-		  if (colsLocalOther.size()) 
-		    MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
-				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
-		}
-		break;
-
-	      case FETI_LUMPED:
-		if (isDual(feSpaces[i], *cursor)) {
-		  int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);		
-		  MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
-			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
-		}			
-		break;
-	      default:
-		break;
-	      }	  
-	    }
-	  } 
-	}
-      }
-
-      // === Start global assembly procedure for preconditioner matrices. ===
 
-      if (fetiPreconditioner != FETI_NONE) {
-	MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
-      }
-
-      if (fetiPreconditioner == FETI_DIRICHLET) {
-	MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
-      
-	MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
-      
-	MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
-      }
-
-      if (printTimings) {
-	MPI::COMM_WORLD.Barrier();
-	MSG("FETI-DP timing 03: %.5f seconds (creation of preconditioner matrices)\n",
-	    MPI::Wtime() - wtime);
-      }
-    }
-
-
-    // For the case, that we want to print the timings, we force the LU 
-    // factorization of the local problems to be done here explicitly.
-    if (printTimings) {
+      // For the case, that we want to print the timings, we force the LU 
+      // factorization of the local problems to be done here explicitly.
       wtime = MPI::Wtime();
       KSPSetUp(subdomain->getSolver());
       KSPSetUpOnBlocks(subdomain->getSolver());
@@ -1530,25 +1364,21 @@ namespace AMDiS {
       MSG("FETI-DP timing 04: %.5f seconds (factorization of subdomain matrices)\n",
 	  MPI::Wtime() - wtime);
     }
+   
+    
+    // === Create matrices for FETI-DP preconditioner. ===    
+    createPreconditionerMatrix(mat);
     
-
     // === Create and fill PETSc matrix for Lagrange constraints. ===
-
     createMatLagrange(feSpaces);
 
-
     // === If required, run debug tests. ===
-
     dbgMatrix();
 
-
-    // === Create PETSc solver for the Schur complement on primal variables. ===
-    
+    // === Create PETSc solver for the Schur complement on primal variables. ===   
     createSchurPrimalKsp(feSpaces);
 
-
     // === Create PETSc solver for the FETI-DP operator. ===
-
     createFetiKsp(feSpaces);
   }
 
@@ -1560,7 +1390,193 @@ namespace AMDiS {
     subdomain->fillPetscRhs(vec);
   }
 
+
+  void PetscSolverFeti::createPreconditionerMatrix(Matrix<DOFMatrix*> *mat)
+  {
+    FUNCNAME("PetscSolverFeti::createPreconditionerMatrix()");
+
+    if (fetiPreconditioner != FETI_NONE)
+      return;
+
+    double wtime = MPI::Wtime();
+    vector<const FiniteElemSpace*> feSpaces = AMDiS::getComponentFeSpaces(*mat);
+    int nRowsDual = dualDofMap.getRankDofs();
+    
+    MatCreateSeqAIJ(PETSC_COMM_SELF,
+		    nRowsDual, nRowsDual, 100, PETSC_NULL,
+		    &mat_duals_duals);
+    MatSetOption(mat_duals_duals, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+    
+    if (fetiPreconditioner == FETI_DIRICHLET) {
+      int nRowsInterior = interiorDofMap.getRankDofs();
+      
+      MatCreateSeqAIJ(PETSC_COMM_SELF,
+		      nRowsInterior, nRowsInterior, 100, PETSC_NULL,
+		      &mat_interior_interior);
+      MatSetOption(mat_interior_interior, 
+		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+      
+      MatCreateSeqAIJ(PETSC_COMM_SELF,
+		      nRowsInterior, nRowsDual, 100, PETSC_NULL,
+		      &mat_interior_duals);
+      MatSetOption(mat_interior_duals, 
+		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+      
+      MatCreateSeqAIJ(PETSC_COMM_SELF,
+		      nRowsDual, nRowsInterior, 100, PETSC_NULL,
+		      &mat_duals_interior);
+      MatSetOption(mat_duals_interior, 
+		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+    }
+    
+    // === Prepare traverse of sequentially created matrices. ===
+    
+    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
+    namespace traits = mtl::traits;
+    typedef DOFMatrix::base_matrix_type Matrix;
+    
+    typedef traits::range_generator<row, Matrix>::type cursor_type;
+    typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+    
+    vector<int> colsLocal, colsLocalOther;
+    vector<double> valuesLocal, valuesLocalOther;
+    colsLocal.reserve(300);
+    colsLocalOther.reserve(300);
+    valuesLocal.reserve(300);
+    valuesLocalOther.reserve(300);
+    
+    
+    // === Traverse all sequentially created matrices and add the values to ===
+    // === the global PETSc matrices.                                       ===
+    
+    int nComponents = mat->getSize();
+    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
+      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
+	DOFMatrix* dofMat = (*mat)[rowComponent][colComponent];
+
+	if (!dofMat)
+	  continue;
+	
+	const FiniteElemSpace *rowFeSpace = feSpaces[rowComponent];
+	const FiniteElemSpace *colFeSpace = feSpaces[colComponent];
+
+	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
+	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
+	
+	// Traverse all rows.
+	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
+	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
+	  
+	  if (isPrimal(rowFeSpace, *cursor))
+	    continue;
+	  
+	  colsLocal.clear();
+	  colsLocalOther.clear();
+	  valuesLocal.clear();
+	  valuesLocalOther.clear();	  
+	  
+	  // Traverse all columns.
+	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
+	       icursor != icend; ++icursor) {
+	    
+	    if (isPrimal(colFeSpace, col(*icursor)))
+	      continue;
+	    
+	    if (!isDual(rowFeSpace, *cursor)) {
+	      if (!isDual(colFeSpace, col(*icursor))) {
+		int colIndex = 
+		  interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
+		colsLocal.push_back(colIndex);
+		valuesLocal.push_back(value(*icursor));
+	      } else {
+		int colIndex = 
+		  dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
+		colsLocalOther.push_back(colIndex);
+		valuesLocalOther.push_back(value(*icursor));
+	      }
+	    } else {
+	      if (!isDual(colFeSpace, col(*icursor))) {
+		int colIndex = 
+		  interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
+		colsLocalOther.push_back(colIndex);
+		valuesLocalOther.push_back(value(*icursor));
+	      } else {
+		int colIndex = 
+		  dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
+		colsLocal.push_back(colIndex);
+		valuesLocal.push_back(value(*icursor));
+	      }
+	    }			    
+	  }  // for each nnz in row
+	  
+	  
+	  // === Set matrix values for preconditioner ===
+	  
+	  switch (fetiPreconditioner) {
+	  case FETI_DIRICHLET:
+	    if (!isDual(rowFeSpace, *cursor)) {
+	      int rowIndex = 
+		interiorDofMap.getLocalMatIndex(rowComponent, *cursor);
+	      MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
+			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+	      
+	      if (colsLocalOther.size()) 
+		MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
+			     &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+	    } else {
+	      int rowIndex = 
+		dualDofMap.getLocalMatIndex(rowComponent, *cursor);
+	      MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
+			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);  
+	      if (colsLocalOther.size()) 
+		MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
+			     &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+	    }
+	    break;
+	    
+	  case FETI_LUMPED:
+	    if (isDual(rowFeSpace, *cursor)) {
+	      int rowIndex = dualDofMap.getLocalMatIndex(rowComponent, *cursor);		
+	      MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
+			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+	    }	
+	    break;
+	  default:
+	    break;
+	  }	  
+	} 
+      }
+    }
+    
+    // === Start global assembly procedure for preconditioner matrices. ===
+    
+    if (fetiPreconditioner == FETI_LUMPED) {
+      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
+    }
+    
+    if (fetiPreconditioner == FETI_DIRICHLET) {
+      MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
+      
+      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
+      
+      MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
+      
+      MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
+    }
+    
+    if (printTimings) {
+      MPI::COMM_WORLD.Barrier();
+      MSG("FETI-DP timing 03: %.5f seconds (creation of preconditioner matrices)\n",
+	  MPI::Wtime() - wtime);
+    }   
+  }
   
+
   void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
   {
     FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index eb86b979..16f3c0d5 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -125,6 +125,9 @@ namespace AMDiS {
     /// to \ref mat_lagrange.
     void createMatLagrange(vector<const FiniteElemSpace*> &feSpaces);
 
+    ///
+    void createPreconditionerMatrix(Matrix<DOFMatrix*> *mat);
+
     /// Creates PETSc KSP solver object for solving the Schur complement
     /// system on the primal variables, \ref ksp_schur_primal
     void createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces);
diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
index 1cf7a349..cd17d10e 100644
--- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc
+++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
@@ -238,8 +238,25 @@ namespace AMDiS {
     PCShellGetContext(pc, &ctx);
     FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);
 
+    Vec xvec, yvec;
+    if (data->enableStokesMode) {
+      Vec x_interface, x_lagrange, y_interface, y_lagrange;    
+      VecNestGetSubVec(x, 0, &x_interface);
+      VecNestGetSubVec(x, 1, &x_lagrange);
+      VecNestGetSubVec(y, 0, &y_interface);
+      VecNestGetSubVec(y, 1, &y_lagrange);
+      
+      xvec = x_lagrange;
+      yvec = y_lagrange;
+
+      VecCopy(x_interface, y_interface);
+    } else {
+      xvec = x;
+      yvec = y;
+    }
+
     // Multiply with scaled Lagrange constraint matrix.
-    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
+    MatMultTranspose(*(data->mat_lagrange_scaled), xvec, data->tmp_vec_b);
 
 
     // === Restriction of the B nodes to the boundary nodes. ===
@@ -280,7 +297,7 @@ namespace AMDiS {
 
 
     // Multiply with scaled Lagrange constraint matrix.
-    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
+    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, yvec);
 
     return 0;
   }
diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h
index 3da72354..2de8c305 100644
--- a/AMDiS/src/parallel/PetscSolverFetiStructs.h
+++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h
@@ -99,6 +99,10 @@ namespace AMDiS {
 
 
   struct FetiLumpedPreconData {
+    /// Defines whether the precondition is applied for a Stokes like FETI-DP
+    /// problem with interface variables which handled in a different way.
+    bool enableStokesMode;
+
     /// Matrix of scaled Lagrange variables.
     Mat *mat_lagrange_scaled;
 
-- 
GitLab