blub bla

AMDiS/src/DOFMatrix.cc

@@ -534,12 +534,10 @@ namespace AMDiS {
 
     // 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 = dirichletDofs.begin(); 
 	 it != dirichletDofs.end(); ++it)
 	ins[*it][*it] = 1.0;
-#endif
   }
 
 

AMDiS/src/Element.h

@@ -438,9 +438,6 @@ namespace AMDiS {
 		       bool baseDofPtr = false,
 		       vector<GeoIndex>* dofGeoIndex = NULL) const = 0;
 
-    virtual void getSubBoundary(BoundaryObject bound, 
-				vector<BoundaryObject> &subBound) const = 0;
-
     /// Combines \ref getNodeDofs and \ref getHigherOrderDofs to one function. 
     /// See parameter description there.
     void getAllDofs(const FiniteElemSpace* feSpace, 
@@ -449,6 +446,8 @@ namespace AMDiS {
 		    bool baseDofPtr = false,
 		    vector<GeoIndex>* dofGeoIndex = NULL);
 
+    virtual void getSubBoundary(BoundaryObject bound, 
+				vector<BoundaryObject> &subBound) const = 0;
 
     /** \} */
 

AMDiS/src/Tetrahedron.cc

@@ -345,8 +345,6 @@ namespace AMDiS {
   {
     FUNCNAME("Tetrahedron::getHigherOrderDofs()");
 
-    TEST_EXIT(dofGeoIndex != NULL)("Not yet implemented!\n");
-
     switch (bound.subObj) {
     case VERTEX:
       return;
@@ -394,14 +392,20 @@ namespace AMDiS {
 	  if (baseDofPtr) {
 	    do {
 	      if (elDofIter.getCurrentPos() == 1 && 
-		  elDofIter.getCurrentElementPos() == bound.ithObj)
+		  elDofIter.getCurrentElementPos() == bound.ithObj) {
 		dofs.push_back(elDofIter.getBaseDof());	
+		if (dofGeoIndex != NULL)
+		  dofGeoIndex->push_back(EDGE);
+	      }
 	    } while (elDofIter.nextStrict());
 	  } else {
 	    do {
 	      if (elDofIter.getCurrentPos() == 1 && 
-		  elDofIter.getCurrentElementPos() == bound.ithObj)
+		  elDofIter.getCurrentElementPos() == bound.ithObj) {
 		dofs.push_back(elDofIter.getDofPtr());	
+		if (dofGeoIndex != NULL)
+		  dofGeoIndex->push_back(EDGE);
+	      }
 	    } while (elDofIter.next());
 	  }
 	}
@@ -432,21 +436,27 @@ namespace AMDiS {
 					 baseDofPtr, dofGeoIndex);
 	  }
 	} else { 
+
 	  ElementDofIterator elDofIter(feSpace, true);
 	  elDofIter.reset(this);
 	  if (baseDofPtr) {
 	    do {
-	      if (elDofIter.getCurrentPos() == 2 && 
+	      if (elDofIter.getCurrentPos() >= 1 && 
 		  elDofIter.getCurrentElementPos() == bound.ithObj)
 		dofs.push_back(elDofIter.getBaseDof());	
+		if (dofGeoIndex != NULL)
+		  dofGeoIndex->push_back(elDofIter.getPosIndex());
 	    } while (elDofIter.nextStrict());
 	  } else {
 	    do {
-	      if (elDofIter.getCurrentPos() == 2 && 
+	      if (elDofIter.getCurrentPos() >= 1 && 
 		  elDofIter.getCurrentElementPos() == bound.ithObj)
 		dofs.push_back(elDofIter.getDofPtr());	
+		if (dofGeoIndex != NULL)
+		  dofGeoIndex->push_back(elDofIter.getPosIndex());
 	    } while (elDofIter.next());      
 	  }
+
 	}
       }
       break;

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -65,7 +65,7 @@ namespace AMDiS {
 
     matAssembly();
 
-    removeDirichletRows(seqMat);
+    //    removeDirichletRows(seqMat);
     
     if (printMatInfo) {
       MatInfo matInfo;
@@ -143,6 +143,8 @@ namespace AMDiS {
 	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
 	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
 
+	std::set<DegreeOfFreedom> &dirichletRows = dofMat->getDirichletRows();
+
 	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
 	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
 	
@@ -151,6 +153,14 @@ namespace AMDiS {
 	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
 
 	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
+
+	  // For the case, this is a dirichlet row we have to check whether the 
+	  // rank is also owner of this row DOF.
+	  if (dirichletRows.count(*cursor)) {
+	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) ||
+		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor)))
+	      continue;
+	  }
   
 	  cols.clear();
 	  colsOther.clear();
@@ -236,7 +246,7 @@ namespace AMDiS {
 
     matAssembly();
 
-    removeDirichletRows(seqMat);
+    //    removeDirichletRows(seqMat);
 
     // === Create solver for the non primal (thus local) variables. ===
 
@@ -279,7 +289,7 @@ namespace AMDiS {
 
     vecRhsAssembly();
 
-    removeDirichletRows(vec);
+    //    removeDirichletRows(vec);
 
     // === For debugging allow to write the rhs vector to a file. ===
 
@@ -575,6 +585,15 @@ namespace AMDiS {
     }
 
     if (hasCoarseSpace()) {
+      MatZeroRows(getMatInteriorCoarse(0), 
+		  dRowsInterior.size(), &(dRowsInterior[0]), 0.0, PETSC_NULL, PETSC_NULL);
+      MatZeroRows(getMatInteriorCoarse(1), 
+		  dRowsInterior.size(), &(dRowsInterior[0]), 0.0, PETSC_NULL, PETSC_NULL);
+
+
+      MatZeroRows(getMatCoarseInterior(0), 
+		  dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0, PETSC_NULL, PETSC_NULL);
+
       Mat mpiMat = getMatCoarse();
       MatZeroRows(mpiMat, dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0,
 		  PETSC_NULL, PETSC_NULL);
@@ -745,7 +764,7 @@ namespace AMDiS {
 
 
   void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
-					     int nRowMat, int nColMat)
+					     int rowComp, int colComp)
   {
     FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");
 
@@ -770,7 +789,8 @@ namespace AMDiS {
 
     // Get periodic mapping object
     PeriodicMap &perMap = meshDistributor->getPeriodicMap();
-
+    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
+      
     const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
     const FiniteElemSpace *colFe = seqMat->getColFeSpace();
 
@@ -781,7 +801,7 @@ namespace AMDiS {
 	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
       // Global index of the current row DOF.
       MultiIndex rowMultiIndex;
-      if ((*interiorMap)[nRowMat].find(*cursor, rowMultiIndex) == false)
+      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
 	continue;
 
       int globalRowDof = rowMultiIndex.global;
@@ -789,11 +809,15 @@ namespace AMDiS {
       // Test if the current row DOF is a periodic DOF.
       bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
 
+      // Dirichlet rows can be set only be the owner ranks.
+      if (dirichletRows.count(*cursor) && !((*interiorMap)[rowComp].isRankDof(*cursor)))
+	continue;
+
       if (!periodicRow) {
 	// === Row DOF index is not periodic. ===
 
 	// Get PETSc's mat row index.
-	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
+	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
 
 	cols.clear();
 	values.clear();
@@ -803,14 +827,14 @@ namespace AMDiS {
 
 	  // Global index of the current column index.
 	  MultiIndex colMultiIndex;
-	  if ((*interiorMap)[nColMat].find(col(*icursor), colMultiIndex) == false)
+	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
 	    continue;
 
 	  int globalColDof = colMultiIndex.global;
 	  // Test if the current col dof is a periodic dof.
 	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
 	  // Get PETSc's mat col index.
-	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
+	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
 
 	  // Ignore all zero entries, expect it is a diagonal entry.
  	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
@@ -839,7 +863,7 @@ namespace AMDiS {
 	    vector<int> newCols;
 	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
 	    for (unsigned int i = 0; i < newCols.size(); i++) {
-	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
+	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
 	      values.push_back(scaledValue);	      
 	    }
 	  }
@@ -861,7 +885,7 @@ namespace AMDiS {
 	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
 	     icursor != icend; ++icursor) {
 	  // Global index of the current column index.
-	  int globalColDof = (*interiorMap)[nColMat][col(*icursor)].global;
+	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
 
 	  // Ignore all zero entries, expect it is a diagonal entry.
  	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
@@ -891,8 +915,8 @@ namespace AMDiS {
 	  // === Translate the matrix entries to PETSc's matrix.
 
 	  for (unsigned int i = 0; i < entry.size(); i++) {
-	    int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first);
-	    int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second);
+	    int rowIdx = interiorMap->getMatIndex(rowComp, entry[i].first);
+	    int colIdx = interiorMap->getMatIndex(colComp, entry[i].second);
 
 	    colsMap[rowIdx].push_back(colIdx);
 	    valsMap[rowIdx].push_back(scaledValue);
@@ -920,7 +944,7 @@ namespace AMDiS {
   void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
 					     Vec vecCoarse,
 					     DOFVector<double>* vec, 
-					     int nRowVec, 
+					     int rowComp, 
 					     bool rankOnly)
   {
     FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");
@@ -929,36 +953,43 @@ namespace AMDiS {
     PeriodicMap &perMap = meshDistributor->getPeriodicMap();
 
     ParallelDofMapping *rowCoarseSpace = 
-      (coarseSpaceMap.size() ? coarseSpaceMap[nRowVec] : NULL);
+      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);
+
+    map<DegreeOfFreedom, double> &dirichletValues = vec->getDirichletValues();
 
     // Traverse all used DOFs in the dof vector.
     DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
     for (dofIt.reset(); !dofIt.end(); ++dofIt) {
+      
+      DegreeOfFreedom dof = dofIt.getDOFIndex();
+      
+      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
+	continue;
 
-      if (rankOnly && !(*interiorMap)[nRowVec].isRankDof(dofIt.getDOFIndex()))
+      // Dirichlet rows can be set only be the owner ranks.
+      if (dirichletValues.count(dof) && !((*interiorMap)[rowComp].isRankDof(dof)))
 	continue;
 
-      if (isCoarseSpace(nRowVec, dofIt.getDOFIndex())) {
+      if (isCoarseSpace(rowComp, dof)) {
 	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");
 
-	int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
+	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
 	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
       } else {
-	if ((*interiorMap)[nRowVec].isSet(dofIt.getDOFIndex()) == false)
+	if ((*interiorMap)[rowComp].isSet(dof) == false)
 	  continue;
 
 	// Calculate global row index of the DOF.
-	DegreeOfFreedom globalRowDof = 
-	  (*interiorMap)[nRowVec][dofIt.getDOFIndex()].global;
+	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
 	
 	// Get PETSc's mat index of the row DOF.
 	int index = 0;
 	if (interiorMap->isMatIndexFromGlobal())
 	  index = 
-	    interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior;
+	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
 	else
 	  index =
-	    interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
+	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
 
 	if (perMap.isPeriodic(feSpace, globalRowDof)) {
 	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
@@ -968,7 +999,7 @@ namespace AMDiS {
 	  for (std::set<int>::iterator perIt = perAsc.begin(); 
 	       perIt != perAsc.end(); ++perIt) {
 	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
-	    int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
+	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
 
 	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
 	  }	  

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -94,19 +94,21 @@ namespace AMDiS {
     virtual void exitPreconditioner(PC pc);
 
     /// Takes a DOF matrix and sends the values to the global PETSc matrix.
-    void setDofMatrix(DOFMatrix* mat, int nRowMat = 0, int nColMat = 0);
+    void setDofMatrix(DOFMatrix* mat, int rowComp = 0, int colComp = 0);
 
     /// Takes a DOF vector and sends its values to a given PETSc vector.
     void setDofVector(Vec vecInterior,
 		      Vec vecCoarse,
 		      DOFVector<double>* vec, 
-		      int nRowVec, bool rankOnly = false);
+		      int rowCompxo, 
+		      bool rankOnly = false);
 
     inline void setDofVector(Vec vecInterior,
 			     DOFVector<double>* vec, 
-			     int nRowVec, bool rankOnly = false)
+			     int rowComp, 
+			     bool rankOnly = false)
     {
-      setDofVector(vecInterior, PETSC_NULL, vec, nRowVec, rankOnly);
+      setDofVector(vecInterior, PETSC_NULL, vec, rowComp, rankOnly);
     }
 
     inline void setDofVector(Vec vecInterior,