* Bug fix

AMDiS/src/Assembler.cc

@@ -1543,11 +1543,8 @@ namespace AMDiS {
     }
 
     elMat->set(0.0);
-      
-    Element *element = elInfo->getElement();
-
-    
-    rowFESpace->getBasisFcts()->getLocalIndicesVec(element,
+         
+    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
 						   rowFESpace->getAdmin(),
 						   &(elMat->rowIndices));
 
@@ -1555,7 +1552,7 @@ namespace AMDiS {
     if (rowFESpace == colFESpace) {
       elMat->colIndices = elMat->rowIndices;
     } else {
-      colFESpace->getBasisFcts()->getLocalIndicesVec(element,
+      colFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
 						     colFESpace->getAdmin(),
 						     &(elMat->colIndices));
     }

AMDiS/src/DOFMatrix.cc

@@ -43,10 +43,12 @@ namespace AMDiS {
       (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
 
     boundaryManager = NEW BoundaryManager;
+    elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
+				      colFESpace->getBasisFcts()->getNumber());
   }
 
   DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
-    : name(rhs.name+"copy")
+    : name(rhs.name + "copy")
   {
     *this = rhs;
     if (rowFESpace && rowFESpace->getAdmin())
@@ -60,13 +62,14 @@ namespace AMDiS {
     if (rowFESpace && rowFESpace->getAdmin()) {
       (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
     }  
+
+    if (elementMatrix)
+      DELETE elementMatrix;
   }
 
   void DOFMatrix::print() const
   {
     FUNCNAME("DOFMatrix::print()");
-    int  i, j, jcol;
-    DOFMatrix::MatrixRow row;
 
     int sizeUsed = rowFESpace->getAdmin()->getUsedSize();
 
@@ -75,19 +78,18 @@ namespace AMDiS {
       return;
     }
 
-    for (i = 0; i < sizeUsed; i++) { 
-      row = matrix[i];
+    for (int i = 0; i < sizeUsed; i++) { 
+      DOFMatrix::MatrixRow row = matrix[i];
       MSG("row %3d:",i);
       int rowSize = static_cast<int>( row.size());
-      for (j = 0; j < rowSize; j++) {
-	jcol = row[j].col;
-	if (entryUsed(i,j)) {
+      for (int j = 0; j < rowSize; j++) {
+	int jcol = row[j].col;
+	if (entryUsed(i, j)) {
 	  Msg::print(" (%3d,%20.17lf)", jcol, row[j].entry);
 	}
       }
       Msg::print("\n");
     }
-    return;
   }
 
   void DOFMatrix::printRow(int i) const
@@ -216,7 +218,13 @@ namespace AMDiS {
     if (rhs.boundaryManager) {
       boundaryManager = new BoundaryManager(*rhs.boundaryManager);
     } else {
-      boundaryManager=NULL;
+      boundaryManager = NULL;
+    }
+    if (rhs.elementMatrix) {
+      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
+					colFESpace->getBasisFcts()->getNumber());
+    } else {
+      elementMatrix = NULL;
     }
 
     return *this;
@@ -230,13 +238,12 @@ namespace AMDiS {
   {
     FUNCNAME("DOFMatrix::addElementMatrix");
 
-    DegreeOfFreedom row, col;
-    double entry;
+    DegreeOfFreedom row;
 
-    int n_row = elMat.rowIndices.getSize();
-    int n_col = elMat.colIndices.getSize();
+    int nRow = elMat.rowIndices.getSize();
+    int nCol = elMat.colIndices.getSize();
 
-    for (int i = 0; i < n_row; i++)  {   // for all rows of element matrix
+    for (int i = 0; i < nRow; i++)  {   // for all rows of element matrix
       row = elMat.rowIndices[i];
 
       BoundaryCondition *condition = 
@@ -252,10 +259,9 @@ namespace AMDiS {
 	  ((*matrixRow)[0]).entry = 1.0;
 	}     
       } else {     
-	for (int j = 0; j < n_col; j++) {  // for all columns
-	  col = elMat.colIndices[j];
-	  entry = elMat[i][j];
-	  addSparseDOFEntry(sign, row, col, entry, add); 
+	for (int j = 0; j < nCol; j++) {  // for all columns
+	  addSparseDOFEntry(sign, row, elMat.colIndices[j], 
+			    elMat[i][j], add); 
 	}
       }
     }
@@ -310,19 +316,18 @@ namespace AMDiS {
       return;
 
     int freeCol = -1;
-    int rowSize = static_cast<int>( row->size());
+    int rowSize = static_cast<int>(row->size());
 
     TEST_EXIT_DBG(jcol >= 0 && 
 		  jcol < colFESpace->getAdmin()->getUsedSize())
       ("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
 
     // first entry is diagonal entry
-    if (rowFESpace == colFESpace)
-      if (rowSize == 0) {
-	MatEntry newEntry = {irow, 0.0};
-	row->push_back(newEntry);
-	rowSize = 1;
-      }
+    if ((rowSize == 0) && (rowFESpace == colFESpace)) {
+      MatEntry newEntry = {irow, 0.0};
+      row->push_back(newEntry);
+      rowSize = 1;
+    }
 
     int i;
     // search jcol
@@ -339,7 +344,7 @@ namespace AMDiS {
       if ((*row)[i].col == NO_MORE_ENTRIES) {
 	freeCol = i;
 	if (rowSize > i+1) {
-	  (*row)[i+1].entry = NO_MORE_ENTRIES;
+	  (*row)[i + 1].entry = NO_MORE_ENTRIES;
 	}
 	break;
       }
@@ -347,8 +352,10 @@ namespace AMDiS {
 
     // jcol found?
     if (i < rowSize) {
-      if (!add) 
+      if (!add) {
 	(*row)[i].entry = 0.0;
+      }
+
       (*row)[i].entry += sign * entry;
     } else {
       if (freeCol == -1) {
@@ -449,9 +456,9 @@ namespace AMDiS {
     
     Operator *operat = op ? op : operators[0];
 
-    // Do no combine the next two lines into one!
-    ElementMatrix *elementMatrix = NULL;
-    elementMatrix = operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+    // !!! Do no combine the next two lines into one !!!
+    //    ElementMatrix *elementMatrix = NULL;
+    operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
     
     if (op) {
       op->getElementMatrix(elInfo, elementMatrix);
@@ -470,8 +477,7 @@ namespace AMDiS {
     }
    
     addElementMatrix(factor, *elementMatrix, bound);   
-
-    DELETE elementMatrix;
+    //    DELETE elementMatrix;
   }
 
   Flag DOFMatrix::getAssembleFlag()

AMDiS/src/DOFMatrix.h

@@ -671,6 +671,11 @@ namespace AMDiS {
      */
     bool coupleMatrix;
 
+    /** \brief
+     * Temporary variable used in assemble()
+     */
+    ElementMatrix *elementMatrix;
+
     friend class DOFAdmin;
     friend class DOFVector<double>;
     friend class DOFVector<unsigned char>;