Fixed problem for Dirichlet boundary conditions on non diagonal blocks in parallel.

AMDiS/src/DirichletBC.h

@@ -108,10 +108,8 @@ namespace AMDiS {
     /// DOFVector containing the boundary values
     DOFVectorBase<double> *dofVec;
 
-    /** \brief
-     * Defines, if the boundary condition must be applied to the matrix. See
-     * comment of \ref BoundaryCondition::applyBoundaryCondition.
-     */
+    /// Defines, if the boundary condition must be applied to the matrix. See
+    /// comment of \ref BoundaryCondition::applyBoundaryCondition.
     bool applyBC;
   };
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -476,33 +476,65 @@ namespace AMDiS {
       return;
 
     int nComponents = seqMat->getSize();
-    vector<PetscInt> zeroRows;
+    vector<int> dirichletRows;
+    vector<int> dirichletCols;
+    vector<double> dirichletValues;
 
     for (int i = 0; i < nComponents; i++) {
-      DOFMatrix *dofMat = (*seqMat)[i][i];
-      if (!dofMat)
+      bool dirichletRow = false;
+      int dirichletMainCol = -1;
+
+      for (int j = 0; j < nComponents; j++) {
+	DOFMatrix *dofMat = (*seqMat)[i][j];
+	if (!dofMat)
+	  continue;
+
+	BoundaryIndexMap& bounds = 
+	  const_cast<BoundaryIndexMap&>(dofMat->getBoundaryManager()->getBoundaryConditionMap());
+
+	for (BoundaryIndexMap::iterator bIt = bounds.begin(); bIt != bounds.end(); ++bIt) {
+	  if (bIt->second && bIt->second->isDirichlet()) {
+	    dirichletRow = true;
+	    if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
+	      TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n");
+	      dirichletMainCol = j;
+	    } 
+	  }
+	}
+      }
+
+      if (!dirichletRow)
 	continue;
 
+      DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol];
       const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
+      TEST_EXIT(dofMat->getRowFeSpace() == dofMat->getColFeSpace())
+	("I have to think about this scenario! Really possible?\n");
 
       std::set<DegreeOfFreedom> &dRows = dofMat->getDirichletRows();
       for (std::set<DegreeOfFreedom>::iterator dofIt = dRows.begin();
 	   dofIt != dRows.end(); ++dofIt) {
 	MultiIndex multiIndex;
 	if (dofMap[feSpace].find(*dofIt, multiIndex) &&
-	    dofMap[feSpace].isRankDof(*dofIt))
-	  zeroRows.push_back(dofMap.getMatIndex(i, multiIndex.global));
+	    dofMap[feSpace].isRankDof(*dofIt)) {
+	  int rowIndex = dofMap.getMatIndex(i, multiIndex.global);
+	  int colIndex = dofMap.getMatIndex(dirichletMainCol, multiIndex.global);
+	  dirichletRows.push_back(rowIndex);
+	  dirichletCols.push_back(colIndex);
+	  dirichletValues.push_back(1.0);
+	}
       }
     }
 
-    MatZeroRows(mpiMat, zeroRows.size(), &(zeroRows[0]), 1.0, 
+    MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 1.0,
 		PETSC_NULL, PETSC_NULL);
 
-    PetscViewer view;
-    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "interior.mat",
-			  FILE_MODE_WRITE, &view);
-    MatView(mpiMat, view);
-    PetscViewerDestroy(&view);
+    for (int i = 0; i < static_cast<int>(dirichletRows.size()); i++) {
+      MatSetValue(mpiMat, dirichletRows[i], dirichletCols[i], dirichletValues[i], INSERT_VALUES);
+    }
+    
+    MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
   }