Blbu

AMDiS/src/parallel/MatrixNnzStructure.cc

@@ -191,13 +191,13 @@ namespace AMDiS {
 		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
 		if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false)
 		  continue;
-
+		
 		// Set of periodic row associations (is empty, if row DOF is not 
 		// periodic.
-		std::set<int> perColAsc;
+		std::set<int> perColAsc = perRowAsc;
 		if (perMap)
 		  perMap->fillAssociations(colFeSpace, colDofIndex.global, elObjDb, perColAsc);
-
+		
 		if (perColAsc.empty()) {
 		  if (colDofMap[colFeSpace].isRankDof(col(*icursor)))
 		    dnnz[localPetscRowIdx]++;
@@ -206,7 +206,7 @@ namespace AMDiS {
 		} else {
 		  vector<int> newCols;
 		  perMap->mapDof(colFeSpace, colDofIndex.global, perColAsc, newCols);
-
+		  
 		  for (int aa = 0; aa < newCols.size(); aa++) {
 		    int petscColIdx = colDofMap.getMatIndex(colComp, newCols[aa]);
 		    
@@ -221,6 +221,14 @@ namespace AMDiS {
 		  }
 		}
 	      }
+	      
+ 	      if (!perRowAsc.empty()) {
+		vector<int> newRows;
+		perMap->mapDof(rowFeSpace, rowIt->second.global, perRowAsc, newRows);
+		
+		dnnz[localPetscRowIdx] += 
+		  (newRows.size() - 1) * (onnz[localPetscRowIdx] + dnnz[localPetscRowIdx]);
+ 	      }
 	    }
 	  } else {
 	    // === The current row DOF is not a rank DOF, i.e., its values   ===

AMDiS/src/parallel/MeshDistributor.cc

@@ -1818,13 +1818,13 @@ namespace AMDiS {
     TEST_EXIT(levelData.getLevelNumber() == 1)
       ("Periodic DOF map does not support multi level domain decomposition!\n");
 
-    MPI::COMM_WORLD.Barrier();
+    //    MPI::COMM_WORLD.Barrier();   [TODO: CHANGE BECAUSE NOT ALL RANKS HAVE PERIODIC MAP!!!]
     double first = MPI::Wtime();
     
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       createPeriodicMap(feSpaces[i]);
 
-    MPI::COMM_WORLD.Barrier();
+    //    MPI::COMM_WORLD.Barrier();
     INFO(info, 8)("Creation of periodic mapping needed %.5f seconds\n", 
 		  MPI::Wtime() - first);
   }

AMDiS/src/parallel/ParallelDofMapping.h

@@ -184,13 +184,18 @@ namespace AMDiS {
       return static_cast<bool>(dofMap.count(dof));
     }
 
-    /// Checks if a given DOF is a rank owned DOF of the DOF mapping. The DOF must
-    /// a DOF of the mapping (this is not checked here), otherwise the result is
-    /// meaningsless.
+    /// Checks if a given DOF is a rank owned DOF of the DOF mapping. The DOF
+    /// must be a DOF of the mapping (this is not checked here), otherwise the
+    /// result is meaningsless.
     bool isRankDof(DegreeOfFreedom dof)
     {
       return !(static_cast<bool>(nonRankDofs.count(dof)));
     }
+
+    bool isRankGlobalDof(int dof)
+    {
+      return (dof >= rStartDofs && dof < rStartDofs + nRankDofs);
+    }
     
     /// Returns number of DOFs in the mapping.
     unsigned int size()

AMDiS/src/parallel/PeriodicMap.cc

@@ -72,6 +72,37 @@ namespace AMDiS {
   }
 
 
+  void PeriodicMap::mapDof(const FiniteElemSpace* rowFeSpace,
+			   const FiniteElemSpace* colFeSpace,
+			   pair<int, int> globalDofIndex, 
+			   const std::set<int>& perAsc, 
+			   vector<pair<int, int> >& mappedDofs)
+  {
+    mappedDofs.clear();
+    mappedDofs.push_back(globalDofIndex);
+
+    for (std::set<int>::iterator it = perAsc.begin(); 
+	 it != perAsc.end(); ++it) {
+      int nDofs = static_cast<int>(mappedDofs.size());
+      for (int i = 0; i < nDofs; i++) {
+	int perRowDof = 0;	
+	if (isPeriodic(rowFeSpace, *it, mappedDofs[i].first))
+	  perRowDof = map(rowFeSpace, *it, mappedDofs[i].first);
+	else
+	  perRowDof = mappedDofs[i].first;
+	
+	int perColDof;
+	if (isPeriodic(colFeSpace, *it, mappedDofs[i].second))
+	  perColDof = map(colFeSpace, *it, mappedDofs[i].second);
+	else
+	  perColDof = mappedDofs[i].second;	      	      
+	
+	mappedDofs.push_back(make_pair(perRowDof, perColDof));
+      }
+    }    
+  }
+
+
   void PeriodicMap::serialize(ostream &out, 
 			      vector<const FiniteElemSpace*> feSpaces)
   {

AMDiS/src/parallel/PeriodicMap.h

@@ -159,6 +159,21 @@ namespace AMDiS {
 		const std::set<int>& perAsc, 
 		vector<int>& mappedDofs);
 
+    /** \brief
+     * Maps a given DOF index pair for all given periodic DOF associations.
+     *
+     * \param[in]   rowFeSpace      feSpace of the DOFs on the first component.
+     * \param[in]   colFeSpace      feSpace of the DOFs on the second component.
+     * \param[in]   globalDofIndex  pair of global index of a DOF.
+     * \param[in]   perAsc          set of periodic associations.
+     * \param[out]  mappedDofs      set of pairs of global DOF indices.
+     */    
+    void mapDof(const FiniteElemSpace* rowFeSpace,
+		const FiniteElemSpace* colFeSpace,
+		pair<int, int> globalDofIndex, 
+		const std::set<int>& perAsc, 
+		vector<pair<int, int> >& mappedDofs);
+
 
     /// Returns true, if the DOF (global index) is a periodic DOF.
     inline bool isPeriodic(const FiniteElemSpace *feSpace, int globalDofIndex)

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -805,32 +805,8 @@ namespace AMDiS {
 	  // === associations of the row and column indices.             ===
 
 	  vector<pair<int, int> > entry;
-	  
-	  // First, add the original entry.
-	  entry.push_back(make_pair(globalRowDof, globalColDof));
-
-	  // Then, traverse the periodic associations of the row and column
-	  // indices and create the corresponding entries.
-	  for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
-	    int nEntry = static_cast<int>(entry.size());
-	    for (int i = 0; i < nEntry; i++) {
-	      int perRowDof = 0;
-
-	      if (perMap.isPeriodic(rowFe, *it, entry[i].first))
-		perRowDof = perMap.map(rowFe, *it, entry[i].first);
-	      else
-		perRowDof = entry[i].first;
-
-	      int perColDof;
-	      if (perMap.isPeriodic(colFe, *it, entry[i].second))
-		perColDof = perMap.map(colFe, *it, entry[i].second);
-	      else
-		perColDof = entry[i].second;	      	      
-	      
-	      entry.push_back(make_pair(perRowDof, perColDof));
-	    }
-	  }
-
+	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
+			perAsc, entry);
 
 	  // === Translate the matrix entries to PETSc's matrix.