Whatever.

AMDiS/src/DOFMatrix.cc

@@ -74,8 +74,6 @@ namespace AMDiS {
     elementMatrix.change_dim(nRow, nCol);
     rowIndices.resize(nRow);
     colIndices.resize(nCol);
-
-    applyDBCs.clear();
   }
 
 
@@ -239,11 +237,10 @@ namespace AMDiS {
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
 	  if (dofMap->isRankDof(rowIndices[i])) {
- 	    applyDBCs.insert(row);
-	    //	    dirichletDofs.push_back(row);
+	    dirichletDofs.insert(row);
 	  }
 #else
- 	  applyDBCs.insert(row);
+	  dirichletDofs.insert(row);
 #endif
 	}
       } else {
@@ -489,20 +486,17 @@ namespace AMDiS {
   }
 
 
-  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
+  void DOFMatrix::clearDirichletRows()
   {      
-    FUNCNAME("DOFMatrix::removeRowsWithDBC()");
+    FUNCNAME("DOFMatrix::clearDirichletRows()");
 
     // Do the following only in sequential code. In parallel mode, the specific
     // solver method must care about dirichlet boundary conditions.
 
-    //#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
     inserter_type &ins = *inserter;  
-    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
+    for (std::set<int>::iterator it = dirichletDofs.begin(); 
+	 it != dirichletDofs.end(); ++it)
       ins[*it][*it] = 1.0;    
-
-    rows->clear();
-    //#endif
   }
 
 
@@ -521,6 +515,8 @@ namespace AMDiS {
     }
 
     inserter = new inserter_type(matrix, nnz_per_row);
+
+    dirichletDofs.clear();
   }
 
 

AMDiS/src/DOFMatrix.h

@@ -242,12 +242,10 @@ namespace AMDiS {
     /// data compressing procedures.
     void finishAssembling();
 
-    /** \brief
-     * Enable insertion for assembly. You can optionally give an upper limit for
-     * entries per row (per column for CCS matrices).  Choosing this parameter
-     * too small can induce perceivable overhead for compressed matrices.  Thus,
-     * it's better to choose a bit too large than too small.
-     */
+    /// Enable insertion for assembly. You can optionally give an upper limit for
+    /// entries per row (per column for CCS matrices).  Choosing this parameter
+    /// too small can induce perceivable overhead for compressed matrices. Thus,
+    /// it's better to choose a bit too large than too small.
     void startInsertion(int nnz_per_row = 10);
 
     /// Finishes insertion. For compressed matrix types, this is where the
@@ -326,7 +324,7 @@ namespace AMDiS {
 			   int irow, int jcol, double entry,
 			   bool add = true);
 
-    void removeRowsWithDBC(std::set<int> *rows);
+    void clearDirichletRows();
 
     /// Prints \ref matrix to stdout
     void print() const;
@@ -362,12 +360,6 @@ namespace AMDiS {
       return boundaryManager; 
     }
 
-    /// Returns a pointer to \ref applyDBCs.
-    std::set<int>* getApplyDBCs() 
-    {
-      return &applyDBCs;
-    }
-
     inline void setBoundaryManager(BoundaryManager *bm) 
     {
       boundaryManager = bm;
@@ -459,8 +451,8 @@ namespace AMDiS {
     /// correspond to a dof at a dirichlet boundary, are ignored and the row is
     /// left blank. After assembling, the diagonal element of the matrix must be
     /// set to 1. The indices of all rows, where this should be done, are stored
-    ///  in this set.
-    std::set<int> applyDBCs;
+    /// in this set.
+    std::set<DegreeOfFreedom> dirichletDofs;
 
     /// Number of non zero entries per row (average). For instationary problems this
     /// information may be used in the next timestep to accelerate insertion of
@@ -472,8 +464,6 @@ namespace AMDiS {
      /// row FE spaces are owned by the rank or not. This is used to ensure that 
      /// Dirichlet BC is handled correctly in parallel computations.
      FeSpaceDofMap *dofMap;
-
-     //     std::set<DegreeOfFreedom> dirichletDofs;
 #endif
 
     /// Inserter object: implemented as pointer, allocated and deallocated as needed

AMDiS/src/DirichletBC.cc

@@ -64,9 +64,26 @@ namespace AMDiS {
     const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
 
     for (int i = 0; i < nBasFcts; i++) {
+
+#if 1
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
       if (vector->isRankDof(dofIndices[i]))
 #endif
+	if (localBound[i] == boundaryType) {
+	  double value = 0.0;
+	  if (f) {
+	    elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords);
+	    value = (*f)(worldCoords);
+	  } else {
+	    if (dofVec)
+	      value = (*dofVec)[dofIndices[i]];
+	    else
+	      ERROR_EXIT("There is something wrong!\n");
+	  }
+
+	  (*vector)[dofIndices[i]] = value;
+	}
+#else
       if (localBound[i] == boundaryType) {
 	double value = 0.0;
 	if (f) {
@@ -76,13 +93,16 @@ namespace AMDiS {
 	if (dofVec)
 	  value = (*dofVec)[dofIndices[i]];
 
-	//#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-	//	vector->setDirichletDofValue(dofIndices[i], value);
-	//#else
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+	vector->setDirichletDofValue(dofIndices[i], value);
+#else
 	(*vector)[dofIndices[i]] = value;
-	//#endif
+#endif
       }
+#endif
+
     }
+
   }
 
 }

AMDiS/src/ProblemStat.cc

@@ -1030,7 +1030,7 @@ namespace AMDiS {
 	assembledMatrix[i][j] = true;
 
 	if (assembleMatrix) {
-	  matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+	  matrix->clearDirichletRows();
 	  matrix->finishInsertion();
 	}
  	if (assembleMatrix && matrix->getBoundaryManager())
@@ -1289,7 +1289,7 @@ namespace AMDiS {
 	if (!matrix)
 	  continue;
 
-	matrix->removeRowsWithDBC(matrix->getApplyDBCs());      
+	matrix->clearDirichletRows();
 	matrix->finishInsertion();
 
  	if (matrix->getBoundaryManager())
@@ -1674,7 +1674,7 @@ namespace AMDiS {
     // == private matrix and vector to the global one.                         ==
 
     if (matrix)
-      matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+      matrix->clearDirichletRows();
 
     if (matrix)
       matrix->finishAssembling();

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -822,35 +822,36 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
 
-    ERROR_EXIT("DO NOT CALL!\n");
-
+#if 0
     vector<int> dofsInterior, dofsCoarse;
 
     int nComponents = mat->getNumRows();
     for (int i = 0; i < nComponents; i++) {
-      if ((*mat)[i][i]) {
-	const FiniteElemSpace *feSpace = (*mat)[i][i]->getRowFeSpace();
-	
-	std::set<DegreeOfFreedom> &dirichletDofs = *((*mat)[i][i]->getApplyDBCs());
-	
-	MSG("DIRICHLET DOFS: %d %d -> %d\n", i, i, dirichletDofs.size());
-	
-	for (std::set<DegreeOfFreedom>::iterator it = dirichletDofs.begin();
-	     it != dirichletDofs.end(); ++it) {
-	  if (isCoarseSpace(feSpace, *it)) {
-	    if ((*coarseSpaceMap)[feSpace].isRankDof(*it)) {
-	      int globalDof = (*coarseSpaceMap)[feSpace][*it].global;
-	      dofsCoarse.push_back(coarseSpaceMap->getMatIndex(i, globalDof));
-	    }
-	  } else {
-	    if ((*interiorMap)[feSpace].isRankDof(*it)) {
-	      int globalDof = (*interiorMap)[feSpace][*it].global;
-	      dofsInterior.push_back(interiorMap->getMatIndex(i, globalDof));
+      for (int j = 0; j < nComponents; j++) {
+	if ((*mat)[i][j]) {
+	  const FiniteElemSpace *feSpace = (*mat)[i][j]->getRowFeSpace();
+	    
+	  std::set<DegreeOfFreedom> &dirichletDofs = *((*mat)[i][j]->getApplyDBCs());
+
+	  MSG("DIRICHLET DOFS: %d %d -> %d\n", i, j, dirichletDofs.size());
+	  
+	  for (std::set<DegreeOfFreedom>::iterator it = dirichletDofs.begin();
+	       it != dirichletDofs.end(); ++it) {
+	    if (isCoarseSpace(feSpace, *it)) {
+	      if ((*coarseSpaceMap)[feSpace].isRankDof(*it)) {
+		int globalDof = (*coarseSpaceMap)[feSpace][*it].global;
+		dofsCoarse.push_back(coarseSpaceMap->getMatIndex(i, globalDof));
+	      }
+	    } else {
+	      if ((*interiorMap)[feSpace].isRankDof(*it)) {
+		int globalDof = (*interiorMap)[feSpace][*it].global;
+		dofsInterior.push_back(interiorMap->getMatIndex(i, globalDof));
+	      }
 	    }
 	  }
+	} else {
+	  MSG("NO MAT DIAG in %d\n", i);
 	}
-      } else {
-	MSG("NO MAT DIAG in %d\n", i);
       }
     }
 
@@ -860,6 +861,7 @@ namespace AMDiS {
     if (coarseSpaceMap != NULL)
       MatZeroRows(matCoarseCoarse, dofsCoarse.size(), &(dofsCoarse[0]), 1.0, 
 		  PETSC_NULL, PETSC_NULL);
+#endif
   }
 
 
@@ -867,8 +869,6 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
 
-    ERROR_EXIT("DO NOT CALL!\n");
-
     int nComponents = vec->getSize();
     for (int i = 0; i < nComponents; i++) {
       const FiniteElemSpace *feSpace = vec->getDOFVector(i)->getFeSpace();