* This and that

AMDiS/compositeFEM/src/SubElementAssembler.cc

@@ -74,8 +74,6 @@ namespace AMDiS {
 						const ElInfo *elInfo, 
 						ElementVector *userVec)
   {
-    double corrFactor;
-
     /** 
      * Manipulate the quadratures of the SubAssemblers for subelement.
      */
@@ -88,7 +86,7 @@ namespace AMDiS {
      * determinant of element has been used instead of the determinant of subelement. Thus
      * the result must be corrected with respect to subelement.
      */
-    corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
+    double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
     for (int i = 0; i < nRow; i++) {
       (*userVec)[i] *= corrFactor;
     }

AMDiS/compositeFEM/src/SubElementAssembler.h

@@ -57,8 +57,8 @@ namespace AMDiS {
     MEMORY_MANAGED(SubElementAssembler);
 
     SubElementAssembler(Operator *op, 
-			const FiniteElemSpace *rowFESpace_,
-			const FiniteElemSpace *colFESpace_=NULL);
+			const FiniteElemSpace *rowFESpace,
+			const FiniteElemSpace *colFESpace = NULL);
 
     virtual ~SubElementAssembler()
     {

AMDiS/src/DOFMatrix.cc

@@ -40,7 +40,7 @@ namespace AMDiS {
     }
 
     if (rowFESpace && rowFESpace->getAdmin())
-      (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
+      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
 
     boundaryManager = NEW BoundaryManager;
     elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
@@ -284,9 +284,12 @@ namespace AMDiS {
   {
     int j;
 
-    for (j=0; j<static_cast<int>(matrix[a].size());j++) 
-      if (b==matrix[a][j].col) break;
-    if (j!=static_cast<int>(matrix[a].size())) matrix[a][j].col=c;
+    for (j = 0; j<static_cast<int>(matrix[a].size());j++) 
+      if (b == matrix[a][j].col) 
+	break;
+
+    if (j != static_cast<int>(matrix[a].size())) 
+      matrix[a][j].col = c;
   }
 
   void DOFMatrix::addSparseDOFEntry(double sign, int irow, 
@@ -386,12 +389,13 @@ namespace AMDiS {
       }
     }
     int usedSize = rowFESpace->getAdmin()->getUsedSize();
-    for(i=0; i < usedSize; i++) {
+    for (i = 0; i < usedSize; i++) {
       row = reinterpret_cast< std::vector<MatEntry>*>(&(matrix[i])); 
       int rowSize = static_cast<int>(row->size());
-      for (j=0; j < rowSize; j++) {
+      for (j = 0; j < rowSize; j++) {
 	col = (*row)[j].col;
-	if (entryUsed(i,j)) (*row)[j].col = newDOF[col];
+	if (entryUsed(i,j)) 
+	  (*row)[j].col = newDOF[col];
       }    
     }
   }

AMDiS/src/DOFMatrix.h

@@ -211,7 +211,7 @@ namespace AMDiS {
      */
     DOFMatrix(const FiniteElemSpace* rowFESpace,
 	      const FiniteElemSpace* colFESpace,
-	      std::string n="");
+	      std::string n = "");
 
     /** \brief
      * Copy-Constructor
@@ -267,18 +267,22 @@ namespace AMDiS {
      * Returns true if matrix[a][b] has entry \ref NO_MORE_ENTRIES
      */
     inline bool noMoreEntries(DegreeOfFreedom a,int b) const {
-      return  (NO_MORE_ENTRIES==((matrix[a])[b]).col);
+      return (NO_MORE_ENTRIES == ((matrix[a])[b]).col);
     };
 
     /** \brief
      * Returns \ref coupleMatrix.
      */
-    inline bool isCoupleMatrix() { return coupleMatrix; };
+    inline bool isCoupleMatrix() { 
+      return coupleMatrix; 
+    };
 
     /** \brief
      * Returns \ref coupleMatrix.
      */
-    inline void setCoupleMatrix(bool c) { coupleMatrix = c; };
+    inline void setCoupleMatrix(bool c) { 
+      coupleMatrix = c; 
+    };
 
     /** \brief
      * Matrix-matrix multiplication.
@@ -437,6 +441,14 @@ namespace AMDiS {
       return matrix.size();
     };
 
+    /** \brief
+     * Returns the number of used rows (equal to number of used DOFs in
+     * the row FE space).
+     */
+    inline int getUsedSize() const {
+      return rowFESpace->getAdmin()->getUsedSize();
+    };
+
     /** \brief
      * Returns number of cols. For that, the function iteratos over all
      * rows and searchs for the entry with the highest col number.
@@ -453,7 +465,8 @@ namespace AMDiS {
 	}
       }
 
-      return max;
+      // Add one to the maximum value, because the indeces start with 0.
+      return max + 1;
     };
 
     /** \brief

AMDiS/src/ProblemVec.cc

@@ -12,6 +12,7 @@
 #include "FileWriter.h"
 #include "CoarseningManager.h"
 #include "RefinementManager.h"
+#include "DualTraverse.h"
 #include "Mesh.h"
 #include "OEMSolver.h"
 #include "Preconditioner.h"
@@ -271,8 +272,6 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemVec::createMatricesAndVectors()");
 
-    int i;
-
     // === create vectors and system matrix ===
 
     systemMatrix_ = NEW Matrix<DOFMatrix*>(nComponents, nComponents);
@@ -282,7 +281,7 @@ namespace AMDiS {
 
     char number[10];
     std::string numberedName;
-    for (i = 0; i < nComponents; i++) {
+    for (int i = 0; i < nComponents; i++) {
       (*systemMatrix_)[i][i] = NEW DOFMatrix(componentSpaces[i], 
 					     componentSpaces[i], "A_ii");
       (*systemMatrix_)[i][i]->setCoupleMatrix(false);
@@ -700,8 +699,6 @@ namespace AMDiS {
 #pragma omp parallel for 
 #endif
     for (i = 0; i < nComponents; i++) {
-      const BasisFunction *basisFcts = componentSpaces[i]->getBasisFcts();
-
       for (int j = 0; j < nComponents; j++) {
 	// Only if this variable is true, the current matrix will be assembled.	
 	bool assembleMatrix = true;
@@ -729,43 +726,21 @@ namespace AMDiS {
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->initMatrix(matrix);
 
-	BoundaryType *bound = NULL;
-	if (useGetBound_) {
-	  bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
-	}
-	
-	TraverseStack stack;
-	ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
-	
-	while (elInfo) {
-	  if (useGetBound_) {
-	    basisFcts->getBound(elInfo, bound);
-	  }
-	  
-	  if (assembleMatrix) {
-	    matrix->assemble(1.0, elInfo, bound);
-	    
-	    if (matrix->getBoundaryManager()) {
-	      matrix->
-		getBoundaryManager()->
-		fillBoundaryConditions(elInfo, matrix);
-	    }		      
-	  }
-	  
-	  if (i == j) {
-	    rhs_->getDOFVector(i)->assemble(1.0, elInfo, bound);
-	  }
-	  
-	  elInfo = stack.traverseNext(elInfo);
+	if (componentSpaces[i] == componentSpaces[j]) {
+	  assembleOnOneMesh(componentSpaces[i],
+			    assembleFlag,
+			    assembleMatrix ? matrix : NULL,
+			    (i == j) ? rhs_->getDOFVector(i) : NULL);
+	} else {
+	  assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
+			      assembleFlag,
+			      assembleMatrix ? matrix : NULL,
+			      (i == j) ? rhs_->getDOFVector(i) : NULL);	  
 	}
-	
+
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	  
 	
-	if (useGetBound_) {
-	  FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
-	}	  
-	
 	assembledMatrix_[i][j] = true;
       }
 
@@ -956,6 +931,112 @@ namespace AMDiS {
 	addBoundaryCondition(periodic);
   }
 
+  void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace, Flag assembleFlag,
+				     DOFMatrix *matrix, DOFVector<double> *vector)
+  {
+    Mesh *mesh = feSpace->getMesh();
+    const BasisFunction *basisFcts = feSpace->getBasisFcts();
+
+    BoundaryType *bound = NULL;
+    if (useGetBound_) {
+      bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
+    }
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+    
+    while (elInfo) {
+      if (useGetBound_) {
+	basisFcts->getBound(elInfo, bound);
+      }
+      
+      if (matrix) {
+	matrix->assemble(1.0, elInfo, bound);
+	
+	if (matrix->getBoundaryManager()) {
+	  matrix->getBoundaryManager()->
+	    fillBoundaryConditions(elInfo, matrix);
+	}		      
+      }
+      
+      if (vector) {
+	vector->assemble(1.0, elInfo, bound);
+      }
+      
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    if (useGetBound_) {
+      FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
+    }	     
+  }
+
+  void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace, FiniteElemSpace *colFeSpace,
+				       Flag assembleFlag,
+				       DOFMatrix *matrix, DOFVector<double> *vector)
+  {
+    Mesh *rowMesh = rowFeSpace->getMesh();
+    Mesh *colMesh = colFeSpace->getMesh();
+
+    const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();
+    BoundaryType *bound = NULL;
+    if (useGetBound_) {
+      bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
+    }
+
+    DualTraverse dualTraverse;
+    ElInfo *rowElInfo, *colElInfo;
+    ElInfo *largeElInfo, *smallElInfo;
+
+    bool cont = dualTraverse.traverseFirst(rowMesh, colMesh, -1, -1,
+					   assembleFlag, assembleFlag,
+					   &rowElInfo, &colElInfo,
+					   &smallElInfo, &largeElInfo);
+    while (cont) {
+      Element *rowElem = rowElInfo->getElement();
+      Element *colElem = colElInfo->getElement();
+
+      if (rowElInfo->getLevel() != colElInfo->getLevel()) {
+	if (smallElInfo == rowElInfo) 
+	  std::cout << "Row = small\n";
+	else
+	  std::cout << "Row = large\n";
+
+	if (largeElInfo == colElInfo) 
+	  std::cout << "Col = large\n";
+	else
+	  std::cout << "Col = small\n";
+
+
+	ERROR_EXIT("KOMMT GLEICH!\n");
+      }
+
+      if (useGetBound_) {
+	basisFcts->getBound(rowElInfo, bound);
+      }
+      
+      if (matrix) {
+	matrix->assemble(1.0, rowElInfo, bound);
+	
+	if (matrix->getBoundaryManager()) {
+	  matrix->getBoundaryManager()->
+	    fillBoundaryConditions(rowElInfo, matrix);
+	}		      
+      }
+      
+      if (vector) {
+	vector->assemble(1.0, rowElInfo, bound);
+      }
+
+      cont = dualTraverse.traverseNext(&rowElInfo, &colElInfo,
+				       &smallElInfo, &largeElInfo);
+    }
+
+    if (useGetBound_) {
+      FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
+    }
+  }
+
   void ProblemVec::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
   {
     FUNCNAME("ProblemVec::writeResidualMesh()");

AMDiS/src/ProblemVec.h

@@ -289,6 +289,17 @@ namespace AMDiS {
     inline void allowFirstRefinement() {
       allowFirstRef_ = true;
     }; 
+
+    /** \brief
+     * This function assembles a DOFMatrix and a DOFVector for the case,
+     * the meshes from row and col FE-space are equal.
+     */
+    void assembleOnOneMesh(FiniteElemSpace *feSpace, Flag assembleFlag,
+			   DOFMatrix *matrix, DOFVector<double> *vector);
+
+    void assembleOnDifMeshes(FiniteElemSpace *rowFeSpace, FiniteElemSpace *colFeSpace,
+			     Flag assembleFlag,
+			     DOFMatrix *matrix, DOFVector<double> *vector);
   
     // ===== getting-methods ======================================================