Added some parallel debuging test for periodic boundary conditions.

AMDiS/src/io/MacroInfo.cc

@@ -525,19 +525,19 @@ namespace AMDiS {
 	      int j1 = mel_vertex[i][(j + 1) % 3];
 	      int j2 = mel_vertex[i][(j + 2) % 3];
 	      
-	      bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
-	      bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
+	      bound[j1] = std::max(bound[j1], (*melIt)->getBoundary(j));
+	      bound[j2] = std::max(bound[j2], (*melIt)->getBoundary(j));
 	    } else if ((*melIt)->getBoundary(j) <= NEUMANN) {
 	      int j1 = mel_vertex[i][(j + 1) % 3];
 	      int j2 = mel_vertex[i][(j + 2) % 3];
 	      
 	      if (bound[j1] != INTERIOR)
-		bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
+		bound[j1] = std::max(bound[j1], (*melIt)->getBoundary(j));
 	      else
 		bound[j1] = (*melIt)->getBoundary(j);
 	      
 	      if (bound[j2] != INTERIOR)
-		bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
+		bound[j2] = std::max(bound[j2], (*melIt)->getBoundary(j));
 	      else
 		bound[j2] = (*melIt)->getBoundary(j);
 	    }

AMDiS/src/io/MacroReader.cc

@@ -495,23 +495,19 @@ namespace AMDiS {
 	      periodic[it->elementSide] = true;
 	}
 
-	for (k = 0; k < mesh->getGeo(EDGE); k++) {      
-	  /*********************************************************************/
-	  /* check for not counted edges                                       */
-	  /*********************************************************************/
+	for (k = 0; k < mesh->getGeo(EDGE); k++) {
+	  // === Check for not counted edges. ===
 	  n_neigh = 1;
 
 	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
 	    mesh->incrementNumberOfEdges(1);
-	    max_n_neigh = max(max_n_neigh, n_neigh);
+	    max_n_neigh = std::max(max_n_neigh, n_neigh);
 	  }
 	}
       
 	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
 	  neigh = (*(mel + i))->getNeighbour(k);
-	  /*********************************************************************/
-	  /* face is counted and dof is added by the element with bigger index */
-	  /*********************************************************************/
+	  // === Face is counted and dof is added by the element with bigger index. ===
 	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))  
 	    continue;
 	

AMDiS/src/io/PovrayWriter.cc

@@ -59,8 +59,8 @@ namespace AMDiS {
       min_value = max_value = (*values)[0];
 
     for (int index = 1; index < values->getSize(); index++) {
-      min_value = min(min_value, (*values)[index]);
-      max_value = max(max_value, (*values)[index]);
+      min_value = std::min(min_value, (*values)[index]);
+      max_value = std::max(max_value, (*values)[index]);
     }
 
     /* map DOFs to values */
@@ -126,13 +126,13 @@ namespace AMDiS {
     for (it.reset(); !it.end(); ++it) {
       // initialize mit and max values with coordinates of first vertex
       std::list<VertexInfo>::iterator it2 = it->begin();
-      bBox->minx = min(bBox->minx, it2->coords[0]);
-      bBox->maxx = max(bBox->maxx, it2->coords[0]);
-      bBox->miny = min(bBox->miny, it2->coords[1]);
-      bBox->maxy = max(bBox->maxy, it2->coords[1]);
+      bBox->minx = std::min(bBox->minx, it2->coords[0]);
+      bBox->maxx = std::max(bBox->maxx, it2->coords[0]);
+      bBox->miny = std::min(bBox->miny, it2->coords[1]);
+      bBox->maxy = std::max(bBox->maxy, it2->coords[1]);
       if (dow == 3) {
-	bBox->minz = min(bBox->minz, it2->coords[2]);
-	bBox->maxz = max(bBox->maxz, it2->coords[2]); 
+	bBox->minz = std::min(bBox->minz, it2->coords[2]);
+	bBox->maxz = std::max(bBox->maxz, it2->coords[2]); 
       }
     }
 
@@ -149,8 +149,8 @@ namespace AMDiS {
   // determines the color for a given value/weight
   string getColorString(double &value) 
   {
-    value = max(0.0, value);
-    value = min(value, 1.0);
+    value = std::max(0.0, value);
+    value = std::min(value, 1.0);
 
     // rot: 1,0,0
     // blau: 0,0,1
@@ -324,8 +324,8 @@ namespace AMDiS {
       for (std::list<VertexInfo>::iterator it2 = it->begin(); it2 != it->end(); ++it2) {
 	// test: use y-coordinate to compute color
 	double redValue = it2->coords[1];
-	redValue = max(0., redValue);
-	redValue = min(redValue, 1.);
+	redValue = std::max(0., redValue);
+	redValue = std::min(redValue, 1.);
 		    
 	out << "\n\t\ttexture{ pigment{ rgb"<< getColorString(redValue) <<" } }" << ",";
       }     

AMDiS/src/parallel/ElementObjectData.cc

@@ -44,6 +44,8 @@ namespace AMDiS {
 	if (mesh->isPeriodicAssociation(elInfo->getBoundary(EDGE, i))) {
 	  // The current element's i-th edge is periodic.
 	  Element *neigh = elInfo->getNeighbour(i);
+
+	  TEST_EXIT_DBG(neigh)("Should not happen!\n");
 	  
 	  DofEdge edge0 = el->getEdge(i);
 	  DofEdge edge1 = neigh->getEdge(elInfo->getOppVertex(i));
@@ -72,6 +74,8 @@ namespace AMDiS {
 	if (mesh->isPeriodicAssociation(elInfo->getBoundary(FACE, i))) {
 	  // The current element's i-th face is periodic.
 	  Element *neigh = elInfo->getNeighbour(i);
+
+	  TEST_EXIT_DBG(neigh)("Should not happen!\n");
 	  
 	  DofFace face0 = el->getFace(i);
 	  DofFace face1 = neigh->getFace(elInfo->getOppVertex(i));
@@ -132,21 +136,8 @@ namespace AMDiS {
     
     if (periodicVertices.size() == 0)
       return;
-   
-
-    // === Search for an unsed boundary index. ===
+ 
 
-    BoundaryType newPeriodicBoundaryType = 0;
-    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
-	 it != mesh->getPeriodicAssociations().end(); ++it)
-      newPeriodicBoundaryType = min(newPeriodicBoundaryType, it->first);
-    
-    TEST_EXIT_DBG(newPeriodicBoundaryType < 0)("Should not happen!\n");
-    newPeriodicBoundaryType--;
-    mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
-      new VertexVector(feSpace->getAdmin(), "");
-    
-    
     // === Get all vertex DOFs that have multiple periodic associations. ===
     
     // We group all vertices together, that have either two or three periodic
@@ -223,10 +214,12 @@ namespace AMDiS {
 	  ("Should not happen!\n");
 	TEST_EXIT_DBG(periodicVertices.count(make_pair(dof3, dof0)) == 0)
 	  ("Should not happen!\n");
-	
-	periodicVertices[make_pair(dof0, dof3)] = newPeriodicBoundaryType;
-	periodicVertices[make_pair(dof3, dof0)] = newPeriodicBoundaryType;
-	
+
+ 	periodicVertices[make_pair(dof0, dof3)] = 
+ 	  provideConnectedPeriodicBoundary(type0, type1);
+ 	periodicVertices[make_pair(dof3, dof0)] = 
+ 	  provideConnectedPeriodicBoundary(type0, type1);
+
 	for (unsigned int j = i + 1; j < multPeriodicDof2.size(); j++)
 	  if (multPeriodicDof2[j] == dof3)
 	    multPeriodicDof2[j] = -1;
@@ -235,6 +228,7 @@ namespace AMDiS {
 
     if (multPeriodicDof3.size() > 0) {
       int nMultPeriodicDofs = multPeriodicDof3.size();
+
       for (int i = 0; i < nMultPeriodicDofs; i++) {
 	for (int j = i + 1; j < nMultPeriodicDofs; j++) {
 	  pair<DegreeOfFreedom, DegreeOfFreedom> perDofs0 = 
@@ -243,20 +237,15 @@ namespace AMDiS {
 	    make_pair(multPeriodicDof3[j], multPeriodicDof3[i]);
 	  
 	  if (periodicVertices.count(perDofs0) == 0) {
-	    TEST_EXIT_DBG(periodicVertices.count(perDofs1) == 0)
-	      ("Should not happen!\n");
-	    
-	    periodicVertices[perDofs0] = newPeriodicBoundaryType;
-	    periodicVertices[perDofs1] = newPeriodicBoundaryType;
-	    newPeriodicBoundaryType--;
-	    mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
-	      new VertexVector(feSpace->getAdmin(), "");
+	    BoundaryType b = getNewBoundaryType();
+	    periodicVertices[perDofs0] = b;
+	    periodicVertices[perDofs1] = b;
 	  }
 	}
       }
     }
-    
-    
+  
+   
     // === Get all edges that have multiple periodic associations (3D only!). ===
     
     for (map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
@@ -265,16 +254,24 @@ namespace AMDiS {
 	TEST_EXIT_DBG(mesh->getDim() == 3)("Should not happen!\n");
 	TEST_EXIT_DBG(it->second.size() == 2)("Should not happen!\n");
 	
-	pair<DofEdge, DofEdge> perEdge0 =
-	  make_pair(*(it->second.begin()), *(++(it->second.begin())));
-	pair<DofEdge, DofEdge> perEdge1 =
-	  make_pair(perEdge0.second, perEdge0.first);
+	DofEdge edge0 = it->first;
+	DofEdge edge1 = *(it->second.begin());
+	DofEdge edge2 = *(++(it->second.begin()));
 	
-	periodicEdges[perEdge0] = newPeriodicBoundaryType;
-	periodicEdges[perEdge1] = newPeriodicBoundaryType;
-	newPeriodicBoundaryType--;
-	mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
-	  new VertexVector(feSpace->getAdmin(), "");      
+	pair<DofEdge, DofEdge> perEdge0 = make_pair(edge1, edge2);
+	pair<DofEdge, DofEdge> perEdge1 = make_pair(edge2, edge1);
+
+	TEST_EXIT_DBG(periodicEdges.count(make_pair(edge0, edge1)) == 1)
+	  ("Should not happen!\n");
+	TEST_EXIT_DBG(periodicEdges.count(make_pair(edge1, edge0)) == 1)
+	  ("Should not happen!\n");
+
+	BoundaryType type0 = periodicEdges[make_pair(edge0, edge1)];
+	BoundaryType type1 = periodicEdges[make_pair(edge0, edge2)];
+	BoundaryType type2 = provideConnectedPeriodicBoundary(type0, type1);
+
+ 	periodicEdges[perEdge0] = type2;
+ 	periodicEdges[perEdge1] = type2;
       }
     }
 
@@ -314,6 +311,59 @@ namespace AMDiS {
   }
 
 
+  BoundaryType ElementObjects::getNewBoundaryType()
+  {
+    FUNCNAME("ElementObjects::getNewBoundaryType()");
+
+    BoundaryType newPeriodicBoundaryType = 0;
+    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
+	 it != mesh->getPeriodicAssociations().end(); ++it)
+      newPeriodicBoundaryType = std::min(newPeriodicBoundaryType, it->first);
+    
+    TEST_EXIT_DBG(newPeriodicBoundaryType < 0)("Should not happen!\n");
+    newPeriodicBoundaryType--;
+    
+    mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
+      new VertexVector(feSpace->getAdmin(), "");
+    
+    return newPeriodicBoundaryType;
+  }
+
+
+  BoundaryType ElementObjects::provideConnectedPeriodicBoundary(BoundaryType b0, 
+								BoundaryType b1)
+  {
+    FUNCNAME("ElementObjects::provideConnectedPeriodicBoundary()");
+
+    std::pair<BoundaryType, BoundaryType> bConn = 
+      (b0 <= b1 ? std::make_pair(b0, b1) : std::make_pair(b1, b0));
+
+    if (bConnMap.count(bConn) == 0) {
+      BoundaryType newPeriodicBoundaryType = getNewBoundaryType();
+
+      VertexVector &vecB0 = mesh->getPeriodicAssociations(b0);
+      VertexVector &vecB1 = mesh->getPeriodicAssociations(b1);
+      VertexVector &vecC = mesh->getPeriodicAssociations(newPeriodicBoundaryType);
+
+      DOFIteratorBase it(const_cast<DOFAdmin*>(feSpace->getAdmin()), USED_DOFS);
+
+      for (it.reset(); !it.end(); ++it) {
+	if (!it.isDofFree()) {
+	  TEST_EXIT_DBG(vecB1[vecB0[it.getDOFIndex()]] == vecB0[vecB1[it.getDOFIndex()]])
+	    ("Should not happen!\n");
+
+	  vecC[it.getDOFIndex()] = vecB1[vecB0[it.getDOFIndex()]];
+	}
+      }
+      
+      bConnMap[bConn] = newPeriodicBoundaryType;
+    }
+
+
+    return bConnMap[bConn];
+  }
+
+
   void ElementObjects::createRankData(map<int, int>& macroElementRankMap) 
   {
     FUNCNAME("ElementObjects::createRankData()");
@@ -329,7 +379,7 @@ namespace AMDiS {
 	if (it2->elIndex > vertexInRank[it->first][elementInRank].elIndex)
 	  vertexInRank[it->first][elementInRank] = *it2;
 
-	vertexOwner[it->first] = max(vertexOwner[it->first], elementInRank);
+	vertexOwner[it->first] = std::max(vertexOwner[it->first], elementInRank);
       }
     }
     
@@ -345,7 +395,7 @@ namespace AMDiS {
 	if (it2->elIndex > edgeInRank[it->first][elementInRank].elIndex)
 	  edgeInRank[it->first][elementInRank] = *it2;
 
-	edgeOwner[it->first] = max(edgeOwner[it->first], elementInRank);
+	edgeOwner[it->first] = std::max(edgeOwner[it->first], elementInRank);
       }
     }
     
@@ -361,7 +411,7 @@ namespace AMDiS {
 	if (it2->elIndex > faceInRank[it->first][elementInRank].elIndex)
 	  faceInRank[it->first][elementInRank] = *it2;
 
-	faceOwner[it->first] = max(faceOwner[it->first], elementInRank);
+	faceOwner[it->first] = std::max(faceOwner[it->first], elementInRank);
       }
     }
   }

AMDiS/src/parallel/ElementObjectData.h

@@ -360,6 +360,21 @@ namespace AMDiS {
       return faceLocalMap[data];
     }
 
+    PerBoundMap<DegreeOfFreedom>::type& getPeriodicVertices()
+    {
+      return periodicVertices;
+    }
+
+    PerBoundMap<DofEdge>::type& getPeriodicEdges()
+    {
+      return periodicEdges;
+    }
+
+    PerBoundMap<DofFace>::type& getPeriodicFaces()
+    {
+      return periodicFaces;
+    }
+
     /// Write the element database to disk.
     void serialize(ostream &out);
     
@@ -400,6 +415,10 @@ namespace AMDiS {
       faceLocalMap[elObj] = face;
     }
 
+    BoundaryType provideConnectedPeriodicBoundary(BoundaryType b0, 
+						  BoundaryType b1);
+
+    BoundaryType getNewBoundaryType();
 
     /// Some auxiliary function to write the element object database to disk.
     void serialize(ostream &out, vector<ElementObjectData>& elVec);
@@ -480,7 +499,8 @@ namespace AMDiS {
     /// iterators, i.e., if no iteration is running.
     GeoIndex iterGeoPos;
 
-  public:
+    std::map<std::pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
+
     // The following three data structures store periodic DOFs, edges and faces.
     PerBoundMap<DegreeOfFreedom>::type periodicVertices;
     PerBoundMap<DofEdge>::type periodicEdges;

AMDiS/src/parallel/InteriorBoundary.cc

@@ -57,12 +57,17 @@ namespace AMDiS {
 	  TEST_EXIT_DBG(localDofs0[el0_v1] == localDofs1[el1_v0] ||
 			localDofs0[el0_v1] == localDofs1[el1_v1])
 	    ("This should not happen!\n");
-	  
+
 	  if (localDofs0[el0_v0] != localDofs1[el1_v0])
 	    otherMode = true; 	
 	} else {
 	  if (mesh->associated(localDofs0[el0_v0], localDofs1[el1_v0]) == false)
 	    otherMode = true;
+
+	  if ((elIndex == 28 && otherBound.elIndex == 41) ||
+	      (elIndex == 41 && otherBound.elIndex == 28)) {
+	    MSG("HERE TEST B (%d %d): %d\n", el0_v0, el1_v0, otherMode);
+	  }
 	}
       }
 

AMDiS/src/parallel/ParMetisPartitioner.cc

@@ -288,7 +288,7 @@ namespace AMDiS {
       if (elementInRank[index]) {
 	// get weight 
 	float wgt = static_cast<float>(elemWeights[index]);
-	maxWgt = max(wgt, maxWgt);
+	maxWgt = std::max(wgt, maxWgt);
 
 	// write float weight
 	TEST_EXIT_DBG(floatWgtsPos < floatWgts.size())("Should not happen!\n");
@@ -352,7 +352,7 @@ namespace AMDiS {
     int kpartMax = 0;
     for (int i = 0; i < nElements; i++)
       if (wgts[i] < kpart)
-	kpartMax = max(kpartMax, wgts[i]);
+	kpartMax = std::max(kpartMax, wgts[i]);
 
     mpi::globalMax(kpartMax);
 
@@ -363,8 +363,8 @@ namespace AMDiS {
     ssum = 0;
 
     for (int i = 0; i < nElements; i++) {
-      wgts[i] = min(wgts[i], kpartMax);
-      wgts[i] = max(wgts[i], 1);
+      wgts[i] = std::min(wgts[i], kpartMax);
+      wgts[i] = std::max(wgts[i], 1);
 
       smin = std::min(smin, wgts[i]);
       smax = std::max(smax, wgts[i]);

AMDiS/src/parallel/ParallelDebug.cc

@@ -143,6 +143,24 @@ namespace AMDiS {
   {
     FUNCNAME("ParallelDebug::testPeriodicBoundary()");
 
+
+    // === 1. check: All periodic DOFs should have at least a correct number ===
+    // === of periodic associations.                                         ===
+
+    for (map<int, std::set<BoundaryType> >::iterator it = pdb.periodicDofAssociations.begin();
+	 it != pdb.periodicDofAssociations.end(); ++it) {
+      WorldVector<double> c;
+      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, c);
+      int nAssoc = it->second.size();
+      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7))
+	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
+	 it->first, c[0], c[1], (pdb.mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
+    }    
+
+
+    // === 2. check: All periodic DOFs must be symmetric, i.e., if A is mapped ===
+    // === to B, then B must be mapped to A.                                   ===
+
     typedef MeshDistributor::PeriodicDofMap PeriodicDofMap;
 
     StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);
@@ -219,6 +237,80 @@ namespace AMDiS {
 
     mpi::globalAdd(foundError);
     TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
+
+
+    // === 3. check: On all edge and face periodic DOFs, at least on coordinate of ===
+    // === each periodic DOF pair must be equal (at least as long we consider      ===
+    // === periodic boundaries only on rectangulars and boxes.                     ===
+
+    RankToCoords sendCoords;
+    std::map<int, std::vector<BoundaryType> > rankToDofType;
+
+    for (InteriorBoundary::RankToBoundMap::iterator it = pdb.periodicBoundary.boundary.begin();
+	 it != pdb.periodicBoundary.boundary.end(); ++it) {
+      if (it->first == pdb.mpiRank)
+	continue;
+
+      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+	   boundIt != it->second.end(); ++boundIt) {	
+	if (boundIt->rankObj.subObj == VERTEX)
+	  continue;
+
+	DofContainer dofs;
+	boundIt->rankObj.el->getVertexDofs(pdb.feSpace, boundIt->rankObj, dofs);
+	boundIt->rankObj.el->getNonVertexDofs(pdb.feSpace, boundIt->rankObj, dofs);
+	
+	for (unsigned int i = 0; i < dofs.size(); i++) {
+	  WorldVector<double> c;
+	  pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpace, c);
+	  sendCoords[it->first].push_back(c);
+	  rankToDofType[it->first].push_back(boundIt->type);
+	}
+      }
+    }
+
+    // Each rank must receive exactly the same number of coordinates as it sends
+    // to another rank.
+    RankToCoords recvCoords;
+    for (RankToCoords::iterator it = sendCoords.begin(); 
+	 it != sendCoords.end(); ++it)
+      recvCoords[it->first].resize(it->second.size());
+
+
+    StdMpi<CoordsVec> stdMpiCoords(pdb.mpiComm, true);
+    stdMpiCoords.send(sendCoords);
+    stdMpiCoords.recv(recvCoords);   
+    stdMpiCoords.startCommunication<double>(MPI_DOUBLE);
+
+
+    for (RankToCoords::iterator it = sendCoords.begin(); 
+	 it != sendCoords.end(); ++it) {
+      for (unsigned int i = 0; i < it->second.size(); i++) {
+	WorldVector<double> &c0 = it->second[i];
+	WorldVector<double> &c1 = stdMpiCoords.getRecvData(it->first)[i];
+
+
+	int nEqual = 0;
+	for (int j = 0; j < pdb.mesh->getDim(); j++)
+	  if (c0[j] == c1[j])
+	    nEqual++;
+
+	if ((rankToDofType[it->first][i] >= -3 && nEqual < 2) ||
+	    (rankToDofType[it->first][i] < -3 && nEqual == 0)) {
+	  MSG("[DBG]  %d-ith periodic DOF in boundary between ranks %d <-> %d is not correct!\n",
+	      i, pdb.mpiRank, it->first);
+	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
+	      pdb.mpiRank, c0[0], c0[1], (pdb.mesh->getDim() == 3 ? c0[2] : 0.0));
+	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
+	      it->first, c1[0], c1[1], (pdb.mesh->getDim() == 3 ? c1[2] : 0.0));
+
+	  foundError = 1;
+	}
+      }
+    }
+
+    mpi::globalAdd(foundError);
+    TEST_EXIT(foundError == 0)("Wrond periodic coordinates found on at least on rank!\n");
   }
 
 
@@ -235,9 +327,6 @@ namespace AMDiS {
       return;
     }
 
-    /// Defines a mapping type from rank numbers to sets of coordinates.
-    typedef std::map<int, std::vector<WorldVector<double> > > RankToCoords;
-
     /// Defines a mapping type from rank numbers to sets of DOFs.
     typedef std::map<int, DofContainer> RankToDofContainer;
 
@@ -316,7 +405,6 @@ namespace AMDiS {
     // === So we can check if also the coordinates of the communicated DOFs are ===
     // === the same on both corresponding ranks.                                ===
 
-    typedef std::vector<WorldVector<double> > CoordsVec;
     StdMpi<CoordsVec> stdMpi(pdb.mpiComm, true);
     stdMpi.send(sendCoords);
     stdMpi.recv(recvCoords);   

AMDiS/src/parallel/ParallelDebug.h

@@ -32,6 +32,11 @@ namespace AMDiS {
   protected:
     typedef MeshDistributor::RankToDofContainer RankToDofContainer;
 
+    typedef std::vector<WorldVector<double> > CoordsVec;
+
+    /// Defines a mapping type from rank numbers to sets of coordinates.
+    typedef std::map<int, CoordsVec> RankToCoords;
+
   public:
     /** \brief
      * Tests the interior and the periodic boundaries on all ranks if their order

AMDiS/src/parallel/PetscSolver.cc

@@ -67,6 +67,8 @@ namespace AMDiS {
     std::vector<double> values;
     cols.reserve(300);
     values.reserve(300);
+    
+    std::vector<int> globalCols;
 
     // === Traverse all rows of the dof matrix and insert row wise the values ===
     // === to the petsc matrix.                                               ===
@@ -123,6 +125,7 @@ namespace AMDiS {
 	    // The col dof index is not periodic, simple add entry.
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
+	    globalCols.push_back(globalColDof);
 	  }
 	}
       }
@@ -136,7 +139,7 @@ namespace AMDiS {
 	std::set<int>& perAsc = meshDistributor->getPerDofAssociations(globalRowDof);
 
 	double scalFactor = 1.0 / (perAsc.size() + 1.0);
-
+	
 	for (unsigned int i = 0; i < values.size(); i++)
 	  values[i] *= scalFactor;	
 
@@ -150,13 +153,15 @@ namespace AMDiS {
 	  perCols.reserve(300);
 	  std::vector<double> perValues;
 	  perValues.reserve(300);
+
 	  for (unsigned int i = 0; i < cols.size(); i++) {
 	    int tmp = (cols[i] - dispAddCol) / dispMult;
 
-	    if (meshDistributor->isPeriodicDof(tmp, *perIt))
+	    if (meshDistributor->isPeriodicDof(tmp, *perIt)) {
 	      perCols.push_back((meshDistributor->getPeriodicMapping(*perIt, tmp) * dispMult) + dispAddCol);
-	    else
+	    } else {
 	      perCols.push_back(cols[i]);
+	    }
 	    
 	    perValues.push_back(values[i]);
 	  }