Fixed several small problems in 3D refinement in parallel code.

AMDiS/libtool

@@ -44,7 +44,7 @@ available_tags=" CXX F77"
 
 # ### BEGIN LIBTOOL CONFIG
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -7263,7 +7263,7 @@ disable_libs=static
 # End:
 # ### BEGIN LIBTOOL TAG CONFIG: CXX
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -7568,7 +7568,7 @@ include_expsyms=""
 
 # ### BEGIN LIBTOOL TAG CONFIG: F77
 
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"

AMDiS/other/include/Makefile_AMDiS.mk

@@ -39,7 +39,7 @@ MPCCI_LIB = -L$(MPCCI_DIR)/lib/linux-x86-glibc22 -lmpcci
 PARMETIS_LIB = -L$(PARMETIS_DIR) -lparmetis -lmetis
 
 LIBS = $(AMDIS_LIB) $(PNG_LIB)
-LIBS += -lboost_iostreams -lboost_filesystem -lboost_system
+LIBS += -lboost_iostreams -lboost_filesystem -lboost_system -lboost_date_time
 
 ifeq ($(strip $(USE_UMFPACK)), 1)
 	LIBS += $(UMFPACK_LIB)

AMDiS/src/ElInfo3d.cc

+#include <typeinfo>
+
 #include "ElInfo3d.h"
 #include "BasisFunction.h"
 #include "Element.h"
@@ -478,6 +480,8 @@ namespace AMDiS {
   {
     FUNCNAME("ElInfo3d::fillElInfo()");
 
+    TEST_EXIT_DBG(elInfoOld)("Missing old elInfo!\n");
+
     int ochild = 0;             /* index of other child = 1-ichild */
     int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
     const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
@@ -497,7 +501,13 @@ namespace AMDiS {
     parent  = el_old;
     level = elInfoOld->level + 1;
     iChild = ichild;
-    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
+    int el_type_local = 0;
+    try {
+      el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
+    } catch (const std::bad_cast& e) {
+      ERROR_EXIT("ElInfo is not of type ElInfo3D but %s!\n", typeid(*elInfoOld).name());
+    }
+
     elType = (el_type_local + 1) % 3;
 
     TEST_EXIT_DBG(element)("missing child %d?\n", ichild);

AMDiS/src/Element.cc

@@ -11,6 +11,7 @@ namespace AMDiS {
 
   std::map<DegreeOfFreedom*, bool> Element::deletedDOFs;
 
+
   int Element::getRegion() const 
   {
     if (!elementData) 
@@ -25,6 +26,7 @@ namespace AMDiS {
     return -1;
   }
 
+
   void Element::setDOFPtrs()
   {
     FUNCNAME("Element::setDOFPtrs()");
@@ -34,6 +36,7 @@ namespace AMDiS {
     dof = mesh->createDOFPtrs();  
   }
 
+
   Element::Element(Mesh *aMesh)
   {
     mesh = aMesh;
@@ -51,6 +54,7 @@ namespace AMDiS {
     }
   }
 
+
   // call destructor through Mesh::freeElement !!!
   Element::~Element()
   {  
@@ -69,6 +73,7 @@ namespace AMDiS {
     }
   }
 
+
   bool Element::deleteElementData(int typeID)
   {
     FUNCNAME("Element::deleteElementData()");
@@ -87,6 +92,7 @@ namespace AMDiS {
     return false;
   }
 
+
   void Element::deleteElementDOFs()
   {
     int dim = mesh->getDim();
@@ -310,6 +316,7 @@ namespace AMDiS {
 #undef CHANGE_DOFS_1
 #undef CHANGE_DOFS_2
 
+
   /****************************************************************************/
   /* opp_vertex checks whether the face with vertices dof[0],..,dof[DIM-1] is */
   /* part of mel's boundary. returns the opposite vertex if true, -1 else     */
@@ -362,13 +369,16 @@ namespace AMDiS {
     }
   }
 
-  void Element::eraseNewCoord() {
+
+  void Element::eraseNewCoord() 
+  {
     if (newCoord != NULL) {
       delete newCoord;
       newCoord = NULL;
     }
   }
  
+
   void Element::serialize(std::ostream &out) 
   {
     // write children
@@ -445,6 +455,7 @@ namespace AMDiS {
     }
   }
 
+
   void Element::deserialize(std::istream &in)
   {
     FUNCNAME("Element::deserialize()");
@@ -555,6 +566,7 @@ namespace AMDiS {
     }
   }
 
+
   int Element::calcMemoryUsage()
   {
     int result = 0;
@@ -568,4 +580,55 @@ namespace AMDiS {
     return result;
   }
 
+
+  void writeDotFile(Element *el, std::string filename, int maxLevels)
+  {
+    std::vector<int> graphNodes;
+    std::vector<std::pair<int, int> > graphEdges;
+
+    std::vector<Element*> traverseElements, nextTraverseElements;
+    traverseElements.push_back(el);
+
+    int nLevels = 0;
+
+    while (!traverseElements.empty() && (maxLevels == -1 || nLevels < maxLevels)) {
+      nextTraverseElements.clear();
+
+      for (unsigned int i = 0; i < traverseElements.size(); i++) {
+	graphNodes.push_back(traverseElements[i]->getIndex());
+
+	if (!traverseElements[i]->isLeaf() && 
+	    (maxLevels == -1 || nLevels + 1 < maxLevels)) {
+	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
+					      traverseElements[i]->getChild(0)->getIndex()));
+	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
+					      traverseElements[i]->getChild(1)->getIndex()));
+	  nextTraverseElements.push_back(traverseElements[i]->getChild(0));
+	  nextTraverseElements.push_back(traverseElements[i]->getChild(1));
+	}
+      }
+
+      traverseElements = nextTraverseElements;
+      nLevels++;
+    }
+
+    std::ofstream file;
+    file.open(filename.c_str());
+
+    file << "digraph G\n";
+    file << "{\n";
+
+    for (unsigned int i = 0; i < graphNodes.size(); i++)
+      file << "   node" << graphNodes[i] 
+	   << " [ label = \"" << graphNodes[i] << "\"];\n";    
+
+    for (unsigned int i = 0; i < graphEdges.size(); i++) 
+      file << "   \"node" << graphEdges[i].first << "\" -> \"node"
+	   << graphEdges[i].second << "\";\n";
+
+    file << "}\n";
+
+    file.close();
+  }
+
 }

AMDiS/src/Element.h

@@ -583,6 +583,11 @@ namespace AMDiS {
 
     friend class Mesh;
   };
+
+
+  /// Writes the element hierarchie to a Graphviz dot-file. Using the dot-tool from
+  /// Graphvis, this dot-file can be converted to a ps-file. Useful for debugging!
+  void writeDotFile(Element *el, std::string filename, int maxLevels = -1);
 }
 
 #endif  // AMDIS_ELEMENT_H

AMDiS/src/FixVec.h

@@ -504,7 +504,7 @@ namespace AMDiS {
 
   /// returns the euclidian distance of a and b
   template<typename T, GeoIndex d>
-  double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);
+  double absteukl(const FixVec<T,d>& a, const FixVec<T,d>& b);
 
   /// FixVec operator for stream output
   template<typename T, GeoIndex d>
@@ -587,6 +587,15 @@ namespace AMDiS {
 
   template<typename T>
     const WorldMatrix<T>& operator+=(WorldMatrix<T>& m1, const WorldMatrix<T>& m2);
+
+
+  inline double norm(const WorldVector<double>& v)
+  {
+    double val = 0.0;
+    for (int i = 0; i < Global::getGeo(WORLD); i++)
+      val += v[i] * v[i];
+    return sqrt(val);
+  }
 }
 
 #include "FixVec.hh"

AMDiS/src/MeshStructure.cc

@@ -330,8 +330,16 @@ namespace AMDiS {
 	cont = nextElement();
       }
     }
-    
-    MSG("Mesh structure code: %s\n", oss.str().c_str());
+   
+    if (oss.str().length() < 255) {
+      MSG("Mesh structure code: %s\n", oss.str().c_str());
+    } else {
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << oss.str() << "\n";
+#else
+      std::cout << "                Mesh structure code: " << oss.str() << "\n";
+#endif
+    }
   }
 
 

AMDiS/src/parallel/InteriorBoundary.cc

@@ -18,6 +18,9 @@ namespace AMDiS {
       break;
 
     case 3:
+      TEST_EXIT_DBG(otherBound.elType == 0)
+	("Only 3D macro elements with level 0 are supported!\n");
+
       if (subObj == EDGE) {	
 	int el0_v0 = el->getVertexOfEdge(ithObj, 0);
 	int el0_v1 = el->getVertexOfEdge(ithObj, 1);
@@ -52,7 +55,7 @@ namespace AMDiS {
 	  otherMode = (localDofs0[0] != localDofs1[0]);
 	  
 	if (ithObj == 0)
-	  otherMode = localDofs0[1] != localDofs1[1];	
+	  otherMode = (localDofs0[1] != localDofs1[1]);
       }
       break;
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -105,10 +105,12 @@ namespace AMDiS {
       int writePartMesh = 1;
       GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
 
-      if (writePartMesh > 0)
+      if (writePartMesh > 0) {
+	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
 	writePartitioningMesh("part.vtu");
-      else 
+      } else {
 	MSG("Skip write part mesh!\n");
+      }
     }
 
     ParallelDebug::testAllElements(*this);
@@ -157,7 +159,6 @@ namespace AMDiS {
 
     createPeriodicMap();    
 
-
     // === Global refinements. ===
 
     int globalRefinement = 0;
@@ -172,7 +173,8 @@ namespace AMDiS {
 
       mesh->dofCompress();
       updateLocalGlobalNumbering();
-      
+
+     
       // === Update periodic mapping, if there are periodic boundaries. ===
       
       createPeriodicMap();
@@ -492,7 +494,7 @@ namespace AMDiS {
       RankToBoundMap allBound;
 
       for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
- 	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
  	  allBound[it.getRank()].push_back(*it);
 
       for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
@@ -500,7 +502,7 @@ namespace AMDiS {
 	  allBound[it.getRank()].push_back(*it);
 
       for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
- 	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
  	  allBound[it.getRank()].push_back(*it);	
 
 
@@ -579,7 +581,7 @@ namespace AMDiS {
 
 	MeshStructure elCode;	
 	elCode.init(boundIt->rankObj);
-	
+
 	if (elCode.getCode() != recvCodes[i].getCode()) {
 	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
 
@@ -1083,7 +1085,7 @@ namespace AMDiS {
 	      
 	      bound.neighObj.el = elIndexMap[it2->second.elIndex];
 	      bound.neighObj.elIndex = it2->second.elIndex;
-	      bound.neighObj.elType = -1;
+	      bound.neighObj.elType = elIndexTypeMap[it2->second.elIndex];
 	      bound.neighObj.subObj = geoIndex;
 	      bound.neighObj.ithObj = it2->second.ithObject;
 	      
@@ -1435,7 +1437,7 @@ namespace AMDiS {
       it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
 
       for (unsigned int i = 0; i < dofs.size(); i++)
-	sendDofs[it.getRank()].push_back(dofs[i]);                  
+	sendDofs[it.getRank()].push_back(dofs[i]);          
     }
 
 
@@ -1448,7 +1450,7 @@ namespace AMDiS {
 	DofContainer::iterator eraseIt = 
 	  find(rankDofs.begin(), rankDofs.end(), dofs[i]);
 	if (eraseIt != rankDofs.end()) 
-	  rankDofs.erase(eraseIt);	 
+	  rankDofs.erase(eraseIt);
 
 	recvDofs[it.getRank()].push_back(dofs[i]);
       }
@@ -1796,7 +1798,7 @@ namespace AMDiS {
 	int *sendbuf = new int[2];
 	sendbuf[0] = periodicDof[type0][it->first];
 	sendbuf[1] = periodicDof[type1][it->first];
-	
+
 	request[requestCounter++] = 
 	  mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0);
 	request[requestCounter++] = 

AMDiS/src/parallel/ParallelDebug.cc

@@ -139,13 +139,13 @@ namespace AMDiS {
     typedef std::map<int, DofContainer> RankToDofContainer;
 
     // Maps to each neighbour rank an array of WorldVectors. This array contains the 
-    // coordinates of all dofs this rank shares on the interior boundary with the 
+    // coordinates of all DOFs this rank shares on the interior boundary with the 
     // neighbour rank. A rank sends the coordinates to another rank, if it owns the
-    // boundarys dof.
+    // boundarys DOFs.
     RankToCoords sendCoords;
 
-    // A rank receives all boundary dofs that are at its interior boundaries but are
-    // not owned by the rank. This map stores for each rank the coordinates of dofs
+    // A rank receives all boundary DOFs that are at its interior boundaries but are
+    // not owned by the rank. This map stores for each rank the coordinates of DOFs
     // this rank expectes to receive from.
     RankToCoords recvCoords;
 
@@ -180,14 +180,16 @@ namespace AMDiS {
       if (i == pdb.mpiRank)
 	continue;
 
-      request[requestCounter++] = pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
+      request[requestCounter++] = 
+	pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
     }   
 
     for (int i = 0; i < pdb.mpiSize; i++) {
       if (i == pdb.mpiRank)
 	continue;
 
-      request[requestCounter++] = pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
+      request[requestCounter++] = 
+	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
     }
 
     MPI::Request::Waitall(requestCounter, request);
@@ -217,26 +219,24 @@ namespace AMDiS {
     stdMpi.recv(recvCoords);   
     stdMpi.startCommunication<double>(MPI_DOUBLE);
 
-    double eps = 1e-13;
     int dimOfWorld = Global::getGeo(WORLD);
 
     // === Compare the received with the expected coordinates. ===
 
     for (RankToCoords::iterator it = stdMpi.getRecvData().begin(); 
 	 it != stdMpi.getRecvData().end(); ++it) {
-      for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
-	for (int j = 0; j < dimOfWorld; j++) {
-	  bool c = fabs((it->second)[i][j] - recvCoords[it->first][i][j]) < eps;
-	  
-	  // === Print error message if the coordinates are not the same. ===
-
-	  if (printCoords && !c) {
-	    MSG("[DBG] i = %d\n", i);
+      for (unsigned int i = 0; i < it->second.size(); i++) {
+	WorldVector<double> tmp = (it->second)[i];
+	tmp -=  recvCoords[it->first][i];
 
+	if (norm(tmp) > 1e-13) {
+	  // === Print error message if the coordinates are not the same. ===
+	  if (printCoords) {
+	    MSG("[DBG] i = %d\n", i);	  
 	    std::stringstream oss;
 	    oss.precision(5);
 	    oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
-		      << " expect coords (";
+		<< " expect coords (";
 	    for (int k = 0; k < dimOfWorld; k++) {
 	      oss << recvCoords[it->first][i][k];
 	      if (k + 1 < dimOfWorld)
@@ -250,16 +250,13 @@ namespace AMDiS {
 	    }
 	    oss << ")";
 	    MSG("%s\n", oss.str().c_str());
-
+	    
 	    debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
 	  }
-
-	  if (!c) {
-	    ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
-	    foundError = 1;
-	  }	 
-	}
-      } 
+	  ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
+	  foundError = 1;
+	}	 
+      }
     }
     mpi::globalAdd(foundError);
     TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
@@ -459,6 +456,14 @@ namespace AMDiS {
       MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
 	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
     }
+
+    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
+      MSG("Periodic boundary with rank %d: \n", it.getRank());
+      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
+	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
+      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
+	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
+    }    
   }