* Several optimizations

AMDiS/src/AdaptInstationary.cc

@@ -205,13 +205,17 @@ namespace AMDiS {
 
       problemTime_->initTimestep(adaptInfo_);
 
-      #ifdef _OPENMP
-#pragma omp parallel sections
-      {
+#ifdef _OPENMP
+      if (problemTime_->existsDelayedCalculation()) {
+#pragma omp parallel sections num_threads(2)
+	{
 #pragma omp section
-	problemTime_->startDelayedTimestepCalculation();
+	  problemTime_->startDelayedTimestepCalculation();
 
 #pragma omp section
+	  oneTimestep();
+	}
+      } else {
 	oneTimestep();
       }
 #else     

AMDiS/src/Assembler.cc

@@ -30,6 +30,7 @@ namespace AMDiS {
     : nRow(0),
       nCol(0),
       coordsAtQPs(NULL),
+      coordsNumAllocated(0),
       quadrature(quadrat),
       psiFast(NULL),
       phiFast(NULL),
@@ -128,16 +129,20 @@ namespace AMDiS {
       return coordsAtQPs;
     }
    
-    // not yet calcualted !
-    if (coordsAtQPs) 
-      DELETE [] coordsAtQPs;
-    // allocate memory
-    coordsAtQPs = NEW WorldVector<double>[numPoints];
+    if (coordsAtQPs)  {
+      if (coordsNumAllocated != numPoints) {
+	DELETE [] coordsAtQPs;
+        coordsAtQPs = NEW WorldVector<double>[numPoints];
+	coordsNumAllocated = numPoints;
+      }
+    } else {
+      coordsAtQPs = NEW WorldVector<double>[numPoints];
+      coordsNumAllocated = numPoints;
+    }
 
     // set new values
-    WorldVector<double>* k = NULL;
-    int l;
-    for (k = &(coordsAtQPs[0]), l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
+    WorldVector<double>* k = &(coordsAtQPs[0]);
+    for (int l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
       elInfo->coordToWorld(localQuad->getLambda(l), k);
     }
 

AMDiS/src/Assembler.h

@@ -245,6 +245,11 @@ namespace AMDiS {
      */
     bool coordsValid;
 
+    /** \brief
+     * Used for \ref getCoordsAtQP(). Stores the number of allocated WorldVectors.
+     */
+    int coordsNumAllocated;
+
     /** \brief
      * Needed Quadrature. Constructed in the constructor of SubAssembler
      */

AMDiS/src/ElInfo.h

@@ -414,7 +414,7 @@ namespace AMDiS {
      * value of the determinant of the affine linear paraetrization's Jacobian.
      * pure virtual => must be overriden in sub-class.
      */
-    virtual double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const = 0;
+    virtual double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) = 0;
 
     /** \brief
      * calculates a normal of the given side (1d,2d: edge, 3d: face) of \ref element.
@@ -422,7 +422,7 @@ namespace AMDiS {
      * transformation to the reference element.
      * pure virtual => must be overriden in sub-class.
      */
-    virtual double getNormal(int side, WorldVector<double> &normal) const = 0;
+    virtual double getNormal(int side, WorldVector<double> &normal) = 0;
 
 
     /** \brief

AMDiS/src/ElInfo1d.cc

@@ -83,7 +83,7 @@ namespace AMDiS {
   /*  value of the determinante from the transformation to the reference      */
   /*  element                                                                 */
   /****************************************************************************/
-  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
+  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
   {
     FUNCNAME("ElInfo1d::calcGrdLambda()");
 
@@ -149,7 +149,7 @@ namespace AMDiS {
   /*  coord; return the absulute value of the determinant from the           */
   /*  transformation to the reference element                                 */
   /****************************************************************************/
-  double ElInfo1d::getNormal(int side, WorldVector<double> &normal) const
+  double ElInfo1d::getNormal(int side, WorldVector<double> &normal)
   {
     normal = coord_[side] - coord_[(side + 1) % 2];
     double det = norm(&normal);

AMDiS/src/ElInfo1d.h

@@ -63,12 +63,12 @@ namespace AMDiS {
     /** \brief
      * 1-dimensional realisation of ElInfo's calcGrdLambda method.
      */
-    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const;
+    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam);
 
     /** \brief
      * 1-dimensional realisation of ElInfo's getNormal method.
      */
-    double getNormal(int side, WorldVector<double> &normal) const;
+    double getNormal(int side, WorldVector<double> &normal);
 
     /** \brief
      * 1-dimensional realisation of ElInfo's getElementNormal method.

AMDiS/src/ElInfo2d.cc

@@ -391,7 +391,7 @@ namespace AMDiS {
     }
   }
 
-  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
+  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
   {
     FUNCNAME("ElInfo2d::calcGrdLambda()");
   
@@ -512,7 +512,7 @@ namespace AMDiS {
   }
 
 
-  double ElInfo2d::getNormal(int side, WorldVector<double> &normal) const
+  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
   {
     FUNCNAME("ElInfo::getNormal()");
 

AMDiS/src/ElInfo2d.h

@@ -63,12 +63,12 @@ namespace AMDiS {
     /** \brief
      * 2-dimensional realisation of ElInfo's calcGrdLambda method.
      */
-    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const;
+    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam);
 
     /** \brief
      * 2-dimensional realisation of ElInfo's getNormal method.
      */
-    double getNormal(int side, WorldVector<double> &normal) const;
+    double getNormal(int side, WorldVector<double> &normal);
 
     /** \brief
      * 2-dimensional realisation of ElInfo's getElementNormal method.

AMDiS/src/ElInfo3d.cc

@@ -107,30 +107,35 @@ namespace AMDiS {
     }
   }
 
-  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
+  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
   {
     FUNCNAME("ElInfo3d::calcGrdLambda()");
 
     TEST_EXIT_DBG(dimOfWorld == 3)
       ("dim != dim_of_world ! use parametric elements!\n");
 
-    WorldVector<double> e1, e2, e3;
-    WorldVector<double> v0;
+    int myRank = omp_get_thread_num();
+
+    WorldVector<double> *e1 = &tmpWorldVecs[myRank][0];
+    WorldVector<double> *e2 = &tmpWorldVecs[myRank][1];
+    WorldVector<double> *e3 = &tmpWorldVecs[myRank][2];
+    WorldVector<double> *v0 = &tmpWorldVecs[myRank][3];
+
     double det, adet;
     double a11, a12, a13, a21, a22, a23, a31, a32, a33;
 
     testFlag(Mesh::FILL_COORDS);
 
-    v0 = coord_[0];
+    *v0 = coord_[0];
     for (int i = 0; i < 3; i++) {
-      e1[i] = coord_[1][i] - v0[i];
-      e2[i] = coord_[2][i] - v0[i];
-      e3[i] = coord_[3][i] - v0[i];
+      (*e1)[i] = coord_[1][i] - (*v0)[i];
+      (*e2)[i] = coord_[2][i] - (*v0)[i];
+      (*e3)[i] = coord_[3][i] - (*v0)[i];
     }
 
-    det =   e1[0] * (e2[1] * e3[2] - e2[2] * e3[1])
-          - e1[1] * (e2[0] * e3[2] - e2[2] * e3[0])
-          + e1[2] * (e2[0] * e3[1] - e2[1] * e3[0]);
+    det = (*e1)[0] * ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1])
+      - (*e1)[1] * ((*e2)[0] * (*e3)[2] - (*e2)[2] * (*e3)[0])
+      + (*e1)[2] * ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]);
 
     adet = abs(det);
 
@@ -141,15 +146,15 @@ namespace AMDiS {
 	  grd_lam[i][j] = 0.0;
     } else {
       det = 1.0 / det;
-      a11 = (e2[1] * e3[2] - e2[2] * e3[1]) * det;    /* (a_ij) = A^{-T} */
-      a12 = (e2[2] * e3[0] - e2[0] * e3[2]) * det;
-      a13 = (e2[0] * e3[1] - e2[1] * e3[0]) * det;
-      a21 = (e1[2] * e3[1] - e1[1] * e3[2]) * det;
-      a22 = (e1[0] * e3[2] - e1[2] * e3[0]) * det;
-      a23 = (e1[1] * e3[0] - e1[0] * e3[1]) * det;
-      a31 = (e1[1] * e2[2] - e1[2] * e2[1]) * det;
-      a32 = (e1[2] * e2[0] - e1[0] * e2[2]) * det;
-      a33 = (e1[0] * e2[1] - e1[1] * e2[0]) * det;
+      a11 = ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1]) * det;    /* (a_ij) = A^{-T} */
+      a12 = ((*e2)[2] * (*e3)[0] - (*e2)[0] * (*e3)[2]) * det;
+      a13 = ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]) * det;
+      a21 = ((*e1)[2] * (*e3)[1] - (*e1)[1] * (*e3)[2]) * det;
+      a22 = ((*e1)[0] * (*e3)[2] - (*e1)[2] * (*e3)[0]) * det;
+      a23 = ((*e1)[1] * (*e3)[0] - (*e1)[0] * (*e3)[1]) * det;
+      a31 = ((*e1)[1] * (*e2)[2] - (*e1)[2] * (*e2)[1]) * det;
+      a32 = ((*e1)[2] * (*e2)[0] - (*e1)[0] * (*e2)[2]) * det;
+      a33 = ((*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0]) * det;
 
       grd_lam[1][0] = a11;
       grd_lam[1][1] = a12;
@@ -161,9 +166,9 @@ namespace AMDiS {
       grd_lam[3][1] = a32;
       grd_lam[3][2] = a33;
 
-      grd_lam[0][0] = -grd_lam[1][0] -grd_lam[2][0] -grd_lam[3][0];
-      grd_lam[0][1] = -grd_lam[1][1] -grd_lam[2][1] -grd_lam[3][1];
-      grd_lam[0][2] = -grd_lam[1][2] -grd_lam[2][2] -grd_lam[3][2];
+      grd_lam[0][0] = -grd_lam[1][0] - grd_lam[2][0] - grd_lam[3][0];
+      grd_lam[0][1] = -grd_lam[1][1] - grd_lam[2][1] - grd_lam[3][1];
+      grd_lam[0][2] = -grd_lam[1][2] - grd_lam[2][2] - grd_lam[3][2];
     }
 
     return adet;
@@ -305,11 +310,17 @@ namespace AMDiS {
   }
 
 
-  double ElInfo3d::getNormal(int face, WorldVector<double> &normal) const
+  double ElInfo3d::getNormal(int face, WorldVector<double> &normal)
   {
-    FUNCNAME("ElInfo3d::getNormal");
+    FUNCNAME("ElInfo3d::getNormal()");
+
     double det = 0.0;
-    WorldVector<double> e0, e1, e2;
+
+    int myRank = omp_get_thread_num();
+
+    WorldVector<double> *e0 = &tmpWorldVecs[myRank][0];
+    WorldVector<double> *e1 = &tmpWorldVecs[myRank][1];
+    WorldVector<double> *e2 = &tmpWorldVecs[myRank][2];
 
     if (dimOfWorld == 3) {
       int i0 = (face + 1) % 4;
@@ -317,14 +328,14 @@ namespace AMDiS {
       int i2 = (face + 3) % 4;
 
       for (int i = 0; i < dimOfWorld; i++) {
-	e0[i] = coord_[i1][i] - coord_[i0][i];
-	e1[i] = coord_[i2][i] - coord_[i0][i];
-	e2[i] = coord_[face][i] - coord_[i0][i];
+	(*e0)[i] = coord_[i1][i] - coord_[i0][i];
+	(*e1)[i] = coord_[i2][i] - coord_[i0][i];
+	(*e2)[i] = coord_[face][i] - coord_[i0][i];
       }
 
-      vectorProduct(e0, e1, normal);
+      vectorProduct(*e0, *e1, normal);
 
-      if ((e2*normal) < 0.0)
+      if ((*e2 * normal) < 0.0)
 	for (int i = 0; i < dimOfWorld; i++)
 	  normal[i] = -normal[i];
 

AMDiS/src/ElInfo3d.h

@@ -24,6 +24,7 @@
 
 #include "ElInfo.h"
 #include "MemoryManager.h"
+#include "OpenMP.h"
 
 namespace AMDiS {
 
@@ -44,7 +45,15 @@ namespace AMDiS {
     /** \brief
      * Constructor. Calls ElInfo's protected Constructor.
      */
-    ElInfo3d(Mesh* aMesh) : ElInfo(aMesh) {};
+    ElInfo3d(Mesh* aMesh) 
+      : ElInfo(aMesh) 
+    {
+      tmpWorldVecs.resize(omp_get_max_threads());
+
+      for (int i = 0; i < static_cast<int>(tmpWorldVecs.size()); i++) {
+	tmpWorldVecs[i].resize(4);
+      }
+    };
 
     /** \brief
      * Assignment operator
@@ -74,12 +83,12 @@ namespace AMDiS {
     /** \brief
      * 3-dimensional realisation of ElInfo's calcGrdLambda method.
      */
-    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const;
+    double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam);
 
     /** \brief
      * 3-dimensional realisation of ElInfo's getNormal method.
      */
-    double getNormal(int side, WorldVector<double> &normal) const;
+    double getNormal(int side, WorldVector<double> &normal);
 
     /** \brief
      * update ElInfo after refinement (of some neighbours). Only in 3d!
@@ -126,6 +135,12 @@ namespace AMDiS {
      * element with vertices (0,0,0), (1,1,1), (1,1,0), (1,0,0).
      */
     signed char orientation;
+
+    /** \brief
+     * Tmp vectors used for calculations in calcGrdLambda and getNormal().
+     * Thread safe!
+     */
+    ::std::vector< ::std::vector< WorldVector<double> > > tmpWorldVecs;
   };
 
 }

AMDiS/src/ElementFileWriter.h

@@ -32,6 +32,10 @@ namespace AMDiS {
 
     void writeDelayedFiles() {};
 
+    bool isWritingDelayed() { 
+      return false; 
+    };
+
   protected:
     /**
      * Writes element data in tecplot format.

AMDiS/src/FileWriter.h

@@ -72,6 +72,8 @@ namespace AMDiS {
 
     virtual void writeDelayedFiles() = 0;
 
+    virtual bool isWritingDelayed() = 0;
+
     void setTraverseProperties(int level, 
 			       Flag flag,
 			       bool (*writeElem)(ElInfo*))
@@ -153,8 +155,18 @@ namespace AMDiS {
 			    Flag traverseFlag = Mesh::CALL_LEAF_EL,
 			    bool (*writeElem)(ElInfo*) = NULL);
 
+    /** \brief
+     * Starts the delayed writing.
+     */
     virtual void writeDelayedFiles();
 
+    /** \brief
+     * Returns true, if the file writer is waiting to start writing.
+     */
+    bool isWritingDelayed() {
+      return writingIsDelayed_;
+    }
+
   protected:
     /** \brief
      * Initialization of the filewriter.

AMDiS/src/Lagrange.cc

@@ -947,13 +947,14 @@ namespace AMDiS {
     
     for (int pos = 0, j = 0; pos <= dim; pos++) {
       posIndex = INDEX_OF_DIM(pos, dim);
-      n0 = admin->getNumberOfPreDOFs(posIndex);
-      node0 = admin->getMesh()->getNode(posIndex);
-      num = Global::getGeo(posIndex, dim);
       nrDOFs = admin->getNumberOfDOFs(posIndex);
 
-      for (int i = 0; i < num; node0++, i++) {
-	if (nrDOFs) {
+      if (nrDOFs) {
+	n0 = admin->getNumberOfPreDOFs(posIndex);
+	node0 = admin->getMesh()->getNode(posIndex);
+	num = Global::getGeo(posIndex, dim);
+
+	for (int i = 0; i < num; node0++, i++) {
 	  indi = orderOfPositionIndices(el, posIndex, i);
 
 	  for (int k = 0; k < nrDOFs; k++)

AMDiS/src/OEMSolver.hh

@@ -96,7 +96,7 @@ namespace AMDiS {
     }
 
     if (res <= tolerance) {
-	INFO(info,6)("finished successfully with %d iterations\n");
+      INFO(info,6)("finished successfully with %d iterations\n", iter);
 	return(1);
     }
 

AMDiS/src/ProblemInstat.cc

@@ -194,4 +194,13 @@ namespace AMDiS {
     problemStat->writeDelayedFiles();
   }
 
+  bool ProblemInstatScal::existsDelayedCalculation()
+  {
+    return problemStat->existsDelayedCalculation();
+  }
+
+  bool ProblemInstatVec::existsDelayedCalculation()
+  {
+    return problemStat->existsDelayedCalculation();
+  }
 }

AMDiS/src/ProblemInstat.h

@@ -206,6 +206,8 @@ namespace AMDiS {
 
     virtual void startDelayedTimestepCalculation();
 
+    virtual bool existsDelayedCalculation();
+
     virtual void serialize(::std::ostream &out) {};
 
 
@@ -295,6 +297,8 @@ namespace AMDiS {
 
     virtual void startDelayedTimestepCalculation();
 
+    virtual bool existsDelayedCalculation();
+
     virtual void serialize(::std::ostream &out) {};
 
     virtual void deserialize(::std::istream &in) {};

AMDiS/src/ProblemScal.cc

@@ -42,6 +42,16 @@ namespace AMDiS {
     }
   }
 
+  bool ProblemScal::existsDelayedCalculation()
+  {
+    for (int i = 0; i < static_cast<int>(fileWriters_.size()); i++) {
+      if (fileWriters_[i]->isWritingDelayed())
+	return true;   
+    }
+
+    return false;
+  }
+
   void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
   {
     solution_->interpol(fct);

AMDiS/src/ProblemScal.h

@@ -200,6 +200,11 @@ namespace AMDiS {
      */
     void writeDelayedFiles();
 
+    /** \brief
+     * Returns true, if there is calculation waiting to be started.
+     */
+    bool existsDelayedCalculation();
+
     /** \brief
      * Interpolates fct to \ref solution.
      */

AMDiS/src/ProblemTimeInterface.h

@@ -75,6 +75,11 @@ namespace AMDiS {
      */
     virtual void startDelayedTimestepCalculation() = 0;
 
+    /** \brief
+     * Returns true, if there is some delayed calculation waiting to be started.
+     */
+    virtual bool existsDelayedCalculation() = 0;
+
     /** \brief
      * Function that serializes the problem plus information about the iteration.
      */

AMDiS/src/ProblemVec.cc

@@ -845,6 +845,16 @@ namespace AMDiS {
     }
   }
 
+  bool ProblemVec::existsDelayedCalculation()
+  {
+    for (int i = 0; i < static_cast<int>(fileWriters_.size()); i++) {
+      if (fileWriters_[i]->isWritingDelayed())
+	return true;   
+    }
+
+    return false;
+  }
+
   void ProblemVec::interpolInitialSolution(::std::vector<AbstractFunction<double, WorldVector<double> >*> *fct) 
   {
     FUNCNAME("ProblemVec::interpolInitialSolution()");