Fixed small problem with periodic boundary conditions in 2D, thanks to Simon for bug reporting.

AMDiS/src/io/FileWriter.cc

@@ -159,14 +159,6 @@ namespace AMDiS {
 //       vtkVectorWriter.setWriteAs3dVector(writeAs3dVector);
 //       vtkVectorWriter.writeFile(fn + paraviewFileExt);
 
-#if HAVE_PARALLEL_DOMAIN_AMDIS
-      ERROR_EXIT("Simon must make some work in this, does not compile in parallel!\n");
-//       if (MPI::COMM_WORLD.Get_rank() == 0)
-// 	vtkVectorWriter.writeParallelFile(paraFilename + paraviewParallelFileExt,
-// 					  MPI::COMM_WORLD.Get_size(),
-// 					  filename, postfix);
-#endif
-
       MSG("ParaView file written to %s\n", (fn + paraviewFileExt).c_str());
     }
 

AMDiS/src/parallel/ElementObjectDatabase.cc

@@ -140,6 +140,12 @@ namespace AMDiS {
     if (periodicVertices.size() == 0)
       return;
  
+    // === Calculate smallest periodic boundary ID in mesh. ===
+
+    smallestPeriodicBcType = 0;
+    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
+	 it != mesh->getPeriodicAssociations().end(); ++it)
+      smallestPeriodicBcType = std::min(smallestPeriodicBcType, it->first);
 
     // === Get all vertex DOFs that have multiple periodic associations. ===
     

AMDiS/src/parallel/ElementObjectDatabase.h

@@ -450,6 +450,14 @@ namespace AMDiS {
       return (periodicVertices.size() != 0);
     }
 
+    /// Returns true if the given boundary type is larger or equal to the smallest
+    /// periodic boundary ID in mesh. See \ref smallestPeriodicBcType for more
+    /// information.
+    bool isValidPeriodicType(BoundaryType t) const
+    {
+      return (t >= smallestPeriodicBcType);
+    }
+
     /// Write the element database to disk.
     void serialize(ostream &out);
     
@@ -574,15 +582,21 @@ namespace AMDiS {
 
     map<pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
 
-    // The following three data structures store periodic DOFs, edges and faces.
+    /// The following three data structures store periodic DOFs, edges and faces.
     PerBoundMap<DegreeOfFreedom>::type periodicVertices;
     PerBoundMap<DofEdge>::type periodicEdges;
     PerBoundMap<DofFace>::type periodicFaces;
 
-    // Stores to each vertex all its periodic associations.
+    /// Defines the smallest boudary ID for periodic boundary conditions. This is
+    /// required to distinguish between "real" periodic boundaries and periodic
+    /// boundary IDs that are set by the parallel algorithm for indirectly 
+    /// connected boundaries.
+    BoundaryType smallestPeriodicBcType;
+
+    /// Stores to each vertex all its periodic associations.
     map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
 
-    // Stores to each edge all its periodic associations.
+    /// Stores to each edge all its periodic associations.
     map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
 
     map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode;

AMDiS/src/parallel/MeshDistributor.cc

@@ -132,6 +132,8 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::initParallelization()");
 
+    MSG("RUN INIT\n");
+
     if (initialized)
       return;
 
@@ -2067,13 +2069,6 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createPeriodicMap()");
 
-
-//     vector<ElementObjectData> objs = elObjDb.getElements(0);
-//     MSG("SIZE OF 0-DOF is: %d\n", objs.size());
-
-//     for (int i = 0; i < objs.size(); i++)
-//       MSG("0-DOF IN EL %d - %d\n", objs[i].elIndex, objs[i].ithObject);
-
     StdMpi<vector<int> > stdMpi(mpiComm, false);
 
     // === Each rank traverse its periodic boundaries and sends the DOF      ===
@@ -2112,17 +2107,10 @@ namespace AMDiS {
 	    DegreeOfFreedom globalDof1 = 
 	      dofFeData[feSpace].mapDofToGlobal[*(dofs1[j])];
 
-	    if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0)) {
-// 	      if (globalDof0 == 0) {
-// 		MSG("A-ADD DOF 0 at BOUND: %d %d %d  <-> %d %d %d \n",
-// 		    bound.rankObj.elIndex, bound.rankObj.subObj, bound.rankObj.ithObj,
-// 		    bound.neighObj.elIndex, bound.neighObj.subObj, bound.neighObj.ithObj);
-// 	      }
-
+	    if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0))
 	      periodicMap.add(feSpace, type, globalDof0, globalDof1);
 	  }
 	}
-	}
 
       } else {
 	// Here we have a periodic boundary between two ranks.
@@ -2132,9 +2120,6 @@ namespace AMDiS {
 	for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
 	     boundIt != it->second.end(); ++boundIt) {
 
-// 	  MSG("  NOW ADD DOFS ON PER BOUND EL %d %d with rank %d\n",
-// 	      boundIt->rankObj.elIndex, boundIt->rankObj.subObj, it->first);
-
 	  int nDofs = dofs.size();
 	  boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs);
 
@@ -2177,16 +2162,10 @@ namespace AMDiS {
 
 	// Check if this global DOF with the corresponding boundary type was
 	// not added before by another periodic boundary from other rank.
-	if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDofIndex)) {
-// 	  MSG("B-ADD DOF %d -> %d at BOUND from rank %d\n", 
-// 	      globalDofIndex, 
-// 	      mapGlobalDofIndex,
-// 	      it->first);
-
+	if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDofIndex))
 	  periodicMap.add(feSpace, type, globalDofIndex, mapGlobalDofIndex);
       }
     }
-    }
 
 
     StdMpi<PeriodicDofMap> stdMpi2(mpiComm);
@@ -2201,7 +2180,9 @@ namespace AMDiS {
 
       for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
 	   boundIt != it->second.end(); ++boundIt) {
-	if (boundIt->type >= -3)
+	// If the boundary type is a normal periodic BC, continue. We just
+	// operate on indirectly connected periodic boundaries.
+	if (elObjDb.isValidPeriodicType(boundIt->type))
 	  continue;
 	
 	DofContainer dofs;
@@ -2217,17 +2198,11 @@ namespace AMDiS {
 
 	  for (std::set<BoundaryType>::iterator perAscIt = assoc.begin();
 	       perAscIt != assoc.end(); ++perAscIt)
-	    if (*perAscIt >= -3) {
-// 	      MSG(" add map DOF %d -> %d (%d)\n", 
-// 		  globalDof,
-// 		  periodicMap.map(feSpace, *perAscIt, globalDof),
-// 		  *perAscIt);
-
+	    if (elObjDb.isValidPeriodicType(*perAscIt))
 	      perObjMap[*perAscIt][globalDof] = 
 		periodicMap.map(feSpace, *perAscIt, globalDof);	    
 	}
       }
-      }
 
       stdMpi2.send(it->first, perObjMap);
       stdMpi2.recv(it->first);
@@ -2240,7 +2215,7 @@ namespace AMDiS {
     //       DOFs. I'm note sure what is correct here.
     for (map<int, PeriodicDofMap>::iterator it = stdMpi2.getRecvData().begin();
 	 it != stdMpi2.getRecvData().end(); ++it)
-      periodicMap.add(feSpace, it->second);
+      periodicMap.add(feSpace, it->second, false);
   }
 
 

AMDiS/src/parallel/MeshDistributor.h

@@ -443,12 +443,14 @@ namespace AMDiS {
     void getAllBoundaryDofs(const FiniteElemSpace *feSpace,
 			    DofContainer& dofs);
 
+    const ElementObjectDatabase& getElementObjectDb() 
+    {
+      return elObjDb;
+    }
 
-  public:
     /// Adds a stationary problem to the global mesh distributor objects.
     static void addProblemStatGlobal(ProblemStatSeq *probStat);
     
-    
   protected:
     void addProblemStat(ProblemStatSeq *probStat);
 

AMDiS/src/parallel/PeriodicMap.cc

@@ -16,14 +16,15 @@
 namespace AMDiS {
 
   void PeriodicMap::add(const FiniteElemSpace *feSpace,
-			PeriodicDofMap &newMap)
+			PeriodicDofMap &newMap,
+			bool a)
   {
     FUNCNAME("PeriodicMap::add()");
     
     for (PeriodicDofMap::iterator it = newMap.begin(); it != newMap.end(); ++it)
       for (DofMap::iterator dofIt =it->second.begin();
 	   dofIt != it->second.end(); ++dofIt)
-	add(feSpace, it->first, dofIt->second, dofIt->first);
+	add(feSpace, it->first, dofIt->second, dofIt->first, a);
   }
 
 

AMDiS/src/parallel/PeriodicMap.h

@@ -101,7 +101,8 @@ namespace AMDiS {
     /// this DOFs that maps to a different DOF index than the given one.
     inline void add(const FiniteElemSpace *feSpace,
 		    BoundaryType type,
-		    DegreeOfFreedom dof0, DegreeOfFreedom dof1)
+		    DegreeOfFreedom dof0, DegreeOfFreedom dof1,
+		    bool a = true)
     {
       FUNCNAME("PeriodicMap::map()");
 
@@ -110,12 +111,13 @@ namespace AMDiS {
 	("Should not happen!\n");
 
       periodicDofMap[feSpace][type][dof0] = dof1;
+      if (a)
 	periodicDofAssociations[feSpace][dof0].insert(type);
     }
 
     
     /// Adds a whole periodic mapping to the current one.
-    void add(const FiniteElemSpace *feSpace, PeriodicDofMap &newMap);
+    void add(const FiniteElemSpace *feSpace, PeriodicDofMap &newMap, bool a = true);
 
 
     /// For a given global DOF index, this function returns the set of periodic

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -307,7 +307,7 @@ namespace AMDiS {
 	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
 	    for (std::set<int>::iterator it = perColAsc.begin(); 
 		 it != perColAsc.end(); ++it)
-	      if (*it >= -3)
+	      if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
 		perAsc.insert(*it);
     
 	    // Scale value to the number of periodic associations of the column index.
@@ -377,14 +377,14 @@ namespace AMDiS {
 	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
 	    for (std::set<int>::iterator it = perColAsc.begin(); 
 		 it != perColAsc.end(); ++it)
-	      if (*it >= -3)
+	      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 (*it >= -3)
+	    if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
 	      perAsc.insert(*it);
 
 	  // Scale the value with respect to the number of periodic associations.