Fixed serialization problem for parallel computations, work on periodic bc in parallel code.

AMDiS/configure

@@ -2074,7 +2074,7 @@ fi;
 MPI_DIR=$mpidir
 if test $mpidir != no ; then
-  CXX=$mpidir/bin/mpiCC
+  CXX=$mpidir/bin/mpicxx
   CC=$mpidir/bin/mpicc

AMDiS/configure.ac

@@ -50,7 +50,7 @@ fi
 AC_ARG_WITH(mpi, [  --with-mpi=MPI_DIR], mpidir=$withval, mpidir=no)
 AC_SUBST(MPI_DIR, $mpidir)
 if test $mpidir != no ; then
-  AC_SUBST(CXX, $mpidir/bin/mpiCC)
+  AC_SUBST(CXX, $mpidir/bin/mpicxx)
   AC_SUBST(CC, $mpidir/bin/mpicc)
 fi
 

AMDiS/libtool

@@ -44,7 +44,7 @@ available_tags=" CXX F77"
 
 # ### BEGIN LIBTOOL CONFIG
 
-# Libtool was configured on host p2q084:
+# Libtool was configured on host p1q020:
 
 # 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 p2q084:
+# Libtool was configured on host p1q020:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -6804,7 +6804,7 @@ LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
 LTCFLAGS="-g -O2"
 
 # A language-specific compiler.
-CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpiCC"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicxx"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -7065,7 +7065,7 @@ include_expsyms=""
 
 # ### BEGIN LIBTOOL TAG CONFIG: F77
 
-# Libtool was configured on host p2q084:
+# Libtool was configured on host p1q020:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"

AMDiS/src/MeshStructure.cc

@@ -63,6 +63,8 @@ namespace AMDiS {
 
     Element *el = bound.el;
 
+    TEST_EXIT_DBG(el)("No element!\n");
+
     clear();
 
     int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);

AMDiS/src/Parameters.cc

@@ -33,7 +33,7 @@ namespace AMDiS {
 
     if (Parameters::singlett->paramInfo) {
       if (Parameters::singlett->paramInfo > 1)
-	info = max(info, Parameters::singlett->paramInfo-1);
+	info = max(info, Parameters::singlett->paramInfo - 1);
     } else {
       info = 0;
     }
@@ -47,7 +47,8 @@ namespace AMDiS {
       }
     }
 
-    i = std::find(Parameters::singlett->allParam.begin(),Parameters::singlett->allParam.end(),tParam);
+    i = std::find(Parameters::singlett->allParam.begin(),
+		  Parameters::singlett->allParam.end(), tParam);
 
     if (i == Parameters::singlett->allParam.end()) {
       if (funcName != funcName2) {
@@ -186,6 +187,7 @@ namespace AMDiS {
     return count;
   }
 
+
   int Parameters::getGlobalParameter(int flag, const std::string& key, 
 				     std::string* param)
   {
@@ -197,6 +199,7 @@ namespace AMDiS {
     return result;
   }
 
+
   void Parameters::read(const std::string& aFilename, 
 			const std::string& maxKey)
   {
@@ -238,6 +241,7 @@ namespace AMDiS {
     singlett->inputFile.close();
   }
 
+
   const std::string& Parameters::getActFile(const std::string& aFilename)
   {
     FUNCNAME("Parameters::getActFile()");
@@ -255,6 +259,7 @@ namespace AMDiS {
     }
   }
 
+
   void Parameters::swap(int i, int j)
   {
     param tmp(allParam[i]);
@@ -262,6 +267,7 @@ namespace AMDiS {
     allParam[j] = tmp;
   }
 
+
   void Parameters::qsort(int left, int right)
   {
     if (left >= right) 

AMDiS/src/Parameters.h

@@ -171,9 +171,9 @@ namespace AMDiS {
      * specified in a parameter file.
      * \see GET_PARAMETER
      */
-    static int getGlobalParameter(int                 flag,
-				  const std::string&  key,
-				  const char         *format, 
+    static int getGlobalParameter(int flag,
+				  const std::string& key,
+				  const char *format, 
 				  ...);  
 
     /// Like getGlobalParameter(flag, key, "%s", param->c_str()).
@@ -245,7 +245,7 @@ namespace AMDiS {
 	    const std::string& nparameters, 
 	    const std::string& nfilename,
 	    const std::string& nfuncName, 
-	    int                line) 
+	    int line) 
 	: key(nkey),
 	  parameters(nparameters),
 	  filename(nfilename),
@@ -327,15 +327,14 @@ namespace AMDiS {
   private:
     std::string cppFlags;
     std::ifstream inputFile;
-    std::list< std::string> filenames;
+    std::list<std::string> filenames;
     size_t maxFiles;
     std::string filename;
     int paramInfo;
     int msgInfo;
     int msgWait;
     std::vector<param> allParam;
- 
-  };
+   };
 
 
   /** \brief

AMDiS/src/ProblemVec.cc

@@ -184,6 +184,8 @@ namespace AMDiS {
 
 	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
 #endif
+
+	deserialized = true;
       } else {
 	int globalRefinements = 0;
 
@@ -476,12 +478,8 @@ namespace AMDiS {
 
     GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
 
-    // Serialization is not allowed to be done by the problem, if its part of a parallel
-    // problem definition. Than, the parallel problem object must be serialized.
-#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
     if (writeSerialization)
       fileWriters.push_back(new Serializer<ProblemVec>(this));
-#endif
   }
 
 

AMDiS/src/ProblemVec.h

@@ -63,6 +63,7 @@ namespace AMDiS {
 	systemMatrix(NULL),
 	useGetBound(true),
 	info(10),
+	deserialized(false),
 	computeExactError(false),
 	boundaryConditionSet(false),
 	writeAsmInfo(false)
@@ -411,6 +412,12 @@ namespace AMDiS {
       return info;
     }
 
+    /// Returns \ref deserialized;
+    bool isDeserialized()
+    {
+      return deserialized;
+    }
+
     /** \} */
 
     /** \name setting methods
@@ -601,6 +608,9 @@ namespace AMDiS {
     /// Info level.
     int info;
 
+    /// If true, the stationary problem was deserialized from some serialization file.
+    bool deserialized;
+
     /** \brief
      * This vectors stores pointers to functions defining the exact solution of
      * the problem. This may be used to compute the real error of the computed

AMDiS/src/Serializer.h

@@ -46,13 +46,30 @@ namespace AMDiS {
       GET_PARAMETER(0, problem->getName() + "->output->serialization filename", &name);
       GET_PARAMETER(0, problem->getName() + "->output->write every i-th timestep", 
 		    "%d", &tsModulo);
-      TEST_EXIT(name != "")("no filename\n");
+      TEST_EXIT(name != "")("No filename!\n");
 
 #if HAVE_PARALLEL_DOMAIN_AMDIS
       name += ".p" + boost::lexical_cast<std::string>(MPI::COMM_WORLD.Get_rank());
 #endif
     }
 
+
+    Serializer(ProblemType *prob, std::string filename, int writeEveryIth)
+      : name(filename),
+	problem(prob),
+	tsModulo(writeEveryIth),
+	timestepNumber(-1)
+    {
+      FUNCNAME("Serializer::Serializer()");
+
+      TEST_EXIT(name != "")("No filename!\n");
+
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+      name += ".p" + boost::lexical_cast<std::string>(MPI::COMM_WORLD.Get_rank());
+#endif      
+    }
+
+
     virtual ~Serializer() {}
 
     virtual void writeFiles(AdaptInfo *adaptInfo, 

AMDiS/src/parallel/InteriorBoundary.cc

@@ -120,6 +120,8 @@ namespace AMDiS {
   void InteriorBoundary::deserialize(std::istream &in, 
 				     std::map<int, Element*> &elIndexMap)
   {
+    FUNCNAME("InteriorBoundary::deserialize()");
+
     int mSize = 0;
     SerUtil::deserialize(in, mSize);
     for (int i = 0; i < mSize; i++) {
@@ -146,6 +148,10 @@ namespace AMDiS {
 	SerUtil::deserialize(in, bound.neighObj.reverseMode);
 	deserializeExcludeList(in, bound.neighObj.excludedSubstructures);
 
+	TEST_EXIT_DBG(elIndexMap.count(bound.rankObj.elIndex) == 1)
+	  ("Cannot find element with index %d for deserialization!\n", 
+	   bound.rankObj.elIndex);
+
 	bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
 	bound.neighObj.el = NULL;
       }

AMDiS/src/parallel/ParallelDomainBase.cc

@@ -2,6 +2,8 @@
 #include <iostream>
 #include <fstream>
 #include <boost/lexical_cast.hpp>
+#include <boost/filesystem.hpp>
+
 #include "parallel/ParallelDomainBase.h"
 #include "parallel/ParallelDomainDbg.h"
 #include "parallel/StdMpi.h"
@@ -29,6 +31,7 @@
 namespace AMDiS {
 
   using boost::lexical_cast;
+  using namespace boost::filesystem;
 
   inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
   {
@@ -49,6 +52,7 @@ namespace AMDiS {
       nOverallDofs(0),
       rstart(0),
       deserialized(false),
+      writeSerializationFile(false),
       lastMeshChangeIndex(0),
       macroElementStructureConsisten(false)
   {
@@ -253,23 +257,49 @@ namespace AMDiS {
 
     // Create parallel serialization file writer, if needed.
     int writeSerialization = 0;
-    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", &writeSerialization);
-    if (writeSerialization)
-      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this));
+    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
+		  &writeSerialization);
+    if (writeSerialization && !writeSerializationFile) {
+      std::string f = "";
+      GET_PARAMETER(0, probVec->getName() + "->output->serialization filename", &f);
+      path myPath(f);
+      std::string meshFilename = 
+	myPath.parent_path().directory_string() + "/meshDistributor.ser";
+
+
+      int tsModulo = -1;
+      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
+		    "%d", &tsModulo);
+      
+      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, meshFilename, tsModulo));
+      writeSerializationFile = true;
+    }    
 
     int readSerialization = 0;
-    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", &readSerialization);
+    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
+		  &readSerialization);
     if (readSerialization) {
-      ERROR_EXIT("Must be reimplemented!\n");
-#if 0      
       std::string filename = "";
       GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
       filename += ".p" + lexical_cast<std::string>(mpiRank);
-      MSG("Start serialization with %s\n", filename.c_str());
+      MSG("Start deserialization with %s\n", filename.c_str());
       std::ifstream in(filename.c_str());
-      deserialize(in);
+      probVec->deserialize(in);
       in.close();
-#endif
+      MSG("Deserialization from file: %s\n", filename.c_str());
+
+      if (!deserialized) {
+	GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
+	path myPath(filename);
+	std::string meshFilename = 
+	  myPath.parent_path().directory_string() + "/meshDistributor.ser.p" + 
+	  lexical_cast<std::string>(mpiRank);
+	std::ifstream in(meshFilename.c_str());
+	deserialize(in);
+	in.close();
+	MSG("Deserializtion of mesh distributor from file: %s\n", meshFilename.c_str());
+	deserialized = true;
+      }
     }
 
     probStat.push_back(probVec);
@@ -475,6 +505,10 @@ namespace AMDiS {
 
     mesh->dofCompress();
     updateLocalGlobalNumbering();
+
+    // === Update periodic mapping, if there are periodic boundaries. ===
+
+    createPeriodicMap();
   }
 
   
@@ -990,10 +1024,14 @@ namespace AMDiS {
 	    // === to \ref otherIntBoundary. It dependes on which rank is respon-   ===
 	    // === sible for this boundary.                                         ===
 
-	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
+	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {
 	      // We have found an element that is at an interior, periodic boundary.
 	      AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
 	      b = bound;
+	      if (mpiRank > otherElementRank)
+		b.rankObj.setReverseMode(b.neighObj, feSpace);
+	      else
+		b.neighObj.setReverseMode(b.rankObj, feSpace);
 	    } else {
 	      // We have found an element that is at an interior, non-periodic boundary.
 	      if (mpiRank > otherElementRank) {
@@ -1927,7 +1965,7 @@ namespace AMDiS {
 
       // Send the global indices to the rank on the other side.
       stdMpi.getSendData(it->first).reserve(dofs.size());
-      for (int i = 0; i < static_cast<int>(dofs.size()); i++)
+      for (unsigned int i = 0; i < dofs.size(); i++)
 	stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]);
 
       // Receive from this rank the same number of dofs.
@@ -1955,14 +1993,14 @@ namespace AMDiS {
 	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, tmpdofs);
 	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, tmpdofs);
 
-	for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
-	  dofs.push_back(tmpdofs[j]);	
+	for (unsigned int i = 0; i < tmpdofs.size(); i++)
+  	  dofs.push_back(tmpdofs[i]);	
       }
 
       // Added the received dofs to the mapping.
-      for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
-	int globalDofIndex = mapLocalGlobalDofs[*(dofs[j])];
-	periodicDof[globalDofIndex].insert(stdMpi.getRecvData(it->first)[j]);
+      for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
+	int globalDofIndex = mapLocalGlobalDofs[*(dofs[i])];
+	periodicDof[globalDofIndex].insert(stdMpi.getRecvData(it->first)[i]);
 	dofFromRank[globalDofIndex].insert(it->first);
       }
     }
@@ -2028,9 +2066,11 @@ namespace AMDiS {
     SerUtil::serialize(out, partitionVec);
     SerUtil::serialize(out, oldPartitionVec);
     SerUtil::serialize(out, nRankDofs);
+    SerUtil::serialize(out, nOverallDofs);
     
     myIntBoundary.serialize(out);
     otherIntBoundary.serialize(out);
+    periodicBoundary.serialize(out);
 
     serialize(out, sendDofs);
     serialize(out, recvDofs);
@@ -2043,6 +2083,7 @@ namespace AMDiS {
     serialize(out, periodicDof);
 
     SerUtil::serialize(out, rstart);
+    SerUtil::serialize(out, macroElementStructureConsisten);
   }
 
 
@@ -2053,6 +2094,7 @@ namespace AMDiS {
     SerUtil::deserialize(in, partitionVec);
     SerUtil::deserialize(in, oldPartitionVec);
     SerUtil::deserialize(in, nRankDofs);    
+    SerUtil::deserialize(in, nOverallDofs);
 
     // Create two maps: one from from element indices to the corresponding element 
     // pointers, and one map from Dof indices to the corresponding dof pointers.
@@ -2060,20 +2102,24 @@ namespace AMDiS {
     std::map<int, const DegreeOfFreedom*> dofMap;
     ElementDofIterator elDofIter(feSpace);
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
     while (elInfo) {
       Element *el = elInfo->getElement();
       elIndexMap[el->getIndex()] = el;
-      elInfo = stack.traverseNext(elInfo);
 
-      elDofIter.reset(el);
-      do {
-	dofMap[elDofIter.getDof()] = elDofIter.getDofPtr();
-      } while (elDofIter.next());      
+      if (el->isLeaf()) {
+	elDofIter.reset(el);
+	do {
+	  dofMap[elDofIter.getDof()] = elDofIter.getDofPtr();
+	} while (elDofIter.next());      
+      }
+
+      elInfo = stack.traverseNext(elInfo);
     }
 
     myIntBoundary.deserialize(in, elIndexMap);
     otherIntBoundary.deserialize(in, elIndexMap);
+    periodicBoundary.deserialize(in, elIndexMap);
 
     deserialize(in, sendDofs, dofMap);
     deserialize(in, recvDofs, dofMap);
@@ -2085,6 +2131,7 @@ namespace AMDiS {
     deserialize(in, vertexDof, dofMap);
 
     SerUtil::deserialize(in, rstart);
+    SerUtil::deserialize(in, macroElementStructureConsisten);
 
     deserialized = true;
   }

AMDiS/src/parallel/ParallelDomainBase.h

@@ -372,12 +372,6 @@ namespace AMDiS {
     }
 		        
   protected:
-    ///
-    ProblemIterationInterface *iterationIF;
-
-    ///
-    ProblemTimeInterface *timeIF;
-
     ///
     std::vector<ProblemVec*> probStat;
 
@@ -511,6 +505,9 @@ namespace AMDiS {
      */
     bool deserialized;
 
+    /// Denotes whether there exists a filewriter for this object.
+    bool writeSerializationFile;
+
     /** \brief
      * Stores the mesh change index. This is used to recognize changes in the mesh 
      * structure (e.g. through refinement or coarsening managers).