Fixed problem with dirichlet boundary dofs in FETI-DP code.

AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc

@@ -154,14 +154,13 @@ namespace AMDiS {
       MatCreateSeqAIJ(mpiCommLocal, nRankInteriorRows, nRankInteriorRows,
 		      0, nnz[0][0].dnnz,
 		      &mat[0][0]);
-      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     } else {
       MatCreateAIJ(mpiCommGlobal, nRankInteriorRows, nRankInteriorRows, 
 		   nOverallInteriorRows, nOverallInteriorRows,
 		   0, nnz[0][0].dnnz, 0, nnz[0][0].onnz,		   
 		   &mat[0][0]);    
-      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
     }
+    MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
     VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecSol[0]);
     VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecRhs[0]);
@@ -178,6 +177,7 @@ namespace AMDiS {
 		   nOverallCoarseRows, nOverallCoarseRows,
 		   0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
 		   &mat[i + 1][i + 1]);
+      MatSetOption(mat[i + 1][i + 1], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
       cMap->createVec(vecSol[i + 1]);
       cMap->createVec(vecRhs[i + 1]);
     }

AMDiS/src/parallel/ParallelCoarseSpaceSolver.h

@@ -234,6 +234,12 @@ namespace AMDiS {
     }
 
 
+    /// Returns whether the solver has a coarse grid.
+    inline bool hasCoarseSpace() 
+    {
+      return (!coarseSpaceMap.empty());
+    }
+
   protected:
     /// Prepare internal data structures. First, it create \ref uniqueCoarseMap
     /// and \ref componentIthCoarseMap . Both are used to create the correct 

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -108,7 +108,8 @@ namespace AMDiS {
     if (subdomain == NULL) {
       subdomain = new PetscSolverGlobalMatrix("");
       subdomain->setSymmetric(isSymmetric);
-      subdomain->setHandleDirichletRows(false);
+      //      subdomain->setHandleDirichletRows(false);
+      subdomain->setHandleDirichletRows(true);
 
       if (meshLevel == 0) {
 	subdomain->setMeshDistributor(meshDistributor, 
@@ -142,6 +143,8 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverFeti::createDirichletData()");
 
+    return;
+
     int nComponents = mat.getSize();
     for (int component = 0; component < nComponents; component++) {
       DOFMatrix* dofMat = mat[component][component];

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -17,6 +17,17 @@
 
 namespace AMDiS {
 
+
+  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
+    : PetscSolver(name),
+      zeroStartVector(false),
+      printMatInfo(false)
+  {
+    Parameters::get("parallel->use zero start vector", zeroStartVector);
+    Parameters::get("parallel->print matrix info", printMatInfo);
+  }
+
+
   void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
   {
     FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");
@@ -54,7 +65,7 @@ namespace AMDiS {
 
     matAssembly();
 
-    removeDirichletRows(seqMat, dofMap, getMatInterior());
+    removeDirichletRows(seqMat);
     
     if (printMatInfo) {
       MatInfo matInfo;
@@ -225,6 +236,8 @@ namespace AMDiS {
 
     matAssembly();
 
+    removeDirichletRows(seqMat);
+
     // === Create solver for the non primal (thus local) variables. ===
 
     KSPCreate(mpiCommLocal, &kspInterior);
@@ -266,7 +279,7 @@ namespace AMDiS {
 
     vecRhsAssembly();
 
-    removeDirichletRows(vec, dofMap, getVecRhsInterior());
+    removeDirichletRows(vec);
 
     // === For debugging allow to write the rhs vector to a file. ===
 
@@ -456,26 +469,25 @@ namespace AMDiS {
   }
 
 
-  void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat, 
-						    ParallelDofMapping &dofMap,
-						    Mat mpiMat)
+  void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat)					       
   {
     FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");
 
     if (!handleDirichletRows)
       return;
 
+    vector<int> dRowsInterior, dRowsCoarse;
+    vector<int> dColsInterior, dColsCoarse;
+    vector<double> dValuesInterior, dValuesCoarse;
+
     int nComponents = seqMat->getSize();
-    vector<int> dirichletRows;
-    vector<int> dirichletCols;
-    vector<double> dirichletValues;
 
-    for (int i = 0; i < nComponents; i++) {
+    for (int rowComp = 0; rowComp < nComponents; rowComp++) {
       bool dirichletRow = false;
       int dirichletMainCol = -1;
 
-      for (int j = 0; j < nComponents; j++) {
-	DOFMatrix *dofMat = (*seqMat)[i][j];
+      for (int colComp = 0; colComp < nComponents; colComp++) {
+	DOFMatrix *dofMat = (*seqMat)[rowComp][colComp];
 	if (!dofMat)
 	  continue;
 
@@ -486,7 +498,7 @@ namespace AMDiS {
 	  if (bIt->second && bIt->second->isDirichlet()) {
 	    dirichletRow = true;
 	    if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
-	      dirichletMainCol = j;
+	      dirichletMainCol = colComp;
 	      break;
 	    } 
 	  }
@@ -496,65 +508,152 @@ namespace AMDiS {
       if (!dirichletRow)
 	continue;
 
-      DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol];
-      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
+      DOFMatrix *dofMat = (*seqMat)[rowComp][dirichletMainCol];
       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) {
+	if (hasCoarseSpace()) {
+	  bool isRowCoarse = isCoarseSpace(rowComp, *dofIt);
+	  bool isColCoarse = isCoarseSpace(dirichletMainCol, *dofIt);
+	  TEST_EXIT(isRowCoarse == isColCoarse)
+	    ("Really possible? So reimplement AMDiS from the scratch!\n");
+
+	  if (isRowCoarse) {
+	    ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComp];
+	    ParallelDofMapping *colCoarseSpace = coarseSpaceMap[dirichletMainCol];
+
 	    MultiIndex multiIndex;
-	if (dofMap[feSpace].find(*dofIt, multiIndex) &&
-	    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);
+	    if ((*rowCoarseSpace)[rowComp].find(*dofIt, multiIndex) &&
+		(*rowCoarseSpace)[rowComp].isRankDof(*dofIt)) {
+	      int rowIndex = rowCoarseSpace->getLocalMatIndex(rowComp, *dofIt);
+	      int colIndex = colCoarseSpace->getLocalMatIndex(dirichletMainCol, *dofIt);
+	      dRowsCoarse.push_back(rowIndex);
+	      dColsCoarse.push_back(colIndex);
+	      dValuesCoarse.push_back(1.0);
+	    }
+	  } else {
+
+	    MultiIndex multiIndex;
+	    if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
+		(*interiorMap)[rowComp].isRankDof(*dofIt)) {
+	      int rowIndex = interiorMap->getLocalMatIndex(rowComp, *dofIt);
+	      int colIndex = interiorMap->getLocalMatIndex(dirichletMainCol, *dofIt);
+	      dRowsInterior.push_back(rowIndex);
+	      dColsInterior.push_back(colIndex);
+	      dValuesInterior.push_back(1.0);
+	    }
+	  }
+	} else {
+	  MultiIndex multiIndex;
+	  if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
+	      (*interiorMap)[rowComp].isRankDof(*dofIt)) {
+	    int rowIndex = interiorMap->getMatIndex(rowComp, multiIndex.global);
+	    int colIndex = interiorMap->getMatIndex(dirichletMainCol, multiIndex.global);
+	    dRowsInterior.push_back(rowIndex);
+	    dColsInterior.push_back(colIndex);
+	    dValuesInterior.push_back(1.0);
+	  }
 	}
       }
     }
 
-    MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 0.0,
+    {
+      Mat mpiMat = getMatInterior();
+      MatZeroRows(mpiMat, dRowsInterior.size(), &(dRowsInterior[0]), 0.0,
 		  PETSC_NULL, PETSC_NULL);
       
-    for (int i = 0; i < static_cast<int>(dirichletRows.size()); i++)
-      MatSetValue(mpiMat, dirichletRows[i], dirichletCols[i], dirichletValues[i], INSERT_VALUES);    
+      for (int i = 0; i < static_cast<int>(dRowsInterior.size()); i++)
+	MatSetValue(mpiMat, dRowsInterior[i], dColsInterior[i], 
+		    dValuesInterior[i], INSERT_VALUES);    
       
       MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
       MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
     }
 
+    if (hasCoarseSpace()) {
+      Mat mpiMat = getMatCoarse();
+      MatZeroRows(mpiMat, dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0,
+		  PETSC_NULL, PETSC_NULL);
+      
+      for (int i = 0; i < static_cast<int>(dRowsCoarse.size()); i++)
+	MatSetValue(mpiMat, dRowsCoarse[i], dColsCoarse[i], 
+		    dValuesCoarse[i], INSERT_VALUES);    
       
-  void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec, 
-						    ParallelDofMapping &dofMap,
-						    Vec mpiVec)
+      MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
+    }
+  }
+
+  
+  void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec)
   {
     FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");
 
     if (!handleDirichletRows)
       return;
 
-    int nComponents = seqVec->getSize();
-
-    for (int i = 0; i < nComponents; i++) {
-      DOFVector<double> *dofVec = seqVec->getDOFVector(i);
-      const FiniteElemSpace *feSpace = dofVec->getFeSpace();
+    int cInterior = 0;
+    int cCoarse = 0;
 
+    int nComponents = seqVec->getSize();
+    for (int component = 0; component < nComponents; component++) {
+      DOFVector<double> *dofVec = seqVec->getDOFVector(component);
       map<DegreeOfFreedom, double>& dValues = dofVec->getDirichletValues();
+
       for (map<DegreeOfFreedom, double>::iterator dofIt = dValues.begin();
 	   dofIt != dValues.end(); ++dofIt) {
+	if (hasCoarseSpace()) {
+	  if (isCoarseSpace(component, dofIt->first)) {
+	    ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[component];
+	    MultiIndex multiIndex;
+	    if ((*rowCoarseSpace)[component].find(dofIt->first, multiIndex) &&
+		(*rowCoarseSpace)[component].isRankDof(dofIt->first)) {
+	      if (dofIt->second > 0.5)
+		cCoarse++;
+	      VecSetValue(getVecRhsCoarse(), 
+			  rowCoarseSpace->getMatIndex(component, multiIndex.global),
+			  dofIt->second, INSERT_VALUES);
+	    }
+	  } else {
+	    MultiIndex multiIndex;
+	    if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
+		(*interiorMap)[component].isRankDof(dofIt->first)) {
+	      if (dofIt->second > 0.5)
+		cInterior++;
+	      VecSetValue(getVecRhsInterior(), 
+			  interiorMap->getLocalMatIndex(component, dofIt->first),
+			  dofIt->second, INSERT_VALUES);
+	    }
+	  }
+	} else {
 	  MultiIndex multiIndex;
-	if (dofMap[feSpace].find(dofIt->first, multiIndex) &&
-	    dofMap[feSpace].isRankDof(dofIt->first))
-	  VecSetValue(mpiVec, dofMap.getMatIndex(i, multiIndex.global),
+	  if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
+	      (*interiorMap)[component].isRankDof(dofIt->first)) {
+	    if (dofIt->second > 0.5)
+	      cInterior++;
+	    VecSetValue(getVecRhsInterior(), 
+			interiorMap->getMatIndex(component, multiIndex.global),
 			dofIt->second, INSERT_VALUES);
 	  }
 	}
+      }
+    }
+
+    VecAssemblyBegin(getVecRhsInterior());
+    VecAssemblyEnd(getVecRhsInterior());
+
+    if (hasCoarseSpace()) {
+      VecAssemblyBegin(getVecRhsCoarse());
+      VecAssemblyEnd(getVecRhsCoarse());
+    }
 
-    VecAssemblyBegin(mpiVec);
-    VecAssemblyEnd(mpiVec);
+    mpi::globalAdd(cInterior);
+    mpi::globalAdd(cCoarse);
+    MSG("WORKED ON DIRICHLET DOFS: %d %d\n", cInterior, cCoarse);
   }
 
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -39,14 +39,7 @@ namespace AMDiS {
   class PetscSolverGlobalMatrix : public PetscSolver
   {
   public:
-    PetscSolverGlobalMatrix(string name)
-      : PetscSolver(name),
-	zeroStartVector(false),
-	printMatInfo(false)
-    {
-      Parameters::get("parallel->use zero start vector", zeroStartVector);
-      Parameters::get("parallel->print matrix info", printMatInfo);
-    }
+    PetscSolverGlobalMatrix(string name);
 
     void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
 
@@ -63,13 +56,9 @@ namespace AMDiS {
     void destroyVectorData();
 
   protected:
-    void removeDirichletRows(Matrix<DOFMatrix*> *seqMat, 
-			     ParallelDofMapping &dofMap,
-			     Mat mpiMat);
+    void removeDirichletRows(Matrix<DOFMatrix*> *seqMat);
 
-    void removeDirichletRows(SystemVector *seqVec, 
-			     ParallelDofMapping &dofMap,
-			     Vec mpiVec);
+    void removeDirichletRows(SystemVector *seqVec);
 
     /// Reads field split information and creats a splitting based on 
     /// component numbers.