-> juropa.

AMDiS/src/DOFMatrix.cc

@@ -280,7 +280,8 @@ namespace AMDiS {
     if (factor != 1.0)
       elementMatrix *= factor;
 
-    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
+    if (operators.size())
+      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
   }
 
 
@@ -289,17 +290,17 @@ namespace AMDiS {
 			   const BoundaryType *bound,
 			   Operator *op)
   {
-      FUNCNAME("DOFMatrix::assemble()");
-
-      TEST_EXIT_DBG(op)("No operator!\n");
-
-      set_to_zero(elementMatrix);
-      op->getElementMatrix(elInfo, elementMatrix, factor);
-
-      if (factor != 1.0)
+    FUNCNAME("DOFMatrix::assemble()");
+    
+    TEST_EXIT_DBG(op)("No operator!\n");
+    
+    set_to_zero(elementMatrix);
+    op->getElementMatrix(elInfo, elementMatrix, factor);
+    
+    if (factor != 1.0)
 	elementMatrix *= factor;
-
-      addElementMatrix(elementMatrix, bound, elInfo, NULL);
+    
+    addElementMatrix(elementMatrix, bound, elInfo, NULL);
   }
 
 
@@ -334,7 +335,8 @@ namespace AMDiS {
     if (factor != 1.0)
       elementMatrix *= factor;
 
-    addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
+    if (op || operators.size())
+      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
   }
 
 

AMDiS/src/Operator.h

@@ -261,7 +261,8 @@ namespace AMDiS {
 	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, LALt);
     }
   
-    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the results to Lb.
+    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the 
+    /// results to Lb.
     void getLbGrdPsi(const ElInfo *elInfo, 
 		     std::vector<mtl::dense_vector<double> >& Lb) const
     {      
@@ -270,7 +271,8 @@ namespace AMDiS {
 	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, Lb);
     }
 
-    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the results to Lb.
+    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the 
+    /// results to Lb.
     void getLbGrdPhi(const ElInfo *elInfo, 
 		     std::vector<mtl::dense_vector<double> >& Lb) const
     {      
@@ -293,13 +295,15 @@ namespace AMDiS {
       return secondOrder.size() != 0;
     }
 
-    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
+    /// Returns true, if there are first order terms (grdPsi). Returns 
+    /// false otherwise.
     inline bool firstOrderTermsGrdPsi() 
     {
       return firstOrderGrdPsi.size() != 0;
     }
 
-    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
+    /// Returns true, if there are first order terms (grdPhi). Returns 
+    /// false otherwise.
     inline bool firstOrderTermsGrdPhi() 
     {
       return firstOrderGrdPhi.size() != 0;
@@ -340,11 +344,9 @@ namespace AMDiS {
     /// If true, the operator needs a dual traverse over two meshes for assembling.
     bool needDualTraverse;
 
-    /** \brief
-     * Calculates the element matrix and/or the element vector. It is
-     * created especially for this Operator, when \ref getElementMatrix()
-     * or \ref getElementVector is called for the first time.
-     */
+    /// Calculates the element matrix and/or the element vector. It is
+    /// created especially for this Operator, when \ref getElementMatrix()
+    /// or \ref getElementVector is called for the first time.
     Assembler* assembler;
 
     /// List of all second order terms
@@ -359,10 +361,9 @@ namespace AMDiS {
     /// List of all zero order terms
     vector<OperatorTerm*> zeroOrder;
 
-    /** \brief
-     * Pointer to the solution of the last timestep. Can be used for a more efficient
-     * assemblage of the element vector when the element matrix was already computed.
-     */
+    /// Pointer to the solution of the last timestep. Can be used for a more
+    /// efficient assemblage of the element vector when the element matrix 
+    /// was already computed.
     const DOFVectorBase<double> *uhOld;
 
     /// Spezifies whether optimized assemblers are used or not.

AMDiS/src/ProblemStat.cc

@@ -95,6 +95,8 @@ namespace AMDiS {
     for (int i = 0; i < nComponents; i++)
       Parameters::get(name + "->name[" + lexical_cast<string>(i) + "]",
 		      componentNames[i]);
+
+    Parameters::get("debug->write asm info", writeAsmInfo);
   }
 
   void ProblemStatSeq::initialize(Flag initFlag,
@@ -757,26 +759,28 @@ namespace AMDiS {
     // Used to calculate the overall number of non zero entries.
     int nnz = 0;
 
-    for (int i = 0; i < nComponents; i++) {
+    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
 
       MSG("%d DOFs for %s\n", 
-	  componentSpaces[i]->getAdmin()->getUsedDofs(), 
-	  componentSpaces[i]->getName().c_str());
+	  componentSpaces[rowComponent]->getAdmin()->getUsedDofs(), 
+	  componentSpaces[rowComponent]->getName().c_str());
 
-      rhs->getDOFVector(i)->set(0.0);
+      rhs->getDOFVector(rowComponent)->set(0.0);
 
-      for (int j = 0; j < nComponents; j++) {
+      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
 
 	if (writeAsmInfo)
-	  MSG("--------- %d %d -------------\n", i, j);
+	  MSG("--------- %d %d -------------\n", rowComponent, colComponent);
 
 	// Only if this variable is true, the current matrix will be assembled.	
 	bool assembleMatrix = true;
 	// The DOFMatrix which should be assembled (or not, if assembleMatrix
 	// will be set to false).
-	DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
+	DOFMatrix *matrix = 
+	  (asmMatrix ? (*systemMatrix)[rowComponent][colComponent] : NULL);
 
 	if (writeAsmInfo && matrix) {
+	  MSG("  -> matrix has %d operators\n", matrix->getOperators().size());
 	  for (vector<Operator*>::iterator it = matrix->getOperatorsBegin();
 	       it != matrix->getOperatorsEnd(); ++it) {
 	    Assembler *assembler = (*it)->getAssembler();
@@ -796,12 +800,13 @@ namespace AMDiS {
 	
 	// If the matrix was assembled before and it is marked to be assembled
 	// only once, it will not be assembled.
-	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
+	if (assembleMatrixOnlyOnce[rowComponent][colComponent] && 
+	    assembledMatrix[rowComponent][colComponent]) {
 	  assembleMatrix = false;
 	} else if (matrix) {
 	  matrix->getBaseMatrix().
-	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
-		       componentSpaces[j]->getAdmin()->getUsedSize());
+	    change_dim(componentSpaces[rowComponent]->getAdmin()->getUsedSize(), 
+		       componentSpaces[colComponent]->getAdmin()->getUsedSize());
 
 	  set_to_zero(matrix->getBaseMatrix());
 	}
@@ -813,7 +818,7 @@ namespace AMDiS {
 	// If the matrix should not be assembled, the rhs vector has to be considered.
 	// This will be only done, if i == j. So, if both is not true, we can jump
 	// to the next matrix.
-	if (!assembleMatrix && i != j) {
+	if (!assembleMatrix && rowComponent != colComponent) {
 	  if (matrix)
 	    nnz += matrix->getBaseMatrix().nnz();
 
@@ -825,12 +830,14 @@ namespace AMDiS {
 
 	// The simplest case: either the right hand side has no operaters, no aux
 	// fe spaces, or all aux fe spaces are equal to the row and col fe space.
-	assembleOnOneMesh(componentSpaces[i],
+	assembleOnOneMesh(componentSpaces[rowComponent],
 			  assembleFlag,
 			  assembleMatrix ? matrix : NULL,
-			  ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+			  ((rowComponent == colComponent) && asmVector) ? 
+			  rhs->getDOFVector(rowComponent) : 
+			  NULL);
       
-	assembledMatrix[i][j] = true;
+	assembledMatrix[rowComponent][colComponent] = true;
 
 	if (assembleMatrix)
 	  matrix->finishInsertion();
@@ -840,14 +847,13 @@ namespace AMDiS {
 
 	if (matrix)
 	  nnz += matrix->getBaseMatrix().nnz();
-      }
-      
+      }   
 
 
       // And now assemble boundary conditions on the vectors
-      assembleBoundaryConditions(rhs->getDOFVector(i),
- 				 solution->getDOFVector(i),
- 				 componentMeshes[i],
+      assembleBoundaryConditions(rhs->getDOFVector(rowComponent),
+ 				 solution->getDOFVector(rowComponent),
+ 				 componentMeshes[rowComponent],
  				 assembleFlag);     
     }  
 
@@ -1671,11 +1677,10 @@ namespace AMDiS {
     // == Finally, if we have assembled in parallel, we have to add the thread ==
     // == private matrix and vector to the global one.                         ==
 
-    if (matrix)
+    if (matrix) {
       matrix->clearDirichletRows();
-
-    if (matrix)
       matrix->finishAssembling();
+    }
 
     if (vector)
       vector->finishAssembling();

AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc

@@ -186,6 +186,7 @@ namespace AMDiS {
       int nRankCoarseRows = cMap->getRankDofs();
       int nOverallCoarseRows = cMap->getOverallDofs();
       
+      MSG("COARSE MAT %d SIZE %d %d\n", i + 1, nOverallCoarseRows, nOverallCoarseRows);
       MatCreateAIJ(mpiCommGlobal,
 		   nRankCoarseRows, nRankCoarseRows,
 		   nOverallCoarseRows, nOverallCoarseRows,
@@ -211,15 +212,16 @@ namespace AMDiS {
 	int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs());
 	int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs());
 	
+	MSG("COARSE MAT %d %d SIZE %d %d\n", i, j, nRowsOverallMat, nColsOverallMat);
 	MatCreateAIJ(mpiCommGlobal,
 		     nRowsRankMat, nColsRankMat,
 		     nRowsOverallMat, nColsOverallMat,
 		     0, nnz[i][j].dnnz, 0, nnz[i][j].onnz,
 		     &mat[i][j]);
-      MSG("REMOVE THIS!\n");
-      MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+	MSG("REMOVE THIS!\n");
+	MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 	  
-	
+	MSG("COARSE MAT %d %d SIZE %d %d\n", j, i, nRowsOverallMat, nColsOverallMat);
 	MatCreateAIJ(mpiCommGlobal,
 		     nColsRankMat, nRowsRankMat,
 		     nColsOverallMat, nRowsOverallMat,

AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h

@@ -132,10 +132,11 @@ namespace AMDiS {
     }
 
     /// Get the coarse space matrix of some system component.
-    inline Mat& getMatCoarseByComponent(int comp)
+    inline Mat& getMatCoarseByComponent(int rowComp, int colComp = -1)
     {
-      int matIndex = componentIthCoarseMap[comp] + 1;
-      return mat[matIndex][matIndex];
+      int rowMatIndex = componentIthCoarseMap[rowComp] + 1;
+      int colMatIndex = componentIthCoarseMap[(colComp == -1 ? rowComp : colComp)] + 1;
+      return mat[rowMatIndex][colMatIndex];
     }
 
     /// Get coupling matrix of the interior and the coarse space of a 

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -300,7 +300,6 @@ namespace AMDiS {
 
     if (fullInterface != NULL)
       interfaceDofMap.init(levelData, feSpaces, uniqueFe);
-      //      interfaceDofMap.init(levelData, fullInterface);
 
     if (fetiPreconditioner != FETI_NONE) {
       TEST_EXIT(meshLevel == 0)
@@ -531,10 +530,13 @@ namespace AMDiS {
 
     for (DofContainer::iterator it = allBoundaryDofs.begin();
 	 it != allBoundaryDofs.end(); ++it)
-      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
+      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) {
+	MSG("INSERT INTERFACE RANK DOF %d\n", **it);
 	interfaceDofMap[feSpace].insertRankDof(**it);
-      else
+      } else {
+	MSG("INSERT INTERFACE NON-RANK DOF %d\n", **it);
 	interfaceDofMap[feSpace].insertNonRankDof(**it);
+      }
   }
 
 
@@ -691,18 +693,10 @@ namespace AMDiS {
 
     int nLocalInterior = 0;    
     for (int i = 0; i < admin->getUsedSize(); i++) {
-      if (i == 0) {
-	MSG("JETZT: %d %d %d %d\n", admin->isDofFree(i), isPrimal(feSpace, i), isDual(feSpace, i), isInterface(feSpace, i));
-      }
-
       if (admin->isDofFree(i) == false && 
 	  isPrimal(feSpace, i) == false &&
 	  isDual(feSpace, i) == false &&
 	  isInterface(feSpace, i) == false) {
-	if (i == 0) {
-	  MSG("DRIN\n");
-	}
-
 	if (meshLevel == 0) {
 	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -136,6 +136,8 @@ namespace AMDiS {
 	if (!dofMat)
 	  continue;
 
+	MSG("=============== COMP %d %d ================= \n", rowComponent, colComponent);
+
 	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
 	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
 
@@ -161,6 +163,10 @@ namespace AMDiS {
 	    bool isColCoarse = 
 	      isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
 
+	    if (colComponent == 2 && col(*icursor) == 19) {
+	      MSG("HERE WE TEST: %d %d\n", isColCoarse, isRowCoarse);
+	    }
+
 	    if (isColCoarse) {
 	      if (isRowCoarse) {
 		cols.push_back(col(*icursor));
@@ -184,21 +190,26 @@ namespace AMDiS {
 	  // === Set matrix values. ===
 
 	  if (isRowCoarse) {
+	    MSG("A-GET MAT INDEX %d\n", rowComponent);
 	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
+	    MSG("B-GET MAT INDEX %d\n", colComponent);
 	    for (unsigned int i = 0; i < cols.size(); i++)
 	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
 
 	    // matCoarseCoarse
-	    MatSetValues(getMatCoarseByComponent(rowComponent),
+	    //MSG("SET COARSE MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
+	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
 			 1, &rowIndex, cols.size(),
 			 &(cols[0]), &(values[0]), ADD_VALUES);
 
 	    if (colsOther.size()) {
 	      if (subdomainLevel == 0) {
+		MSG("C0-GET MAT INDEX %d\n", colComponent);
 		for (unsigned int i = 0; i < colsOther.size(); i++)
 		  colsOther[i] = 
 		    interiorMap->getMatIndex(colComponent, colsOther[i]);
 	      } else {
+		MSG("C1-GET MAT INDEX %d\n", colComponent);
 		for (unsigned int i = 0; i < colsOther.size(); i++)
 		  colsOther[i] = 
 		    interiorMap->getMatIndex(colComponent, colsOther[i]) + 
@@ -206,38 +217,45 @@ namespace AMDiS {
 	      }
 
 	      // matCoarseInt
+	      //    MSG("SET COARSE-INT MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
 	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
 			   1, &rowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	    }
 	  } else {
+	    MSG("D-GET MAT INDEX %d\n", rowComponent);
 	    int localRowIndex = 
 	      (subdomainLevel == 0 ? 
 	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
 	       interiorMap->getMatIndex(rowComponent, *cursor));
 
+	    MSG("E-GET MAT INDEX %d\n", colComponent);
 	    for (unsigned int i = 0; i < cols.size(); i++) {
 	      if (subdomainLevel == 0)
 		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
 	      else
 		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
 	    }
+	    MSG("DONE\n");
 	    
   	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
   			 &(cols[0]), &(values[0]), ADD_VALUES);
-
+	    
 	    if (colsOther.size()) {
+	      MSG("F-GET MAT INDEX %d\n", colComponent);
 	      int globalRowIndex = interiorMap->getMatIndex(rowComponent, *cursor);
 
 	      if (subdomainLevel != 0)
 		globalRowIndex += rStartInterior;
-
+	      
+	      MSG("G-GET MAT INDEX %d\n", colComponent);
 	      for (unsigned int i = 0; i < colsOther.size(); i++)
 		colsOther[i] = 
 		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
 
 	      // matIntCoarse
-  	      MatSetValues(getMatInteriorCoarseByComponent(rowComponent), 
+	      //	      MSG("SET INT-COARSE MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
+  	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
 			   1, &globalRowIndex, colsOther.size(),
   			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
 	    }