New handling of dirichlet bcs in dof matrix.

AMDiS/src/DOFMatrix.cc

@@ -147,8 +147,6 @@ namespace AMDiS {
     }
   }
 
-
-  // ok
   DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
   {
     TEST_EXIT(rhs.inserter == 0 && inserter == 0)("Cannot copy during insertion");
@@ -158,6 +156,7 @@ namespace AMDiS {
     operators = rhs.operators;
     operatorFactor = rhs.operatorFactor;
     matrix = rhs.matrix;
+    coupleMatrix = rhs.coupleMatrix;
     if (rhs.boundaryManager) {
       boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
     } else {
@@ -173,7 +172,6 @@ namespace AMDiS {
     return *this;
   }
 
-
   void DOFMatrix::addElementMatrix(double sign, 
 				   const ElementMatrix &elMat, 
 				   const BoundaryType *bound,
@@ -198,7 +196,7 @@ namespace AMDiS {
 
       if (condition && condition->isDirichlet()) {
 	if (!coupleMatrix)
-	  ins[row][row]= 1.0;
+	  applyDBCs.insert(static_cast<int>(row));
       } else
 	for (int j = 0; j < n_col; j++) {  // for all columns
 	  col = elMat.colIndices[j];
@@ -211,7 +209,7 @@ namespace AMDiS {
     }
   }
 
-  double DOFMatrix::logAcc(DegreeOfFreedom a,DegreeOfFreedom b) const
+  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
   {
     return matrix[a][b];
   }
@@ -366,20 +364,9 @@ namespace AMDiS {
     return fillFlag;
   }
 
-  void DOFMatrix::axpy(double a, 
-		       const DOFMatrix& x,
-		       const DOFMatrix& y)
+  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
   {
-    FUNCNAME("DOFMatrix::axpy");
-    TEST_EXIT(x.getRowFESpace() == y.getRowFESpace() &&
-	      rowFESpace == x.getRowFESpace())
-      ("row fe-spaces not equal\n");
-    TEST_EXIT(x.getColFESpace() == y.getColFESpace() &&
-	      colFESpace == x.getColFESpace())
-      ("col fe-spaces not equal\n");
-
     matrix+= a * x.matrix + y.matrix;
-
   }
 
 
@@ -402,21 +389,12 @@ namespace AMDiS {
 
   void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
   {      
-      ERROR_EXIT("TODO: removeRowsWithDBC");
-
-//     for (std::set<int>::iterator it = rows->begin();
-// 	 it != rows->end(); 
-// 	 ++it) {
-//       if (coupleMatrix) {
-// 	matrix[*it].resize(0);
-//       } else {
-// 	matrix[*it].resize(1);
-// 	matrix[*it][0].col = *it;
-// 	matrix[*it][0].entry = 1.0;
-//       }
-//     }
-
-//     rows->clear();
+    inserter_type &ins= *inserter;
+   
+    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
+      ins[*it][*it] = 1.0;
+
+    rows->clear();
   }
 
   void DOFMatrix::createPictureFile(const char* filename, int dim)

AMDiS/src/DOFMatrix.h

@@ -256,38 +256,26 @@ namespace AMDiS {
     const bool& operator[](int i) const {ERROR_EXIT("Shouldn't be used, only fake."); return dummy;}
 #endif
  
-    /** \brief
-     * DOFMatrix does not need to be compressed before assembling, when using MTL4.
-     */
+    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
     void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
 
-    /** \brief
-     * Implements DOFIndexedBase::freeDOFContent()
-     */ 
+    /// Implements DOFIndexedBase::freeDOFContent()
     virtual void freeDOFContent(int index);
 
-    /** \brief
-     * Returns \ref coupleMatrix.
-     */
+    /// Returns \ref coupleMatrix.
     inline bool isCoupleMatrix() { 
       return coupleMatrix; 
     }
 
-    /** \brief
-     * Returns \ref coupleMatrix.
-     */
+    /// Returns \ref coupleMatrix.
     inline void setCoupleMatrix(bool c) { 
       coupleMatrix = c; 
     }
 
-      /** \brief
-     * a*x + y
-     */
+    /// a * x + y
     void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
 
-    /** \brief
-     * Multiplication with a scalar.
-     */
+    /// Multiplication with a scalar.
     void scal(double s);
 
     /** \brief
@@ -427,7 +415,7 @@ namespace AMDiS {
 
     /// Returns number of rows (\ref matrix.size())
     inline int getSize() const { 
-	return num_rows(matrix);
+      return num_rows(matrix);
     }
 
     /** \brief
@@ -438,7 +426,6 @@ namespace AMDiS {
       return rowFESpace->getAdmin()->getUsedSize();
     }
 
-
     // Only fake, shouldn't be called
     /** \brief
      * Returns number of cols. For that, the function iteratos over all
@@ -510,6 +497,7 @@ namespace AMDiS {
       return boundaryManager; 
     }
 
+    /// Returns a pointer to \ref applyDBCs.
     std::set<int>* getApplyDBCs() {
       return &applyDBCs;
     }
@@ -621,7 +609,7 @@ namespace AMDiS {
     std::string name;
 
     /// Sparse matrix, type is a template parameter by default compressed2D<double>
-    base_matrix_type         matrix;
+    base_matrix_type matrix;
 
     /// Used while mesh traversal
     static DOFMatrix *traversePtr;
@@ -644,13 +632,22 @@ namespace AMDiS {
     ///
     BoundaryManager *boundaryManager;
 
-    ///
+    /** \brief
+     * If false, the matrix is a diagonal matrix within a matrix of DOF matrices.
+     * Otherwise the value is true, and the matrix is an off-diagonal matrix.
+     */
     bool coupleMatrix;
 
     /// Temporary variable used in assemble()
     ElementMatrix *elementMatrix;
 
-    ///
+    /* \brief
+     * A set of row indices. When assembling the DOF matrix, all rows, that
+     * 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;
 
 #ifdef HAVE_PARALLEL_AMDIS
@@ -659,7 +656,7 @@ namespace AMDiS {
 #endif
 
     /// Inserter object: implemented as pointer, allocated and deallocated as needed
-    inserter_type      *inserter;
+    inserter_type *inserter;
       
     friend class DOFAdmin;
     friend class DOFVector<double>;