From f1d2d34b7cebf0ba4ef655e64ad5ee5d01d089c5 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 29 Jun 2009 15:04:52 +0000
Subject: [PATCH] Work on pdd.

---
 AMDiS/src/DOFMatrix.cc             |  15 +-
 AMDiS/src/ElementDofIterator.cc    |   3 +-
 AMDiS/src/ElementDofIterator.h     |  15 +-
 AMDiS/src/Mesh.cc                  |  13 +-
 AMDiS/src/ParallelDomainProblem.cc | 320 ++++++++++++++++++++++-------
 AMDiS/src/ParallelDomainProblem.h  |  72 ++++---
 6 files changed, 322 insertions(+), 116 deletions(-)

diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 847eb1bf..6de6a388 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -15,6 +15,8 @@
 #include "BoundaryManager.h"
 #include "Assembler.h"
 
+#include "mpi.h"
+
 namespace AMDiS {
 
   using namespace mtl;
@@ -221,18 +223,7 @@ namespace AMDiS {
 	for (int j = 0; j < nCol; j++) {  // for all columns
 	  DegreeOfFreedom col = colIndices[j];
 	  double entry = elMat[i][j];
-
-	  /*
-	  if (MPI::COMM_WORLD.Get_rank() == 0 && (row >= 156 || col >= 156))
-	    std::cout << "PROB 0: " << row << " " << col << std::endl;
-	  if (MPI::COMM_WORLD.Get_rank() == 1 && (row >= 151 || col >= 151))
-	    std::cout << "PROB 1: " << row << " " << col << std::endl;
-	  if (MPI::COMM_WORLD.Get_rank() == 2 && (row >= 146 || col >= 146))
-	    std::cout << "PROB 2: " << row << " " << col << std::endl;
-	  if (MPI::COMM_WORLD.Get_rank() == 3 && (row >= 153 || col >= 153))
-	    std::cout << "PROB 3: " << row << " " << col << std::endl;
-	  */
-
+  
 	  if (add)
 	    ins[row][col] += sign * entry;
 	  else
diff --git a/AMDiS/src/ElementDofIterator.cc b/AMDiS/src/ElementDofIterator.cc
index 338e37cd..78c1f3b6 100644
--- a/AMDiS/src/ElementDofIterator.cc
+++ b/AMDiS/src/ElementDofIterator.cc
@@ -82,5 +82,6 @@ namespace AMDiS {
     }
 
     return true;
-  } 
+  }
+
 }
diff --git a/AMDiS/src/ElementDofIterator.h b/AMDiS/src/ElementDofIterator.h
index 60e376f1..7b571f79 100644
--- a/AMDiS/src/ElementDofIterator.h
+++ b/AMDiS/src/ElementDofIterator.h
@@ -66,6 +66,19 @@ namespace AMDiS {
       return &dofs[node0 + elementPos][n0 + dofPos];
     }
 
+    /// Returns \ref pos, the current position (vertex, edge, face) of the traverse.
+    int getCurrentPos()
+    {
+      return pos;
+    }
+
+    /// Returns \ref elementPos, the number of vertex, edge or face that is traversed.
+    int getCurrentElementPos()
+    {
+      return elementPos;
+    }
+
+      
   protected:
     /// The dof admin for which dofs should be traversed.
     const DOFAdmin* admin;
@@ -108,7 +121,7 @@ namespace AMDiS {
      */
     int nElements;
 
-    /// Current element (vertex, edge, face) that is traversed.
+    /// Current element, i.e., ith vertex, edge or face, that is traversed.
     int elementPos;
 
     /// Currrent dof that is traversed on the current element;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index b3093503..ce15ab50 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -8,6 +8,7 @@
 #include "AdaptInstationary.h"
 #include "FiniteElemSpace.h"
 #include "ElementData.h"
+#include "ElementDofIterator.h"
 #include "MacroElement.h"
 #include "MacroReader.h"
 #include "Mesh.h"
@@ -867,19 +868,23 @@ namespace AMDiS {
   {
     DimVec<double>* baryCoords;
     bool found = false;
+    ElementDofIterator elDofIter(feSpace->getAdmin());
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(this, -1, 
 					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
     while (elInfo) {
-      for (int i = 0; i < nDOFEl; i++) {
-	if (elInfo->getElement()->getDOF(i) == dof) {
+      elDofIter.reset(elInfo->getElement());
+      int i = 0;
+      do {
+	if (elDofIter.getDofPtr() == dof) {
 	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
 	  elInfo->coordToWorld(*baryCoords, coords);
 	  found = true;
 	  break;
 	}
-      }
-      
+	i++;
+      } while (elDofIter.next());
+
       if (found)
 	break;
 	  
diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc
index e457e100..1427fb61 100644
--- a/AMDiS/src/ParallelDomainProblem.cc
+++ b/AMDiS/src/ParallelDomainProblem.cc
@@ -38,6 +38,13 @@ namespace AMDiS {
       initialPartitionMesh(true),
       nRankDOFs(0)
   {
+    FUNCNAME("ParallelDomainBase::ParalleDomainBase()");
+
+    TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
+      ("Only meshes with one DOFAdmin are supported!\n");
+    TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
+      ("Wrong pre dof number for DOFAdmin!\n");
+
     mpiRank = MPI::COMM_WORLD.Get_rank();
     mpiSize = MPI::COMM_WORLD.Get_size();
     mpiComm = MPI::COMM_WORLD;
@@ -100,8 +107,6 @@ namespace AMDiS {
       admin.setFirstHole(mapLocalGlobalDOFs.size());
     }
 
-    exit(0);
-
     // === Global refinements. ===
 
     int globalRefinement = 0;
@@ -114,9 +119,9 @@ namespace AMDiS {
     }
 
 
-// #if (DEBUG != 0)
-//     testInteriorBoundary();
-// #endif
+#if (DEBUG != 0)
+    testInteriorBoundary();
+#endif
 
     // === Create petsc matrix. ===
 
@@ -182,7 +187,7 @@ namespace AMDiS {
 	     icend = end<nz>(cursor); icursor != icend; ++icursor)
 	if (value(*icursor) != 0.0) {
 	  int r = mapLocalGlobalDOFs[row(*icursor)];
-	  int c = mapLocalGlobalDOFs[col(*icursor)];
+	  int c = mapLocalGlobalDOFs[col(*icursor)];  
 	  double v = value(*icursor);
 
 	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
@@ -234,9 +239,8 @@ namespace AMDiS {
     int requestCounter = 0;
     
     int i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   sendIt = sendDofs.begin();
-	 sendIt != sendDofs.end();
+    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
+	 sendIt != sendDofs.end(); 
 	 ++sendIt, i++) {
       int nSendDOFs = sendIt->second.size();
       sendBuffers[i] = new double[nSendDOFs];
@@ -249,9 +253,8 @@ namespace AMDiS {
     }
 
     i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   recvIt = recvDofs.begin();
-	 recvIt != recvDofs.end();
+    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
+	 recvIt != recvDofs.end(); 
 	 ++recvIt, i++) {
       int nRecvDOFs = recvIt->second.size();
       recvBuffers[i] = new double[nRecvDOFs];
@@ -265,12 +268,12 @@ namespace AMDiS {
 
     
     i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   recvIt = recvDofs.begin();
+    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
 	 recvIt != recvDofs.end();
 	 ++recvIt, i++) {
-      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
-	(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
+      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
+	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
+      }
 
       delete [] recvBuffers[i];
     }
@@ -330,8 +333,8 @@ namespace AMDiS {
   }
 
 
-  void ParallelDomainBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
-						      std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
+  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs,
+						      DofToRank& boundaryDOFs)
   {
     FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
 
@@ -499,14 +502,14 @@ namespace AMDiS {
   }
 
 
-  void ParallelDomainBase::createLocalGlobalNumbering(
-                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
-				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
-				 int& nRankDOFs, 
-				 int& nOverallDOFs)
+  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
+						      DofToRank& boundaryDOFs,
+						      int& nRankDOFs, 
+						      int& nOverallDOFs)
   {
     FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
 
+
     // === Get rank information about DOFs. ===
 
     // Stores to each DOF pointer the set of ranks the DOF is part of.
@@ -537,7 +540,7 @@ namespace AMDiS {
     // another rank.
     std::map<int, int> recvNewDofs;
 
-    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
+    for (DofToRank::iterator it = boundaryDOFs.begin();
 	 it != boundaryDOFs.end();
 	 ++it) {
 
@@ -570,7 +573,6 @@ namespace AMDiS {
       }
     }
 
-
     // === Send and receive the dof indices at boundary. ===
 
     std::vector<int*> sendBuffers(sendNewDofs.size());
@@ -646,6 +648,20 @@ namespace AMDiS {
 
     // === Change dof indices at boundary from other ranks. ===
 
+    // Within this small data structure we track which dof index was already changed.
+    // This is used to avoid the following situation: Assume, there are two dof indices
+    // a and b in boundaryDOFs. Then we have to change index a to b and b to c. When
+    // the second rule applies, we have to avoid that not the first b, resulted from
+    // changing a to b, is set to c, but the second one. Therefore, after the first
+    // rule was applied, the dof pointer is set to false in this data structure and 
+    // is not allowed to be changed anymore.
+    std::map<const DegreeOfFreedom*, bool> dofChanged;
+    for (DofToRank::iterator dofIt = boundaryDOFs.begin();
+	 dofIt != boundaryDOFs.end();
+	 ++dofIt) {
+      dofChanged[dofIt->first] = false;
+    }
+
     i = 0;
     for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
 	 recvIt != recvNewDofs.end();
@@ -659,13 +675,18 @@ namespace AMDiS {
 
 	bool found = false;
 
-	for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDOFs.begin();
+	// Iterate over all boundary dofs to find the dof, which index we have to change.
+
+	for (DofToRank::iterator dofIt = boundaryDOFs.begin();
 	     dofIt != boundaryDOFs.end();
 	     ++dofIt) {
 
-	  if ((dofIt->first)[0] == oldDof) {
+	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
+	    dofChanged[dofIt->first] = true;
+
 	    recvDofs[recvIt->first].push_back(dofIt->first);
-	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
+	    *(const_cast<DegreeOfFreedom*>(dofIt->first)) = newLocalDof;
+
 	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
 	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
 	    isRankDOF[newLocalDof] = false;
@@ -688,27 +709,31 @@ namespace AMDiS {
     FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
 
     std::set<const DegreeOfFreedom*> rankDOFs;
-    std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
+    DofToRank newBoundaryDOFs;
     
+
     // === Get all DOFs in ranks partition. ===
 
+    ElementDofIterator elDofIt(feSpace->getAdmin());
+
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
     while (elInfo) {
       Element *element = elInfo->getElement();
       
-      for (int i = 0; i < 3; i++) 
-	rankDOFs.insert(element->getDOF(i));
+      elDofIt.reset(element);
+      do {
+	rankDOFs.insert(elDofIt.getDofPtr());
+      } while(elDofIt.next());
 
       elInfo = stack.traverseNext(elInfo);
     }
 
-
     // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
     // === rankDOFs to boundaryDOFs.                                               ===
 
-    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
-    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
+    RankToDofContainer sendNewDofs;
+    RankToDofContainer recvNewDofs;
 
     for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
 	   myIntBoundary.boundary.begin();
@@ -744,10 +769,22 @@ namespace AMDiS {
 
 	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
 	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
+	
+	if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof1) ==
+	    sendNewDofs[it->first].end())
+	  sendNewDofs[it->first].push_back(dof1);
+	if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof2) ==
+	    sendNewDofs[it->first].end())
+	  sendNewDofs[it->first].push_back(dof2);
+	
+	DofContainer boundDOFs;
+	addAllVertexDOFs(boundIt->rankObject.el, 
+			 boundIt->rankObject.ithObjAtBoundary,
+			 boundDOFs);
+	addAllEdgeDOFs(boundIt->rankObject.el, 
+		       boundIt->rankObject.ithObjAtBoundary,
+		       boundDOFs);
 
-	std::vector<const DegreeOfFreedom*> boundDOFs;
-	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
-		   boundDOFs);
 	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
 	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
 	  sendNewDofs[it->first].push_back(boundDOFs[i]);
@@ -755,7 +792,6 @@ namespace AMDiS {
       }
     }    
 
-
     for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
 	   otherIntBoundary.boundary.begin();
 	 it != otherIntBoundary.boundary.end();
@@ -793,9 +829,21 @@ namespace AMDiS {
 	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
 	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
 
-	std::vector<const DegreeOfFreedom*> boundDOFs;
-	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
-		   boundDOFs);
+	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof2) ==
+	    recvNewDofs[it->first].end())
+	  recvNewDofs[it->first].push_back(dof2);
+	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof1) ==
+	    recvNewDofs[it->first].end())
+	  recvNewDofs[it->first].push_back(dof1);
+
+	DofContainer boundDOFs;
+	addAllEdgeDOFs(boundIt->rankObject.el, 
+		       boundIt->rankObject.ithObjAtBoundary,
+		       boundDOFs);
+	addAllVertexDOFs(boundIt->rankObject.el, 
+			 boundIt->rankObject.ithObjAtBoundary,
+			 boundDOFs);
+
 	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
 	//	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
 	  rankDOFs.erase(boundDOFs[i]);
@@ -803,7 +851,6 @@ namespace AMDiS {
 	  recvNewDofs[it->first].push_back(boundDOFs[i]);
 	}
       }
-
     }
 
     nRankDOFs = rankDOFs.size();
@@ -846,14 +893,13 @@ namespace AMDiS {
     int requestCounter = 0;
 
     i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   sendIt = sendNewDofs.begin();
+    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
 	 sendIt != sendNewDofs.end(); 
 	 ++sendIt, i++) {
       int nSendDofs = sendIt->second.size();
       sendBuffers[i] = new int[nSendDofs];
       int c = 0;
-      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = sendIt->second.begin();
+      for (DofContainer::iterator dofIt = sendIt->second.begin();
 	   dofIt != sendIt->second.end();
 	   ++dofIt)
 	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
@@ -863,9 +909,9 @@ namespace AMDiS {
     }
 
     i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   recvIt = recvNewDofs.begin();
-	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
+    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end(); 
+	 ++recvIt, i++) {
       int nRecvDofs = recvIt->second.size();
       recvBuffers[i] = new int[nRecvDofs];
 	
@@ -873,28 +919,25 @@ namespace AMDiS {
 	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
     }
 
-
     MPI::Request::Waitall(requestCounter, request);
 
 
     i = 0;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   sendIt = sendNewDofs.begin();
+    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
 	 sendIt != sendNewDofs.end(); 
 	 ++sendIt)
       delete [] sendBuffers[i++];
 
     i = 0;
     int newDofIndex = nRankDOFs;
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   recvIt = recvNewDofs.begin();
-	 recvIt != recvNewDofs.end(); ++recvIt) {
-      
+    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end(); 
+	 ++recvIt) {      
       int j = 0;
-      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = recvIt->second.begin();
+      for (DofContainer::iterator dofIt = recvIt->second.begin();
 	   dofIt != recvIt->second.end();
 	   ++dofIt) {
-	const_cast<DegreeOfFreedom*>(*dofIt)[0] = newDofIndex;
+	(*const_cast<DegreeOfFreedom*>(*dofIt)) = newDofIndex;
 	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
 	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
 	isRankDOF[newDofIndex] = false;
@@ -912,31 +955,31 @@ namespace AMDiS {
   }
 
 
-  void ParallelDomainBase::addAllDOFs(Element *el, int ithEdge, 
-				      std::vector<const DegreeOfFreedom*>& dofs)
+  void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge, 
+					    DofContainer& dofs)
   {
-    FUNCNAME("ParallelDomainBase::addAllDOFs()");
+    FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
 
     switch (ithEdge) {
     case 0:
       if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
-	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
+	addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
 	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
-	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
+	addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
       }
       break;
     case 1:
       if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
-	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
+	addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
 	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
-	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
+	addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
       }
       break;
     case 2:      
       if (el->getFirstChild()) {
-	addAllDOFs(el->getFirstChild(), 0, dofs);
+	addAllVertexDOFs(el->getFirstChild(), 0, dofs);
 	dofs.push_back(el->getFirstChild()->getDOF(2));
-	addAllDOFs(el->getSecondChild(), 1, dofs);
+	addAllVertexDOFs(el->getSecondChild(), 1, dofs);
       }
       break;      
     default:
@@ -945,10 +988,57 @@ namespace AMDiS {
   }
 
 
-  void ParallelDomainBase::createDOFMemberInfo(
-		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
-		       std::vector<const DegreeOfFreedom*>& rankDOFs,
-		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
+  void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge,
+					  DofContainer& dofs)
+  {
+    FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()");
+
+    bool addThisEdge = false;
+
+    switch (ithEdge) {
+    case 0:
+      if (el->getSecondChild())
+	addAllEdgeDOFs(el->getSecondChild(), 2, dofs);
+      else
+	addThisEdge = true;
+
+      break;
+    case 1:
+      if (el->getFirstChild())
+	addAllEdgeDOFs(el->getFirstChild(), 2, dofs);
+      else
+	addThisEdge = true;
+
+      break;
+    case 2:
+      if (el->getFirstChild()) {
+	addAllEdgeDOFs(el->getFirstChild(), 1, dofs);
+	addAllEdgeDOFs(el->getFirstChild(), 0, dofs);
+      } else {
+	addThisEdge = true;
+      } 
+
+      break;
+    default:
+      ERROR_EXIT("Should never happen!\n");
+    }
+
+    if (addThisEdge) {
+      ElementDofIterator elDofIter(feSpace->getAdmin());
+      elDofIter.reset(el);
+      do {
+	if (elDofIter.getCurrentPos() == 1 &&
+	    elDofIter.getCurrentElementPos() == ithEdge) {
+	  dofs.push_back(elDofIter.getDofPtr());
+	}
+      } while(elDofIter.next());      
+    }
+  }
+
+
+  void ParallelDomainBase::createDOFMemberInfo(DofToPartitions& partitionDOFs,
+					       DofContainer& rankDOFs,
+					       DofToRank& boundaryDOFs)
   {
     // === Determine to each dof the set of partitions the dof belongs to. ===
 
@@ -972,7 +1062,7 @@ namespace AMDiS {
     // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
 
     // iterate over all DOFs
-    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
+    for (DofToPartitions::iterator it = partitionDOFs.begin();
 	 it != partitionDOFs.end();
 	 ++it) {
 
@@ -1052,13 +1142,99 @@ namespace AMDiS {
       }
     }
 
+    std::vector<int> sendSize(mpiSize, 0);
+    std::vector<int> recvSize(mpiSize, 0);
+    std::vector<int> recvSizeBuffer(mpiSize, 0);
+    MPI::Request request[(mpiSize - 1) * 2];
+    int requestCounter = 0;
+
     for (std::map<int, std::vector<WorldVector<double> > >::iterator it 
-	   = sendCoords.begin();
+	   = sendCoords.begin(); 
 	 it != sendCoords.end(); ++it) {
+      sendSize[it->first] = it->second.size();
     }
 
     for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin();
 	 it != recvCoords.end(); ++it) {
+      recvSize[it->first] = it->second.size();
+    }
+
+    for (int i = 0; i < mpiSize; i++) {
+      if (i == mpiRank)
+	continue;
+
+      request[requestCounter++] = mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
+    }   
+
+    for (int i = 0; i < mpiSize; i++) {
+      if (i == mpiRank)
+	continue;
+
+      request[requestCounter++] = mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
+    }
+
+    MPI::Request::Waitall(requestCounter, request);
+
+
+    for (int i = 0; i < mpiSize; i++) {
+      if (i == mpiRank)
+	continue;
+
+      if (recvSize[i] != recvSizeBuffer[i]) {
+	std::cout << "Error: MPI rank " << mpiRank << " expectes to receive " 
+		  << recvSize[i] << " DOFs from rank " << i << ". But this rank sends "
+		  << recvSizeBuffer[i] << " DOFs!" << std::endl;
+
+	ERROR_EXIT("Should not happen!\n");
+      }
+    }
+
+    std::map<int, double*> sendCoordsBuffer;
+    std::map<int, double*> recvCoordsBuffer;
+    requestCounter = 0;
+
+    for (std::map<int, std::vector<WorldVector<double> > >::iterator it 
+	   = sendCoords.begin(); 
+	 it != sendCoords.end(); ++it) {
+      sendCoordsBuffer[it->first] = new double[it->second.size() * 2];
+
+      for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
+	sendCoordsBuffer[it->first][i * 2] = (it->second)[i][0];
+	sendCoordsBuffer[it->first][i * 2 + 1] = (it->second)[i][1];
+      }
+
+      request[requestCounter++] = mpiComm.Isend(sendCoordsBuffer[it->first],
+						it->second.size() * 2,
+						MPI_DOUBLE, it->first, 0);
+    }
+
+    for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin();
+	 it != recvCoords.end(); ++it) {
+      recvCoordsBuffer[it->first] = new double[it->second.size() * 2];
+
+      request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first], 
+						it->second.size() * 2,
+						MPI_DOUBLE, it->first, 0);
+    }
+
+    MPI::Request::Waitall(requestCounter, request);    
+
+    for (std::map<int, std::vector<WorldVector<double> > >::iterator it 
+	   = sendCoords.begin(); 
+	 it != sendCoords.end(); ++it) {
+      delete [] sendCoordsBuffer[it->first];
+    }
+
+    for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin();
+	 it != recvCoords.end(); ++it) {
+      for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
+
+	if ((it->second)[i][0] != recvCoordsBuffer[it->first][i * 2] ||
+	    (it->second)[i][1] != recvCoordsBuffer[it->first][i * 2 + 1])
+	  ERROR_EXIT("WRONG COORDS!\n");
+      } 
+
+      delete [] recvCoordsBuffer[it->first];
     }
   }
 
diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h
index 56d3d030..462c189a 100644
--- a/AMDiS/src/ParallelDomainProblem.h
+++ b/AMDiS/src/ParallelDomainProblem.h
@@ -45,12 +45,31 @@ namespace AMDiS {
   class ParallelDomainBase : public ProblemIterationInterface,
                                     public ProblemTimeInterface
   {
+  private:
+    /// Defines type for a vector of DOFs.
+    typedef std::vector<const DegreeOfFreedom*> DofContainer;
+
+    /// Defines a mapping type from DOFs to rank numbers.
+    typedef std::map<const DegreeOfFreedom*, int> DofToRank;
+
+    /// Defines a mapping type from DOFs to a set of rank numbers.
+    typedef std::map<const DegreeOfFreedom*, std::set<int> > DofToPartitions;
+
+    /// Defines a mapping type from rank numbers to sets of DOFs.
+    typedef std::map<int, DofContainer> RankToDofContainer;
+
+    /// Defines a mapping type from DOF indices to DOF indices.
+    typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
+
+    /// Defines a mapping type from DOF indices to boolean values.
+    typedef std::map<DegreeOfFreedom, bool> DofToBool;
+
   public:
     ParallelDomainBase(const std::string& name,
-			      ProblemIterationInterface *iterationIF,
-			      ProblemTimeInterface *timeIF,
-			      FiniteElemSpace *feSpace,
-			      RefinementManager *refineManager);
+		       ProblemIterationInterface *iterationIF,
+		       ProblemTimeInterface *timeIF,
+		       FiniteElemSpace *feSpace,
+		       RefinementManager *refineManager);
 
     virtual ~ParallelDomainBase() {}
 
@@ -132,8 +151,8 @@ namespace AMDiS {
      * Determine the interior boundaries, i.e. boundaries between ranks, and store
      * all information about them in \ref interiorBoundary.
      */
-    void createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
-				    std::map<const DegreeOfFreedom*, int>& boundaryDOFs);
+    void createInteriorBoundaryInfo(DofContainer& rankDOFs,
+				    DofToRank& boundaryDOFs);
 
     /// Removes all macro elements from the mesh that are not part of ranks partition.
     void removeMacroElements();
@@ -142,23 +161,24 @@ namespace AMDiS {
     /** \brief
      * Creates from a macro mesh a correct local and global DOF index numbering.
      *
-     * @param[out] rankDOFs      Returns all DOFs from the macro mesh, which are owned
+     * \param[out] rankDOFs      Returns all DOFs from the macro mesh, which are owned
      *                           by the rank after partitioning the macro mesh.
-     * @param[out] boundaryDOFs  Returns all DOFs from the macro mesh, which lies at
+     * \param[out] boundaryDOFs  Returns all DOFs from the macro mesh, which lies at
      *                           an interior boundary of the rank. This object maps
      *                           each such DOF to the rank that owns this DOF.
-     * @param[out] nRankDOFs     Number of DOFs owned by rank.
-     * @param[out] nOverallDOFs  Number of all DOFs in macro mesh.
+     * \param[out] nRankDOFs     Number of DOFs owned by rank.
+     * \param[out] nOverallDOFs  Number of all DOFs in macro mesh.
      */
-    void createLocalGlobalNumbering(std::vector<const DegreeOfFreedom*>& rankDOFs,
-				    std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
+    void createLocalGlobalNumbering(DofContainer& rankDOFs,
+				    DofToRank& boundaryDOFs,
 				    int& nRankDOFs, 
 				    int& nOverallDOFs);
 
     void updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs);
 
-    void addAllDOFs(Element *el, int ithEdge, 
-		    std::vector<const DegreeOfFreedom*>& dofs);
+    void addAllVertexDOFs(Element *el, int ithEdge, DofContainer& dofs);
+
+    void addAllEdgeDOFs(Element *el, int ithEdge, DofContainer& dofs);
 
     /** \brief
      * This function traverses the whole mesh, i.e. before it is really partitioned,
@@ -166,15 +186,15 @@ namespace AMDiS {
      * be used, if \ref partitionVec is set correctly. This is only the case, when
      * the macro mesh is partitioned.
      *
-     * @param[out] partionDOFs   Stores to each DOF pointer the set of ranks the DOF is
+     * \param[out] partionDOFs   Stores to each DOF pointer the set of ranks the DOF is
      *                           part of.
-     * @param[out] rankDOFs      Stores all rank DOFs.
-     * @param[out] boundaryDOFs  Stores all DOFs in ranks partition that are on an 
+     * \param[out] rankDOFs      Stores all rank DOFs.
+     * \param[out] boundaryDOFs  Stores all DOFs in ranks partition that are on an 
      *                           interior boundary but correspond to another rank.
      */
-    void createDOFMemberInfo(std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
-			     std::vector<const DegreeOfFreedom*>& rankDOFs,
-			     std::map<const DegreeOfFreedom*, int>& boundaryDOFs);
+    void createDOFMemberInfo(DofToPartitions& partitionDOFs,
+			     DofContainer& rankDOFs,
+			     DofToRank& boundaryDOFs);
 
      
     /** \brief
@@ -254,7 +274,7 @@ namespace AMDiS {
      * Set of all interior boundary DOFs in ranks partition. The object maps to
      * each such DOF to the number of the rank that owns this DOF.
      */
-    std::map<const DegreeOfFreedom*, int> boundaryDOFs;
+    DofToRank boundaryDOFs;
 
     /** \brief 
      * Defines the interior boundaries of the domain that result from partitioning
@@ -276,26 +296,26 @@ namespace AMDiS {
      * This map contains for each rank the list of dofs the current rank must send
      * to exchange solution dofs at the interior boundaries.
      */
-    std::map<int, std::vector<const DegreeOfFreedom*> > sendDofs;
+    RankToDofContainer sendDofs;
 
     /** \brief
      * This map contains for each rank the list of dofs from which the current rank 
      * must receive solution values of dofs at the interior boundaries.
      */
-    std::map<int, std::vector<const DegreeOfFreedom*> > recvDofs;
+    RankToDofContainer recvDofs;
 
     /// Maps local to global dof indices.
-    std::map<DegreeOfFreedom, DegreeOfFreedom> mapLocalGlobalDOFs;
+    DofMapping mapLocalGlobalDOFs;
 
     /// Maps global to local dof indices.
-    std::map<DegreeOfFreedom, DegreeOfFreedom> mapGlobalLocalDOFs;
+    DofMapping mapGlobalLocalDOFs;
 
     /** \brief
      * Maps all DOFs in ranks partition to a bool value. If it is true, the DOF is
      * owned by the rank. Otherwise, its an interior boundary DOF that is owned by
      * another rank.
      */
-    std::map<DegreeOfFreedom, bool> isRankDOF;
+    DofToBool isRankDOF;
   };
 
   class ParallelDomainScal : public ParallelDomainBase
-- 
GitLab