FETI-DP Dirichlet precondition works now again.

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -88,8 +88,9 @@ namespace AMDiS {
     VecGetArray(data->tmp_vec_b, &local_b);
     VecGetArray(data->tmp_vec_duals0, &local_duals);
 
-    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
-      local_duals[j] = local_b[i];
+    for (map<int, int>::iterator it = data->localToDualMap.begin();
+	 it != data->localToDualMap.end(); ++it)
+      local_duals[it->second] = local_b[it->first];
 
     VecRestoreArray(data->tmp_vec_b, &local_b);
     VecRestoreArray(data->tmp_vec_duals0, &local_duals);
@@ -111,8 +112,9 @@ namespace AMDiS {
     VecGetArray(data->tmp_vec_b, &local_b);
     VecGetArray(data->tmp_vec_duals0, &local_duals);
 
-    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
-      local_b[i] = local_duals[j];
+    for (map<int, int>::iterator it = data->localToDualMap.begin();
+	 it != data->localToDualMap.end(); ++it)
+      local_b[it->first] = local_duals[it->second];
 
     VecRestoreArray(data->tmp_vec_b, &local_b);
     VecRestoreArray(data->tmp_vec_duals0, &local_duals);
@@ -674,7 +676,19 @@ namespace AMDiS {
 		 &(fetiDirichletPreconData.tmp_vec_duals1));
       MatGetVecs(mat_interior_interior, PETSC_NULL, 
 		 &(fetiDirichletPreconData.tmp_vec_interior));
-      
+
+      for (unsigned int i = 0; i < feSpaces.size(); i++) {
+	map<DegreeOfFreedom, MultiIndex> &dualMap = 
+	  dualDofMap[feSpaces[i]].getMap();
+	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
+	     it != dualMap.end(); ++it) {
+	  DegreeOfFreedom d = it->first;
+	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
+	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
+	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
+	}
+      }
+
       KSPGetPC(ksp_feti, &precon_feti);
       PCSetType(precon_feti, PCSHELL);
       PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));

AMDiS/src/parallel/PetscSolverFetiStructs.h

@@ -23,8 +23,12 @@
 #ifndef AMDIS_PETSC_SOLVER_FETI_STRUCTS_H
 #define AMDIS_PETSC_SOLVER_FETI_STRUCTS_H
 
+#include <map>
+
 namespace AMDiS {
 
+  using namespace std;
+
   class PetscSolverFeti;
 
   /** \brief
@@ -92,7 +96,9 @@ namespace AMDiS {
     Vec tmp_vec_duals0, tmp_vec_duals1;
     
     /// Temporal vector on the interior variables.
-    Vec tmp_vec_interior;    
+    Vec tmp_vec_interior;
+
+    map<int, int> localToDualMap;
   };