correction of MeshManipulation

AMDiS/src/AbstractFunction.h

@@ -180,7 +180,7 @@ namespace AMDiS {
   #define ABSTRACT_FUNCTION_MACRO(z, n, _) \
    template< typename ReturnType, \
       BOOST_PP_ENUM_PARAMS_Z(z, n, typename T) \
-   > class AbstractFunction ## n { \
+   > struct AbstractFunction ## n { \
     AbstractFunction ## n (int degree = 0) : \
       degree_(degree) \
     {} \

AMDiS/src/DOFVector.hh

@@ -551,7 +551,7 @@ namespace AMDiS {
     TOut tmp; nullify(tmp);
       for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
 	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
- 	tmp += fastQuad->getWeight(iq) * (*fct)(coords);
+ 	tmp += (*fct)(coords) * fastQuad->getWeight(iq);
       }
 
       value += tmp * elInfo->getDet();

AMDiS/src/FixVec.h

@@ -468,6 +468,9 @@ namespace AMDiS {
       TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
       this->set(ini);
     }
+    
+    // import assignment operator from base class
+    using Matrix<T>::operator=;
   
     /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
     bool isDiagMatrix() const;

AMDiS/src/Global.h

@@ -197,8 +197,13 @@ namespace AMDiS {
 
   template<typename T> inline void nullify(WorldVector<T> &a)
   {
-    T null;
-    nullify(null);
+    T null; nullify(null);
+    a = null;
+  }
+
+  template<typename T> inline void nullify(WorldMatrix<T> &a)
+  {
+    T null; nullify(null);
     a = null;
   }
 

AMDiS/src/MatrixVector.h

@@ -296,6 +296,18 @@ namespace AMDiS {
       return this->valArray + cols * i;
     }
 
+    /// Access to i-th matrix row.
+    inline T& operator()(int i, int j) 
+    {
+      return this->operator[](i)[j];
+    }
+
+    /// Access to i-th matrix row for constant matrices.
+    inline const T& operator()(int i, int j) const 
+    {
+      return this->operator[](i)[j];
+    }
+
     /// Returns \ref rows.
     inline int getNumRows() const 
     { 

AMDiS/src/ProblemInstat.h

@@ -122,6 +122,10 @@ namespace AMDiS {
       return &invTau;
     }
 
+    virtual void serialize(std::ostream &out) {}
+
+    virtual void deserialize(std::istream &in) {}
+
   protected:
     /// Name of the problem.
     std::string name;
@@ -195,10 +199,6 @@ namespace AMDiS {
     /// Used by \ref problemInitial
     virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
 
-    virtual void serialize(std::ostream &out) {}
-
-    virtual void deserialize(std::istream &in) {}
-
   protected:
     /// Space problem solved in each timestep.
     ProblemStatSeq* problemStat;

AMDiS/src/Timer.cc

@@ -22,29 +22,34 @@
 #include <mpi.h>
 #endif
 
+using namespace boost::posix_time;
+
 namespace AMDiS {
-	Timer::Timer():
-		first_seq(clock())
+  
+  Timer::Timer() : 
+     first_seq(microsec_clock::local_time())
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-		,first_mpi(MPI::Wtime())
+    ,first_mpi(MPI::Wtime())
 #endif
-	{}
+  {}
 
-	void Timer::reset()
-	{
+  void Timer::reset()
+  {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-		first_mpi = MPI::Wtime();
+    first_mpi = MPI::Wtime();
 #else
-		first_seq = clock();
+    first_seq = microsec_clock::local_time();
 #endif
-	}
+  }
 
-	double Timer::elapsed()
-	{
+  double Timer::elapsed()
+  {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-		return MPI::Wtime() - first_mpi;
+    return MPI::Wtime() - first_mpi;
 #else
-		return ((double)(clock() - first_seq))/((double)CLOCKS_PER_SEC);
+    time_duration td = microsec_clock::local_time()-first_seq;
+    return static_cast<double>(td.total_milliseconds())*1.e-3;
 #endif
-	}
-}
+  }
+  
+} // end namespace AMDiS

AMDiS/src/Timer.h

@@ -20,26 +20,30 @@
 #ifndef AMDIS_TIMER_H
 #define AMDIS_TIMER_H
 
-#include <time.h>
-namespace AMDiS {
-	/// Helper class to destinct between different time measurement methods
-	class Timer {
-		/// begin value for sequentiell measurement
-		clock_t first_seq;
+#include "boost/date_time/posix_time/posix_time.hpp"
 
-		/// begin value for parallel measurement
-		double first_mpi;
+namespace AMDiS {
+  
+  /// Helper class to distinguish between different time measurement methods
+  class Timer
+  {
+  private:
+    /// begin value for sequentiell measurement
+    boost::posix_time::ptime first_seq;
 
-		public:
+    /// begin value for parallel measurement
+    double first_mpi;
 
-		///initializes the timer with current time
-		Timer();
+  public:
+    /// initializes the timer with current time
+    Timer();
 
-		/// resets the timer to current time
-		void reset();
+    /// resets the timer to current time
+    void reset();
 
-		/// returns the elapsed time (from construction or last reset) to now in seconds
-		double elapsed();
-	};
+    /// returns the elapsed time (from construction or last reset) to now in seconds
+    double elapsed();
+  };
 }
-#endif
+
+#endif // AMDIS_TIMER_H

AMDiS/src/parallel/MeshManipulation.cc

@@ -318,10 +318,14 @@ namespace AMDiS { namespace Parallel {
 #if (DEBUG != 0)
       ElInfo *elInfo = 
 	stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
+	
+      TEST_EXIT(elInfo->getElement() == boundEl.el)
+	("This should not happen!\n");
+#else
+      // stack must be initialized before passed to fitElementToMeshCode()
+      stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
 #endif
 
-      TEST_EXIT_DBG(elInfo->getElement() == boundEl.el)
-	("This should not happen!\n");
 
       pcode = &code;
       pstack = &stack;

AMDiS/src/parallel/MeshPartitioner.cc

@@ -68,9 +68,10 @@ namespace AMDiS { namespace Parallel {
 		      arhMetaFile);
       bool partitioningArhBased = (arhMetaFile != "");
       if (partitioningArhBased) {
+	MSG("Read Meta-Arh partitioning file: %s\n", arhMetaFile.c_str());
 	int nProc = ArhReader::readMetaData(arhMetaFile, mapElInRank, arhElCodeSize);
 	if (nProc != mpiSize)
-	  partitioningArhBased = false;
+	  useInitialPartitioning = false;
 	else
 	  useInitialPartitioning = true;
       }

AMDiS/src/parallel/MpiHelper.h

@@ -39,7 +39,7 @@ namespace AMDiS { namespace Parallel {
     {
       globalAdd(MPI::COMM_WORLD, value);
     }
-
+    
     void globalAdd(MPI::Intracomm &mpiComm, int &value);
 
     inline void globalAdd(int &value)
@@ -53,6 +53,14 @@ namespace AMDiS { namespace Parallel {
       WARNING("Unknown type for globalAdd. Can not sum up the values of all processors!\n");
     }
 
+    template<typename VectorType>
+    inline typename boost::enable_if< mtl::traits::is_vector<VectorType>, void >
+    globalAdd(VectorType &value)
+    {
+      for (size_t i = 0; i < num_rows(value); i++)
+	globalAdd(value[i]);
+    }
+
     void globalMin(double &value);
 
     void globalMin(int &value);

AMDiS/src/parallel/ParallelProblemStat.cc

@@ -179,6 +179,7 @@ namespace AMDiS { namespace Parallel {
       
     creator = new PetscSolverNavierStokes::Creator;
     CreatorMap< LinearSolver >::addCreator("p_petsc_petsc-navierstokes", creator);
+    CreatorMap< LinearSolver >::addCreator("petsc-navierstokes", creator);
     
 #elif defined  HAVE_BDDC_ML    
     creator = new BddcMlSolver::Creator;

AMDiS/src/solver/MatrixStreams.h

@@ -156,7 +156,7 @@ namespace AMDiS {
     bool initMatrix = false;
     if (mat == PETSC_NULL) {
       std::vector<PetscInt> nnz(rhs.getSize());
-      for (size_t k = 0; k < rhs.getSize(); k++)
+      for (size_t k = 0; k < static_cast<size_t>(rhs.getSize()); k++)
 	nnz[k] = rhs.getBaseMatrix().nnz_local(k);	  
       
       MatCreateSeqAIJ(PETSC_COMM_SELF, rhs.getSize(), rhs.getSize(), 0, &(nnz[0]), &mat);
@@ -184,13 +184,12 @@ namespace AMDiS {
   inline void operator<<(PetscMatrix& mat, const SolverMatrix<Matrix<DOFMatrix*> >& Asolver)
   {
     const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
-    int ns = A.getSize();
     
     std::vector<PetscInt> nnz;      
     mat.nestMat.resize(A.getNumRows() * A.getNumCols());
     
-    for (size_t i = 0; i < A.getNumRows(); i++) {
-      for (size_t j = 0; j < A.getNumCols(); j++) {
+    for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
+      for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
 	size_t idx = i * A.getNumCols() + j;
 	
 	if (A[i][j] == NULL 
@@ -244,14 +243,14 @@ namespace AMDiS {
     VecGetType(petscVec, &vecType);
     
     if (strcmp(vecType, VECNEST) == 0) {
-      for (size_t i = 0; i < vec.getSize(); i++) {  
+      for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {  
 	Vec v;
 	VecNestGetSubVec(petscVec, i, &v);
 	v << *(vec.getDOFVector(i));
       }
     } else {
       PetscInt index = 0;
-      for (size_t i = 0; i < vec.getSize(); i++) {
+      for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
 	DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
 	for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
 	  double value = *dofIt;
@@ -269,14 +268,14 @@ namespace AMDiS {
     VecGetType(petscVec, &vecType);
     
     if (strcmp(vecType, VECNEST) == 0) {
-      for (size_t i = 0; i < vec.getSize(); i++) {  
+      for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {  
 	Vec v;
 	VecNestGetSubVec(petscVec, i, &v);
 	v >> *(vec.getDOFVector(i));
       }
     } else {
-      int n = 0, ni;
-      for (size_t i = 0; i < vec.getSize(); i++)
+      PetscInt n = 0;
+      for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++)
 	n += vec.getDOFVector(i)->getUsedSize();
       
       PetscInt N = 0;
@@ -284,7 +283,7 @@ namespace AMDiS {
       assert(n == N);
       
       PetscInt index = 0;
-      for (size_t i = 0; i < vec.getSize(); i++) {
+      for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
 	DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
 	for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
 	  double value = 0.0;
@@ -317,7 +316,7 @@ namespace AMDiS {
       // create nested vector from vector of block vectors
       VecCreateNest(PETSC_COMM_SELF, nComponents, PETSC_NULL, &(petscVec.nestVec[0]), &(petscVec.vector));
     } else {
-      int n = 0, ni;
+      PetscInt n = 0;
       for (size_t i = 0; i < nComponents; i++)
 	n += rhs.getDOFVector(i)->getUsedSize();