Better matrix nnz computation when using periodic boundary conditions.

AMDiS/src/parallel/MatrixNnzStructure.cc

@@ -37,6 +37,7 @@ namespace AMDiS {
 				  ParallelDofMapping &rowDofMap,
 				  ParallelDofMapping &colDofMap,
 				  PeriodicMap *perMap,
+				  const ElementObjectDatabase& elObjDb,
 				  bool localMatrix)
   {
     FUNCNAME("MatrixNnzStructure::create()");
@@ -147,8 +148,12 @@ namespace AMDiS {
 	      petscRowIdx = rowDofMap.getMatIndex(rowComp, *cursor);
 	  }
 
-	  bool periodicRow = 
-	    (perMap && perMap->isPeriodic(rowFeSpace, rowIt->second.global));
+
+	  // Set of periodic row associations (is empty, if row DOF is not 
+	  // periodic.
+	  std::set<int> perRowAsc;
+	  if (perMap)
+	    perMap->fillAssociations(rowFeSpace, rowIt->second.global, elObjDb, perRowAsc);
 
 
 	  if (localMatrix || rowDofMap[rowFeSpace].isRankDof(*cursor)) {
@@ -186,24 +191,35 @@ namespace AMDiS {
 		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
 		if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false)
 		  continue;
-		
-		int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? 
-				   colDofMap.getMatIndex(colComp, colDofIndex.global) :
-				   colDofMap.getMatIndex(colComp, col(*icursor)));
-	
-		// The row DOF is a rank DOF, if also the column is a rank DOF, 
-		// increment the diagNnz values for this row, 
-		// otherwise the offdiagNnz value.
-		if (petscColIdx >= rankStartColIndex && 
-		    petscColIdx < rankStartColIndex + nRankCols)
-		  dnnz[localPetscRowIdx]++;
-		else
-		  onnz[localPetscRowIdx]++;
-	      }
 
-	      if (periodicRow) {
-		dnnz[localPetscRowIdx] *= 4;
-		onnz[localPetscRowIdx] *= 4;
+		// Set of periodic row associations (is empty, if row DOF is not 
+		// periodic.
+		std::set<int> perColAsc;
+		if (perMap)
+		  perMap->fillAssociations(colFeSpace, colDofIndex.global, elObjDb, perColAsc);
+
+		if (perColAsc.empty()) {
+		  if (colDofMap[colFeSpace].isRankDof(col(*icursor)))
+		    dnnz[localPetscRowIdx]++;
+		  else
+		    onnz[localPetscRowIdx]++;
+		} 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]);
+		    
+		    // The row DOF is a rank DOF, if also the column is a rank DOF, 
+		    // increment the diagNnz values for this row, 
+		    // otherwise the offdiagNnz value.
+		    if (petscColIdx >= rankStartColIndex && 
+			petscColIdx < rankStartColIndex + nRankCols)
+		      dnnz[localPetscRowIdx]++;
+		    else
+		      onnz[localPetscRowIdx]++;
+		  }
+		}
 	      }
 	    }
 	  } else {

AMDiS/src/parallel/MatrixNnzStructure.h

@@ -44,16 +44,18 @@ namespace AMDiS {
 		MPI::Intracomm mpiComm,
 		ParallelDofMapping &rowDofMap,
 		ParallelDofMapping &colDofMap,
-		PeriodicMap *perMap = NULL,
+		PeriodicMap *perMap,
+		const ElementObjectDatabase& elObjDb,
 		bool localMatrix = false);
 
     void create(Matrix<DOFMatrix*> *mat,
 		MPI::Intracomm mpiComm,
 		ParallelDofMapping &dofMap,
-		PeriodicMap *perMap = NULL,
+		PeriodicMap *perMap,
+		const ElementObjectDatabase& elObjDb,
 		bool localMatrix = false)
     {
-      create(mat, mpiComm, dofMap, dofMap, perMap, localMatrix);
+      create(mat, mpiComm, dofMap, dofMap, perMap, elObjDb, localMatrix);
     }
 
   protected:

AMDiS/src/parallel/PeriodicMap.cc

@@ -13,6 +13,7 @@
 #include <map>
 #include <fstream>
 #include "Global.h"
+#include "parallel/ElementObjectDatabase.h"
 #include "parallel/PeriodicMap.h"
 #include "parallel/ParallelTypes.h"
 
@@ -33,6 +34,44 @@ namespace AMDiS {
   }
 
 
+  void PeriodicMap::fillAssociations(const FiniteElemSpace* feSpace,
+				     int globalDofIndex,
+				     const ElementObjectDatabase &elObjDb,
+				     std::set<int>& perAsc)
+  {
+    if (isPeriodic(feSpace, globalDofIndex) == false)
+      return;
+
+    std::set<int>& tmpAsc = getAssociations(feSpace, globalDofIndex);
+    for (std::set<int>::iterator it = tmpAsc.begin(); it != tmpAsc.end(); ++it)
+      if (elObjDb.isValidPeriodicType(*it))
+	perAsc.insert(*it);
+  }
+
+
+  void PeriodicMap::mapDof(const FiniteElemSpace* feSpace,
+			   int globalDofIndex, 
+			   const std::set<int>& perAsc, 
+			   vector<int>& mappedDofs)
+  {
+    mappedDofs.clear();
+    mappedDofs.push_back(globalDofIndex);
+    
+    for (std::set<int>::const_iterator it = perAsc.begin();
+	 it != perAsc.end(); ++it) {
+      int nDofs = static_cast<int>(mappedDofs.size());
+      
+      for (int i = 0; i < nDofs; i++) {
+	TEST_EXIT_DBG(isPeriodic(feSpace, *it, mappedDofs[i]))
+	  ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
+	   *it, mappedDofs[i]);
+	
+	mappedDofs.push_back(map(feSpace, *it, mappedDofs[i]));
+      }
+    }
+  }
+
+
   void PeriodicMap::serialize(ostream &out, 
 			      vector<const FiniteElemSpace*> feSpaces)
   {

AMDiS/src/parallel/PeriodicMap.h

@@ -26,6 +26,7 @@
 #include <map>
 #include <set>
 #include <vector>
+#include "AMDiS_fwd.h"
 #include "parallel/ParallelTypes.h"
 #include "Boundary.h"
 #include "Serializer.h"
@@ -126,6 +127,39 @@ namespace AMDiS {
     }
 
 
+    /** \brief
+     * Given a global DOF index in some given finite element space, this
+     * function creates all valid periodic associations of this DOF. If the
+     * DOF is not periodic, this function does not fail but has no effect on
+     * the output data.
+     *
+     * \param[in]  feSpace         feSpace on which the function should work on.
+     * \param[in]  globalDofIndex  global index of a DOF.
+     * \param[in]  elObjDb         Element object database that is used to check
+     *                             if a given periodic index is valid.
+     * \param[out] perAsc          set of periodic associations, is not deleted
+     *                             by this function
+     */
+    void fillAssociations(const FiniteElemSpace* feSpace,
+			  int globalDofIndex,
+			  const ElementObjectDatabase &elObjDb,
+			  std::set<int>& perAsc);
+
+    
+    /** \brief
+     * Maps a given DOF index for all given periodic DOF associations.
+     *
+     * \param[in]   feSpace         feSpace on which the function should work on.
+     * \param[in]   globalDofIndex  global index of a DOF.
+     * \param[in]   perAsc          set of periodic associations.
+     * \param[out]  mappedDofs      set of global DOF indices.
+     */
+    void mapDof(const FiniteElemSpace* feSpace,
+		int globalDofIndex, 
+		const std::set<int>& perAsc, 
+		vector<int>& mappedDofs);
+
+
     /// Returns true, if the DOF (global index) is a periodic DOF.
     inline bool isPeriodic(const FiniteElemSpace *feSpace, int globalDofIndex)
     {

AMDiS/src/parallel/PetscSolver.h

@@ -113,7 +113,6 @@ namespace AMDiS {
 
     void setKspPrefix(std::string s)
     {
-      MSG("======>>>>>> SET TO %s\n", s.c_str());
       kspPrefix = s;
     }
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -37,7 +37,8 @@ namespace AMDiS {
 
     if (checkMeshChange(mat))
       nnzInterior.create(mat, mpiCommGlobal, *interiorMap, 
-			 &(meshDistributor->getPeriodicMap()));
+			 &(meshDistributor->getPeriodicMap()),
+			 meshDistributor->getElementObjectDb());
 
     // === Create PETSc vector (solution and a temporary vector). ===
 
@@ -109,7 +110,6 @@ namespace AMDiS {
     KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); 
     KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
     KSPSetType(kspInterior, KSPBCGS);
-    MSG("SET PREFIX TOOO: %s\n", kspPrefix.c_str());
     KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
     KSPSetFromOptions(kspInterior);
 
@@ -148,12 +148,14 @@ namespace AMDiS {
 
     bool localMatrix = (subdomainLevel == 0);
     if (checkMeshChange(mat, localMatrix)) {
-      nnzInterior.create(mat, mpiCommGlobal, *interiorMap, NULL, localMatrix);
+      nnzInterior.create(mat, mpiCommGlobal, *interiorMap, NULL, 
+			 meshDistributor->getElementObjectDb(),
+			 localMatrix);
 
       if (coarseSpaceMap) {
-	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap);
-	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap);
-	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap);
+	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
+	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap, NULL, meshDistributor->getElementObjectDb());
+	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
       }
     }
 
@@ -743,12 +745,9 @@ namespace AMDiS {
 
 	    // Create set of all periodic associations of the column index.
 	    std::set<int> perAsc;
-	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
-	    for (std::set<int>::iterator it = perColAsc.begin(); 
-		 it != perColAsc.end(); ++it)
-	      if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
-		perAsc.insert(*it);
-    
+	    perMap.fillAssociations(colFe, globalColDof, 
+				    meshDistributor->getElementObjectDb(), perAsc);
+
 	    // Scale value to the number of periodic associations of the column index.
 	    double scaledValue = 
 	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
@@ -758,24 +757,7 @@ namespace AMDiS {
 	    // === associations of the column DOF index.                       ===
 
 	    vector<int> newCols;
-
-	    // First, add the original matrix index.
-	    newCols.push_back(globalColDof);
-
-	    // And add all periodic matrix indices.
-	    for (std::set<int>::iterator it = perAsc.begin(); 
-		 it != perAsc.end(); ++it) {
-	      int nCols = static_cast<int>(newCols.size());
-
-	      for (int i = 0; i < nCols; i++) {
- 		TEST_EXIT_DBG(perMap.isPeriodic(colFe, *it, newCols[i]))
- 		  ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
-		   *it, newCols[i]);
-
-		newCols.push_back(perMap.map(colFe, *it, newCols[i]));
-	      }
-	    }
-
+	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
 	    for (unsigned int i = 0; i < newCols.size(); i++) {
 	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
 	      values.push_back(scaledValue);	      
@@ -809,20 +791,10 @@ namespace AMDiS {
 	  // === indices to the set perAsc.                                    ===
 
 	  std::set<int> perAsc;
-
-	  if (perMap.isPeriodic(colFe, globalColDof)) {
-	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
-	    for (std::set<int>::iterator it = perColAsc.begin(); 
-		 it != perColAsc.end(); ++it)
-	      if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
-		perAsc.insert(*it);
-	  }
-
-	  std::set<int>& perRowAsc = perMap.getAssociations(rowFe, globalRowDof);
-	  for (std::set<int>::iterator it = perRowAsc.begin(); 
-	       it != perRowAsc.end(); ++it)
-	    if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
-	      perAsc.insert(*it);
+	  perMap.fillAssociations(colFe, globalColDof, 
+				  meshDistributor->getElementObjectDb(), perAsc);
+	  perMap.fillAssociations(rowFe, globalRowDof, 
+				  meshDistributor->getElementObjectDb(), perAsc);
 
 	  // Scale the value with respect to the number of periodic associations.
 	  double scaledValue =