Bugfix in domain parallelization.

AMDiS/src/AdaptInstationary.cc

@@ -108,7 +108,7 @@ namespace AMDiS {
 
       INFO(info_,6)("time = %e, try timestep = %e\n",
 		    adaptInfo->getTime(), adaptInfo->getTimestep());
-    
+
       problemIteration_->oneIteration(adaptInfo, NO_ADAPTION);
 
       adaptInfo->incTimestepIteration();

AMDiS/src/ElementFileWriter.cc

@@ -6,7 +6,7 @@
 
 namespace AMDiS {
 
-  ElementFileWriter::ElementFileWriter(const std::string& name_, 
+  ElementFileWriter::ElementFileWriter(std::string name_, 
 				       const FiniteElemSpace *feSpace_,
 				       std::map<int, double> &mapvec)
     : name(name_),
@@ -95,14 +95,14 @@ namespace AMDiS {
 
   void ElementFileWriter::writeFile(std::map<int, double> &vec,
 				    const FiniteElemSpace *feSpace,
-				    const std::string& filename)
+				    std::string filename)
   {
     ElementFileWriter efw("", feSpace, vec);
     efw.writeVtkValues(filename);
   }
   
 
-  void ElementFileWriter::writeTecPlotValues(const std::string &filename)
+  void ElementFileWriter::writeTecPlotValues(std::string filename)
   {
     FUNCNAME("ElementFileWriter::writeTecPlotValues()");
     std::ofstream fout(filename.c_str());
@@ -154,10 +154,9 @@ namespace AMDiS {
     
       // Write coordinates of all element vertices and element value.
       for (int i = 0; i <= dim; ++i) {
-      
-	for (int j = 0; j < dim; ++j) {
+      	for (int j = 0; j < dim; ++j)
 	  fout << elInfo->getCoord(i)[j] << " ";
-	}
+	
 	fout << val << "\n";
       }
     
@@ -182,7 +181,7 @@ namespace AMDiS {
     fout.close();
   }
 
-  void ElementFileWriter::writeMeshDatValues(const std::string &filename, double time)
+  void ElementFileWriter::writeMeshDatValues(std::string filename, double time)
   {
     FUNCNAME("ElementFileWriter::writeMeshDatValues()");
     std::ofstream fout(filename.c_str());
@@ -281,7 +280,7 @@ namespace AMDiS {
     fout.close();
   }
 
-  void ElementFileWriter::writeVtkValues(const std::string &filename)
+  void ElementFileWriter::writeVtkValues(std::string filename)
   {
     FUNCNAME("ElementFileWriter::writeVtkValues()");
     std::ofstream fout(filename.c_str());
@@ -289,6 +288,7 @@ namespace AMDiS {
     TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
 
     int dim = mesh->getDim();
+    int dimOfWorld = Global::getGeo(WORLD);
     int vertices = mesh->getGeo(VERTEX);
     int nElements = mesh->getNumberOfLeaves();
     double val;
@@ -312,12 +312,12 @@ namespace AMDiS {
     while (elInfo) {
       // Write coordinates of all element vertices.
       for (int i = 0; i <= dim; i++) {     
-	for (int j = 0; j < dim; j++) {
+	for (int j = 0; j < dimOfWorld; j++) 
 	  fout << elInfo->getCoord(i)[j] << " ";
-	}
-	for (int j = dim; j < 3; j++) {
+	
+	for (int j = dimOfWorld; j < 3; j++) 
 	  fout << "0.0";
-	}
+	
 	fout << "\n";
       }
    
@@ -329,9 +329,8 @@ namespace AMDiS {
     fout << "      <Cells>" << std::endl;
 
     fout << "        <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
-    for (int i = 0; i < nElements; i++) {
-      fout << " " << (i + 1) * vertices << std::endl;
-    }
+    for (int i = 0; i < nElements; i++) 
+      fout << " " << (i + 1) * vertices << std::endl;    
     fout << "        </DataArray>" << std::endl;
 
 

AMDiS/src/ElementFileWriter.h

@@ -21,7 +21,7 @@ namespace AMDiS {
   {
   public:
     /// Constructor.
-    ElementFileWriter(const std::string& name, 
+    ElementFileWriter(std::string name, 
 		      const FiniteElemSpace *feSpace,
 		      std::map<int, double> &vec);
 
@@ -34,17 +34,17 @@ namespace AMDiS {
     /// Simple writing procedure for one element map.
     static void writeFile(std::map<int, double> &vec,
 			  const FiniteElemSpace *feSpace,
-			  const std::string& filename);
+			  std::string filename);
 
   protected:
     /// Writes element data in tecplot format.
-    void writeTecPlotValues(const std::string &filename);
+    void writeTecPlotValues(std::string filename);
 
     /// Writes element data in AMDiS format (1 file !).
-    void writeMeshDatValues(const std::string &filename, double time);
+    void writeMeshDatValues(std::string filename, double time);
 
     /// Writes element data in VTK format.
-    void writeVtkValues(const std::string &filename);
+    void writeVtkValues(std::string filename);
 
   protected:
     /// Name.

AMDiS/src/InteriorBoundary.h

@@ -35,6 +35,9 @@ namespace AMDiS {
     /// The macro element to which the boundary element corresponds to.
     Element* el;
 
+    /// Index of the macro element.
+    int elIndex;
+
     /** \brief
      * Defines the geometrical object at the boundary. It must be "a part" of the
      * macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 (a face).

AMDiS/src/Lagrange.cc

@@ -48,12 +48,11 @@ namespace AMDiS {
 
   Lagrange::~Lagrange()
   {
-    for (int i = 0; i < static_cast<int>(bary->size()); i++) {
+    for (int i = 0; i < static_cast<int>(bary->size()); i++)
       if ((*bary)[i]) {
 	delete (*bary)[i];
 	(*bary)[i] = NULL;
       }
-    }
   }
 
   Lagrange* Lagrange::getLagrange(int dim, int degree)
@@ -70,13 +69,12 @@ namespace AMDiS {
 
   void Lagrange::clear()
   {    
-    std::list<Lagrange*>::iterator it;
-    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
+    for (std::list<Lagrange*>::iterator it = allBasFcts.begin(); 
+	 it != allBasFcts.end(); it++)
       if (*it) {
 	delete *it;
 	*it = NULL;
       }     
-    }
   }
 
   void Lagrange::setFunctionPointer()

AMDiS/src/ParallelDomainBase.cc

@@ -15,6 +15,7 @@
 #include "ElementDofIterator.h"
 #include "ProblemStatBase.h"
 #include "StandardProblemIteration.h"
+#include "ElementFileWriter.h"
 
 #include "petscksp.h"
 
@@ -22,7 +23,7 @@ namespace AMDiS {
 
   PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
   {    
-    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
+    //    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
       std::cout << "  Iteration " << iter << ": " << rnorm << std::endl;
 
     return 0;
@@ -66,6 +67,23 @@ namespace AMDiS {
     if (mpiSize <= 1)
       return;
 
+#if 0
+    if (mpiRank == 0) {
+      std::map<int, double> vec;
+      
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1,
+					   Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+      while (elInfo) {
+	vec[elInfo->getElement()->getIndex()] = static_cast<double>(elInfo->getElement()->getIndex());
+	elInfo = stack.traverseNext(elInfo);
+      }
+
+      ElementFileWriter::writeFile(vec, feSpace, "test.vtu");
+    }
+#endif
+
+
     // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
     // only for macro meshes, so it will not work yet if the mesh is already refined
     // in some way.
@@ -112,6 +130,41 @@ namespace AMDiS {
 
     updateDofAdmins();
 
+#if 0
+    if (mpiRank == 0) {
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+      while (elInfo) {
+	if (elInfo->getElement()->getIndex() == 4) {
+	  WorldVector<double> x;
+	  mesh->getDofIndexCoords(elInfo->getElement()->getDOF(0), feSpace, x);
+	  std::cout << "FOUND!" << std::endl;
+	  x.print();
+	}
+
+	elInfo = stack.traverseNext(elInfo);
+      }
+    }
+
+
+    if (mpiRank == 1) {
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+      while (elInfo) {
+	if (elInfo->getElement()->getIndex() == 3) {
+	  WorldVector<double> x;
+	  mesh->getDofIndexCoords(elInfo->getElement()->getDOF(2), feSpace, x);
+	  std::cout << "FOUND!" << std::endl;
+	  x.print();
+	}
+
+	elInfo = stack.traverseNext(elInfo);
+      }
+    }
+
+    exit(0);
+#endif
+
 
     // === Global refinements. ===
 
@@ -137,7 +190,7 @@ namespace AMDiS {
 
 
 #if (DEBUG != 0)
-    DbgTestCommonDofs();
+    DbgTestCommonDofs(true);
 #endif
 
 
@@ -293,7 +346,7 @@ namespace AMDiS {
     KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
     KSPGetPC(ksp, &pc);
     //    PCSetType(pc, PCNONE);
-    PCSetType(pc, PCJACOBI);
+    PCSetType(pc, PCILU);
     KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
     KSPSetType(ksp, KSPBCGS);
     //KSPSetType(ksp, KSPCG);
@@ -555,6 +608,24 @@ namespace AMDiS {
 	    bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end());
 	    bool ranksBoundary = isRankDOF1 || isRankDOF2;
 
+#if 0
+	    if (mpiRank == 3 && ranksBoundary && 
+		partitionVec[elInfo->getNeighbour(i)->getIndex()] == 2) {
+	      std::cout << "ADD MY BOUND " << element->getIndex() << "/" << i 
+			<< "  with  " 
+			<< elInfo->getNeighbour(i)->getIndex() << "/"
+			<< elInfo->getSideOfNeighbour(i) << std::endl;
+	    }
+
+	    if (mpiRank == 2 && !ranksBoundary && 
+		partitionVec[elInfo->getNeighbour(i)->getIndex()] == 3) {
+	      std::cout << "ADD OT BOUND " << element->getIndex() << "/" << i 
+			<< "  with  " 
+			<< elInfo->getNeighbour(i)->getIndex() << "/"
+			<< elInfo->getSideOfNeighbour(i) << std::endl;
+	    }
+#endif
+
 	    // === And add the part of the interior boundary. ===
 
 	    AtomicBoundary& bound = 
@@ -563,11 +634,19 @@ namespace AMDiS {
 	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));
 
 	    bound.rankObject.el = element;
+	    bound.rankObject.elIndex = element->getIndex();
 	    bound.rankObject.subObjAtBoundary = EDGE;
 	    bound.rankObject.ithObjAtBoundary = i;
-	    bound.neighbourObject.el = elInfo->getNeighbour(i);
+	    // Do not set a pointer to the element, because if will be deleted from
+	    // mesh after partitioning and the pointer would become invalid.
+	    bound.neighbourObject.el = NULL;
+	    bound.neighbourObject.elIndex = elInfo->getNeighbour(i)->getIndex();
 	    bound.neighbourObject.subObjAtBoundary = EDGE;
-	    bound.neighbourObject.ithObjAtBoundary = -1;
+	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);
+
+	    // i == 2  =>  getSideOfNeighbour(i) == 2
+	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
+	      ("Should not happen!\n");
  	  }
 	}
       }
@@ -589,23 +668,25 @@ namespace AMDiS {
     for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
 	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
       int nSendInt = rankIt->second.size();
-      int* buffer = new int[nSendInt];
-      for (int i = 0; i < nSendInt; i++)
-	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
+      int* buffer = new int[nSendInt * 2];
+      for (int i = 0; i < nSendInt; i++) {
+	buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
+	buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
+      }
       sendBuffers.push_back(buffer);
       
       request[requestCounter++] =
-	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
+	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
     }
 
     for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
 	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
       int nRecvInt = rankIt->second.size();
-      int *buffer = new int[nRecvInt];
+      int *buffer = new int[nRecvInt * 2];
       recvBuffers.push_back(buffer);
 
       request[requestCounter++] = 
-	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
+	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
     }
 
 
@@ -622,10 +703,12 @@ namespace AMDiS {
       
       for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
 	// If the expected object is not at place, search for it.
-	if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
+	if ((rankIt->second)[j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
+	    (rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
 	  int k = j + 1;
 	  for (; k < static_cast<int>(rankIt->second.size()); k++)
-	    if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
+	    if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] && 
+		(rankIt->second)[k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
 	      break;
 
 	  // The element must always be found, because the list is just in another order.
@@ -790,6 +873,11 @@ namespace AMDiS {
 	sendBuffers[i][c++] = *(dofIt->first);
 	sendBuffers[i][c++] = dofIt->second;
 
+#if 0
+	if (mpiRank == 3 && sendIt->first == 2)
+	  std::cout << "SEND DOF: " << dofIt->first << std::endl;
+#endif
+
 	sendDofs[sendIt->first].push_back(dofIt->first);
       }
 
@@ -852,6 +940,11 @@ namespace AMDiS {
 	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
 	    dofChanged[dofIt->first] = true;
 
+#if 0
+	    if (mpiRank == 2 && recvIt->first == 3)
+	      std::cout << "RECV DOF: " << dofIt->first << std::endl;
+#endif
+	
 	    recvDofs[recvIt->first].push_back(dofIt->first);
 	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
 	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
@@ -897,8 +990,15 @@ namespace AMDiS {
     }
 
     DofContainer rankAllDofs;
-    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt)
+    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt) {
+      /*      if (mpiRank == 1) {
+	std::cout << "COORDs of dof = " << **dofIt << std::endl;
+	WorldVector<double> x;
+	mesh->getDofIndexCoords(*dofIt, feSpace, x);
+	x.print();	
+	}*/
       rankAllDofs.push_back(*dofIt);
+    }
     sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
     DofContainer rankDOFs = rankAllDofs;
 
@@ -915,6 +1015,16 @@ namespace AMDiS {
       for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
 	   boundIt != it->second.end(); ++boundIt) {
 
+#if 0
+	if (mpiRank == 3 && it->first == 2) 
+	  std::cout << "GO ON MY BOUND: " << boundIt->rankObject.elIndex 
+		    << "/" << boundIt->rankObject.ithObjAtBoundary << " with "
+		    << boundIt->neighbourObject.elIndex << "/"
+		    << boundIt->neighbourObject.ithObjAtBoundary
+		    << std::endl;
+#endif
+
+
 	DofContainer dofs;
 	DofContainer &dofsToSend = sendDofs[it->first];
 
@@ -935,22 +1045,32 @@ 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 0
+	if (mpiRank == 3 && it->first == 2) {
+	  WorldVector<double> x;
+	  mesh->getDofIndexCoords(dofs[0], feSpace, x);
+	  x.print();
+	  mesh->getDofIndexCoords(dofs[1], feSpace, x);
+	  x.print();
 	}
+#endif
+	
+
+	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt)
+	  if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
+	    dofsToSend.push_back(*dofIt);	
 
 	dofs.clear();
 	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
 			 dofs);