From ed1ab0230d5bf817e01f753e5d5f8ef976d21b5a Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 13 Aug 2009 13:50:02 +0000
Subject: [PATCH] Bugfix for assembling with 8 or more processors.

---
 AMDiS/libtool                   |   6 +-
 AMDiS/src/ElementDofIterator.cc |   2 +-
 AMDiS/src/Lagrange.cc           |  25 ++----
 AMDiS/src/ParallelDomainBase.cc | 150 ++++++++++++++++++++++++--------
 AMDiS/src/ParallelDomainBase.h  |  12 ++-
 5 files changed, 135 insertions(+), 60 deletions(-)

diff --git a/AMDiS/libtool b/AMDiS/libtool
index 4c12058b..f8592f9b 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
 
 # ### BEGIN LIBTOOL CONFIG
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host p2d087:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
 # End:
 # ### BEGIN LIBTOOL TAG CONFIG: CXX
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host p2d087:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -7065,7 +7065,7 @@ include_expsyms=""
 
 # ### BEGIN LIBTOOL TAG CONFIG: F77
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host p2d087:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
diff --git a/AMDiS/src/ElementDofIterator.cc b/AMDiS/src/ElementDofIterator.cc
index d8dc65ca..4e44a572 100644
--- a/AMDiS/src/ElementDofIterator.cc
+++ b/AMDiS/src/ElementDofIterator.cc
@@ -42,7 +42,7 @@ namespace AMDiS {
 
   bool ElementDofIterator::next()
   {
-    // First iteratore over the dofs of one element (vertex, edge, face).
+    // First iterate over the dofs of one element (vertex, edge, face).
     dofPos++;
 
     if (dofPos >= nDofs) {
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index e13fbca3..4a9857d9 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -2982,24 +2982,17 @@ namespace AMDiS {
     if (n < 1) 
       return;
 
-    const Element *el;
-    const DegreeOfFreedom *cd;
-    DegreeOfFreedom  pd[19], cdi;
-    int node0, n0, i, typ, lr_set;
-    const DOFAdmin *admin;
-
-    el = list->getElement(0);
-    typ = list->getType(0);
-
-    admin = drv->getFESpace()->getAdmin();
-
+    DegreeOfFreedom pd[20], cdi;
+    const Element *el = list->getElement(0);
+    int typ = list->getType(0);
+    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
     basFct->getLocalIndices(el, admin, pd);
 
     /****************************************************************************/
     /*  values on child[0]                                                      */
     /****************************************************************************/
 
-    cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+    const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
 
     (*drv)[pd[0]] += 
       (0.0625*(-(*drv)[cd[3]] + (*drv)[cd[12]] + (*drv)[cd[14]]
@@ -3087,15 +3080,15 @@ namespace AMDiS {
     /*   adjust neighbour values                                                */
     /****************************************************************************/
 
-    node0 = drv->getFESpace()->getMesh()->getNode(FACE);
-    n0 = admin->getNumberOfPreDOFs(FACE);
+    int node0 = drv->getFESpace()->getMesh()->getNode(FACE);
+    int n0 = admin->getNumberOfPreDOFs(FACE);
 
-    for (i = 1; i < n; i++) {
+    for (int i = 1; i < n; i++) {
       el = list->getElement(i);
       typ = list->getType(i);
       basFct->getLocalIndices(el, admin, pd);
 
-      lr_set = 0;
+      int lr_set = 0;
       if (list->getNeighbourElement(i, 0) &&  list->getNeighbourNr(i, 0) < i)
 	lr_set = 1;
 
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index 520bbf77..73f31594 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -320,14 +320,22 @@ namespace AMDiS {
 
 	    if (isRankDof[*cursor]) {
 	      r -= rstart * nComponents;
-	
-	      TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n");
+
+#if (DEBUG != 0)    
+	      if (r < 0 || r >= nRankRows) {
+		std::cout << "ERROR in rank: " << mpiRank << std::endl;
+		std::cout << "  Wrong r = " << r << " " << *cursor << " " 
+			  << mapLocalGlobalDOFs[*cursor] << " " 
+			  << nRankRows << std::endl;
+		ERROR_EXIT("Should not happen!\n");
+	      }
+#endif
 	      
 	      for (icursor_type icursor = begin<nz>(cursor), 
 		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
 		if (value(*icursor) != 0.0) {
 		  int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;
-		  
+	  
 		  if (c >= rstart * nComponents && 
 		      c < rstart * nComponents + nRankRows)
 		    d_nnz[r]++;
@@ -409,10 +417,9 @@ namespace AMDiS {
 	 it != recvDofs.end(); ++it) {
       int nSend = sendMatrixEntry[it->first].size();
 
-      if (nSend > 0) {
+      if (nSend > 0)
 	request[requestCounter++] = 
 	  mpiComm.Isend(sendBuffers[i], nSend * 2, MPI_INT, it->first, 0);
-      }
 
       i++;
     }
@@ -588,7 +595,8 @@ namespace AMDiS {
     KSPCreate(PETSC_COMM_WORLD, &solver);
     KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
 
-    KSPSetTolerances(solver, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
+    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
+
     KSPSetType(solver, KSPBCGS);
     KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
     KSPSetFromOptions(solver);
@@ -912,7 +920,7 @@ namespace AMDiS {
     DofContainer rankAllDofs;
     DofToRank boundaryDofs;
 
-    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
+    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs, vertexDof);
 
     nRankDOFs = rankDOFs.size();
     nOverallDOFs = partitionDOFs.size();
@@ -966,7 +974,7 @@ namespace AMDiS {
 
     // Stores to each rank a map from DOF pointers (which are all owned by the rank
     // and lie on an interior boundary) to new global DOF indices.
-    std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> > sendNewDofs;
+    std::map<int, DofIndexMap > sendNewDofs;
 
     // Maps to each rank the number of new DOF indices this rank will receive from.
     // All these DOFs are at an interior boundary on this rank, but are owned by
@@ -1011,18 +1019,14 @@ namespace AMDiS {
     int requestCounter = 0;
 
     i = 0;
-    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
-	   sendIt = sendNewDofs.begin();
-	 sendIt != sendNewDofs.end();
-	 ++sendIt, i++) {
+    for (std::map<int, DofIndexMap >::iterator sendIt = sendNewDofs.begin();
+	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
       int nSendDofs = sendIt->second.size() * 2;
       sendBuffers[i] = new int[nSendDofs];
 
       int c = 0;
-      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
-	     dofIt = sendIt->second.begin();
-	   dofIt != sendIt->second.end();
-	   ++dofIt) {
+      for (DofIndexMap::iterator dofIt = sendIt->second.begin();
+	   dofIt != sendIt->second.end(); ++dofIt) {
 	sendBuffers[i][c++] = *(dofIt->first);
 	sendBuffers[i][c++] = dofIt->second;
 
@@ -1035,8 +1039,7 @@ namespace AMDiS {
 
     i = 0;
     for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
-	 recvIt != recvNewDofs.end();
-	 ++recvIt, i++) {
+	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
       int nRecvDOFs = recvIt->second * 2;
       recvBuffers[i] = new int[nRecvDOFs];
 
@@ -1142,9 +1145,12 @@ namespace AMDiS {
     // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
     // === rankDOFs to boundaryDOFs.                                               ===
 
+    RankToDofContainer oldSendDofs = sendDofs;
+    RankToDofContainer oldRecvDofs = recvDofs;
+
     sendDofs.clear();
     recvDofs.clear();
-
+    
     for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
 	 it != myIntBoundary.boundary.end(); ++it) {
 
@@ -1154,6 +1160,11 @@ namespace AMDiS {
 	DofContainer dofs;
 	DofContainer &dofsToSend = sendDofs[it->first];
 
+	for (DofContainer::iterator iit = oldSendDofs[it->first].begin(); 
+	     iit != oldSendDofs[it->first].end(); ++iit)
+	  if (vertexDof[*iit])
+	    dofsToSend.push_back(*iit);
+	
 	switch (boundIt->rankObject.ithObjAtBoundary) {
 	case 0:
 	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
@@ -1171,21 +1182,27 @@ namespace AMDiS {
 	  ERROR_EXIT("Should never happen!\n");
 	}
 
+	/*
 	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt)
 	  if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
-	    dofsToSend.push_back(*dofIt);	
+	    if (find(oldSendDofs[it->first].begin(), oldSendDofs[it->first].end(), *dofIt) != oldSendDofs[it->first].end())      
+	      dofsToSend.push_back(*dofIt);	
+	*/
 
 	dofs.clear();
 	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
 			 dofs);	
   	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
   		       dofs);
-	for (int i = 0; i < static_cast<int>(dofs.size()); i++)
-	  dofsToSend.push_back(dofs[i]);
+	for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
+	  TEST_EXIT_DBG(find(dofsToSend.begin(), dofsToSend.end(), dofs[i]) == dofsToSend.end())
+	    ("Should not happen!\n");
+
+	  dofsToSend.push_back(dofs[i]);	
+	}
       }
     }
-
-    
+   
 
     for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
 	 it != otherIntBoundary.boundary.end(); ++it) {
@@ -1196,7 +1213,26 @@ namespace AMDiS {
 	DofContainer dofs;
 	DofContainer &dofsToRecv = recvDofs[it->first];
 
-	switch (boundIt->rankObject.ithObjAtBoundary) {
+	/*
+	for (DofContainer::iterator iit = dofsToRecv.begin(); iit != dofsToRecv.end(); ++iit) {
+	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *iit);
+	  if (eraseIt != rankDOFs.end())
+	    rankDOFs.erase(eraseIt);    
+	}
+	*/
+
+	for (DofContainer::iterator iit = oldRecvDofs[it->first].begin(); 
+	     iit != oldRecvDofs[it->first].end(); ++iit)
+	  if (vertexDof[*iit]) {
+	    dofsToRecv.push_back(*iit);
+
+	    DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *iit);
+	    if (eraseIt != rankDOFs.end())
+	      rankDOFs.erase(eraseIt);    
+	  }
+
+
+	/*	switch (boundIt->rankObject.ithObjAtBoundary) {
 	case 0:
 	  if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
 	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
@@ -1225,12 +1261,20 @@ namespace AMDiS {
 
 	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
 	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *dofIt);
-	  if (eraseIt != rankDOFs.end())
-	    rankDOFs.erase(eraseIt);
+	  if (eraseIt != rankDOFs.end()) {
+	    if (mpiRank == 4)
+	      std::cout << "ERASE a: " << **eraseIt << std::endl;
+
+	    rankDOFs.erase(eraseIt);  
+	  }
+	  
 	  if (find(dofsToRecv.begin(), dofsToRecv.end(), *dofIt) == dofsToRecv.end())
-	    dofsToRecv.push_back(*dofIt);
+	    if (find(oldRecvDofs[it->first].begin(), oldRecvDofs[it->first].end(), *dofIt) != oldRecvDofs[it->first].end())
+	      dofsToRecv.push_back(*dofIt);
 	}
 
+	*/
+
 	dofs.clear();
 	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
 		       dofs);
@@ -1238,12 +1282,15 @@ namespace AMDiS {
 			 dofs);
 
 	for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
-	  TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), dofs[i]) != rankDOFs.end())
-	    ("Should never happen!\n");
+// 	  TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), dofs[i]) != rankDOFs.end())
+// 	    ("Should never happen!\n");
+	  TEST_EXIT_DBG(find(dofsToRecv.begin(), dofsToRecv.end(), dofs[i]) == dofsToRecv.end())
+	    ("Should not happen!\n");
+
 
 	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), dofs[i]);
 	  if (eraseIt != rankDOFs.end())
-	    rankDOFs.erase(eraseIt);
+	    rankDOFs.erase(eraseIt);	  
 
 	  dofsToRecv.push_back(dofs[i]);
 	}
@@ -1307,8 +1354,9 @@ namespace AMDiS {
       sendBuffers[i] = new int[nSendDofs];
       int c = 0;
       for (DofContainer::iterator dofIt = sendIt->second.begin();
-	   dofIt != sendIt->second.end(); ++dofIt)
+	   dofIt != sendIt->second.end(); ++dofIt) {
 	sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
+      }
 
       request[requestCounter++] = 
 	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
@@ -1336,7 +1384,7 @@ namespace AMDiS {
       for (DofContainer::iterator dofIt = recvIt->second.begin();
 	   dofIt != recvIt->second.end(); ++dofIt) {
 
-	rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];
+	rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];     
 	isRankDof[rankDofsNewLocalIndex[*dofIt]] = false;
 	j++;
       }
@@ -1455,12 +1503,14 @@ namespace AMDiS {
   void ParallelDomainBase::createDOFMemberInfo(DofToPartitions& partitionDofs,
 					       DofContainer& rankOwnedDofs,
 					       DofContainer& rankAllDofs,
-					       DofToRank& boundaryDofs)
+					       DofToRank& boundaryDofs,
+					       DofToBool& vertexDof)
   {
     partitionDofs.clear();
     rankOwnedDofs.clear();
     rankAllDofs.clear();
     boundaryDofs.clear();
+    vertexDof.clear();
 
     // === Determine to each dof the set of partitions the dof belongs to. ===
 
@@ -1474,6 +1524,11 @@ namespace AMDiS {
       do {
 	// Determine to each dof the partition(s) it corresponds to.
 	partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
+
+	if (elDofIt.getCurrentPos() == 0) 
+	  vertexDof[elDofIt.getDofPtr()] = true;
+	else
+	  vertexDof[elDofIt.getDofPtr()] = false;
       } while(elDofIt.next());
 
       elInfo = stack.traverseNext(elInfo);
@@ -1486,11 +1541,13 @@ namespace AMDiS {
     for (DofToPartitions::iterator it = partitionDofs.begin();
 	 it != partitionDofs.end(); ++it) {
 
+      bool isInRank = false;
+
       // iterate over all partition the current DOF is part of.
       for (std::set<int>::iterator itpart1 = it->second.begin();
 	   itpart1 != it->second.end(); ++itpart1) {
 
-	if (*itpart1 == mpiRank) {
+	if (*itpart1 == mpiRank) {	  
 	  rankAllDofs.push_back(it->first);
 
 	  if (it->second.size() == 1) {
@@ -1516,10 +1573,14 @@ namespace AMDiS {
 	    boundaryDofs[it->first] = highestRank;
 	  }
 
+	  isInRank = true;
 	  break;
 	}
 
       }
+
+      if (!isInRank) 
+	vertexDof.erase(it->first);
     }
 
     sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
@@ -1657,7 +1718,14 @@ namespace AMDiS {
       for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
 	   dofIt != it->second.end(); ++dofIt) {
 	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
-	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
+
+	if (!b) {
+	  std::cout << "ERROR (send): Cannot find DOF in mesh!" << std::endl;
+	  std::cout << "  mpi rank = " << mpiRank 
+		    << "  send to rank = " << it->first 
+		    << "  dof = " << *dofIt << " = " << **dofIt << std::endl;
+	  ERROR_EXIT("Should not happen!\n");
+	}
 
 	sendCoords[it->first].push_back(coords);
       }
@@ -1668,7 +1736,15 @@ namespace AMDiS {
       for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
 	   dofIt != it->second.end(); ++dofIt) {
 	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
-	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
+
+	if (!b) {
+	  std::cout << "ERROR (recv): Cannot find DOF in mesh!" << std::endl;
+	  std::cout << "  mpi rank = " << mpiRank 
+		    << "  send to rank = " << it->first 
+		    << "  dof = " << *dofIt << " = " << **dofIt << std::endl;
+	  ERROR_EXIT("Should not happen!\n");
+	}
+
 	recvCoords[it->first].push_back(coords);
       }
     }
diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h
index 7543ead1..6f6f3ef9 100644
--- a/AMDiS/src/ParallelDomainBase.h
+++ b/AMDiS/src/ParallelDomainBase.h
@@ -64,8 +64,11 @@ namespace AMDiS {
     /// Defines a mapping type from DOF indices to DOF indices.
     typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
 
+    /// Defines a mapping type from DOFs to boolean values.
+    typedef std::map<const DegreeOfFreedom*, bool> DofToBool;
+
     /// Defines a mapping type from DOF indices to boolean values.
-    typedef std::map<DegreeOfFreedom, bool> DofToBool;
+    typedef std::map<DegreeOfFreedom, bool> DofIndexToBool;
 
     /// Defines a mapping type from rank numbers to sets of coordinates.
     typedef std::map<int, std::vector<WorldVector<double> > > RankToCoords;
@@ -244,7 +247,8 @@ namespace AMDiS {
     void createDOFMemberInfo(DofToPartitions& partitionDofs,
 			     DofContainer& rankOwnedDofs,
 			     DofContainer& rankAllDofs,
-			     DofToRank& boundaryDofs);
+			     DofToRank& boundaryDofs,
+			     DofToBool& vertexDof);
 
     void setDofMatrix(DOFMatrix* mat, int dispMult = 1, 
 		      int dispAddRow = 0, int dispAddCol = 0);
@@ -403,7 +407,9 @@ namespace AMDiS {
      * owned by the rank. Otherwise, its an interior boundary DOF that is owned by
      * another rank.
      */
-    DofToBool isRankDof;
+    DofIndexToBool isRankDof;
+
+    DofToBool vertexDof;
 
     int rstart;
 
-- 
GitLab