From 7c35d72dc271127c3d768aa154886ab881b1d86d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 9 Mar 2010 07:22:19 +0000
Subject: [PATCH] Some optimizations.

---
 AMDiS/Reinit/src/ElementLevelSet.h | 179 ++++++++++++-------------
 AMDiS/src/Operator.cc              | 204 +++++++++++++----------------
 AMDiS/src/Operator.h               | 186 +++++++++++++-------------
 3 files changed, 272 insertions(+), 297 deletions(-)

diff --git a/AMDiS/Reinit/src/ElementLevelSet.h b/AMDiS/Reinit/src/ElementLevelSet.h
index a6292e7f..0b291c41 100644
--- a/AMDiS/Reinit/src/ElementLevelSet.h
+++ b/AMDiS/Reinit/src/ElementLevelSet.h
@@ -77,12 +77,11 @@ class ElementLevelSet
     mesh = mesh_;
     dim = mesh->getDim();
 
-    elIntersecPoints = 
-      NEW VectorOfFixVecs<DimVec<double> >(dim, 
-					   MAX_INTERSECTION_POINTS, 
-					   NO_INIT);
-    elVertexStatusVec = new int[dim+1];
-    elVertexLevelSetVec = new double[dim+1];
+    elIntersecPoints = new VectorOfFixVecs<DimVec<double> >(dim, 
+							    MAX_INTERSECTION_POINTS, 
+							    NO_INIT);
+    elVertexStatusVec = new int[dim + 1];
+    elVertexLevelSetVec = new double[dim + 1];
 
     int setElementLevelSetTol = 0;
     GET_PARAMETER(0, name + "->set ElementLevelSet tolerances", "%d", 
@@ -106,7 +105,7 @@ class ElementLevelSet
     if(elVertexLevelSetVec)
       delete [] elVertexLevelSetVec;
     if (elIntersecPoints)
-      DELETE elIntersecPoints;
+      delete elIntersecPoints;
   }
 
   /**
@@ -138,18 +137,15 @@ class ElementLevelSet
   int createElementLevelSet(const ElInfo *elInfo_, 
 			    const bool doCalcIntersecPts_ = true);
 
-  /**
-   * Gets value of level set function at point given in 
-   * barycentric coordinates.
-   */
-  inline double calcLevelSetFct(const DimVec<double>& bary) {
+  /// Gets value of level set function at point given in barycentric coordinates.
+  inline double calcLevelSetFct(const DimVec<double>& bary) 
+  {
     return (*lSFct)(bary);
-  };
+  }
 
-  /**
-   * Resets level set information on element.
-   */
-  inline void resetElement() {
+  /// Resets level set information on element.
+  inline void resetElement() 
+  {
     FUNCNAME("ElementLevelSet::resetElement");
 
     numElVertexInterior = 0;
@@ -157,116 +153,109 @@ class ElementLevelSet
     numElVertexExterior = 0;
     numIntersecPoints = 0;
     elStatus = LEVEL_SET_UNDEFINED;
-  };
+  }
 
-  /**
-   * Defines current element (elInfo).
-   */
-  inline void setElement(const ElInfo *elInfo_) {
+  /// Defines current element (elInfo).
+  inline void setElement(const ElInfo *elInfo_) 
+  {
     elInfo = elInfo_;
     resetElement();
-  };
-
-  /**
-   * Set level_set_domain.
-   */
-  inline void setLevelSetDomain(int status_) {
+  }
 
+  /// Set level_set_domain.
+  inline void setLevelSetDomain(int status_) 
+  {
     TEST_EXIT(status_ == LEVEL_SET_INTERIOR ||
 	      status_ == LEVEL_SET_EXTERIOR ||
 	      status_ == LEVEL_SET_BOUNDARY)("illegal level set status !\n");
     level_set_domain = status_;
-  };
+  }
 
-  /**
-   * Functions to set tolerances for intersection point calculation.
-   */
-  inline void setLsValTol(double tol) {LS_VAL_TOL = tol;};
-  inline void setLsValMin(double min) {LS_VAL_MIN = min;};
-  inline void setSpBaryTol(double tol) {SP_BARY_TOL = tol;};
+  /// Functions to set tolerances for intersection point calculation.
+  inline void setLsValTol(double tol) 
+  {
+    LS_VAL_TOL = tol;
+  }
 
-  /**
-   * Get level_set_domain.
-   */
-  inline const int& getLevelSetDomain() const {
+  inline void setLsValMin(double min) 
+  {
+    LS_VAL_MIN = min;
+  }
+
+  inline void setSpBaryTol(double tol) 
+  {
+    SP_BARY_TOL = tol;
+  }
+
+  /// Get level_set_domain.
+  inline const int& getLevelSetDomain() const 
+  {
     return level_set_domain;
-  };
+  }
 
-  /**
-   * Get LevelSet-Status of element.
-   */
-  inline const int& getElementLevelSetStatus() const {
+  /// Get LevelSet-Status of element.
+  inline const int& getElementLevelSetStatus() const 
+  {
     return elStatus;
-  };
+  }
 
-  /**
-   * Get number of vertices which are intersection points.
-   */
-  inline const int& getNumVertIntPoints() const {
+  /// Get number of vertices which are intersection points.
+  inline const int& getNumVertIntPoints() const 
+  {
     FUNCNAME("ElementLevelSet::getNumVertIntPoints");
     TEST_EXIT(numElVertexBoundary == 0)("numElVertexBoundary should be zero!\n");
     return numElVertexBoundary;
-  }; 
+  };
 
-  /**
-   * Get vector elVertexStatusVec.
-   */
-  inline const int *getElVertStatusVec() const {
+  /// Get vector elVertexStatusVec.
+  inline const int *getElVertStatusVec() const 
+  {
     return elVertexStatusVec;
-  };
+  }
 
-  /**
-   * Get i-th component of vector elVertexLevelSetVec.
-   */
-  inline const double getElVertLevelSetVec(const int i) const {
+  /// Get i-th component of vector elVertexLevelSetVec.
+  inline const double getElVertLevelSetVec(const int i) const 
+  {
     return elVertexLevelSetVec[i];
-  };
+  }
 
-  /**
-   * Get vector elVertexLevelSetVec.
-   */
-  inline const double *getElVertLevelSetVec() const {
+  /// Get vector elVertexLevelSetVec.
+  inline const double *getElVertLevelSetVec() const 
+  {
     return elVertexLevelSetVec;
-  };
+  }
 
-  /**
-   * Get levelSetFct.
-   */
-  inline ElementFunction<double> *getLevelSetFct() const {
+  /// Get levelSetFct.
+  inline ElementFunction<double> *getLevelSetFct() const 
+  {
     return lSFct;
-  };
+  }
 
-  /**
-   * Get mesh.
-   */
-  inline Mesh *getMesh() const {
+  /// Get mesh.
+  inline Mesh *getMesh() const 
+  {
     return mesh;
-  };
+  }
 
-  /**
-   * Get dim.
-   */
-  inline int getDim() const {
+  /// Get dim.
+  inline int getDim() const 
+  {
     return dim;
-  };
+  }
 
-  /**
-   * Get the intersection points.
-   */
-  inline VectorOfFixVecs<DimVec<double> > *getElIntersecPoints() const {
+  /// Get the intersection points.
+  inline VectorOfFixVecs<DimVec<double> > *getElIntersecPoints() const 
+  {
     return elIntersecPoints;
-  };
+  }
 
-  /**
-   * Get number of intersection points.
-   */
-  inline int getNumElIntersecPoints() const {
+  /// Get number of intersection points.
+  inline int getNumElIntersecPoints() const 
+  {
     return numIntersecPoints;
-  };
+  }
 
-  /**
-   * Calculate exterior normal to intersection plane.
-   */
+  /// Calculate exterior normal to intersection plane.
   void calcIntersecNormal(WorldVector<double> &normal);
 
   /**
@@ -287,9 +276,7 @@ class ElementLevelSet
   int getVertexPos(const DimVec<double> barCoords);
 
  protected:
-  /**
-   * Calculates level set value of each vertex of element.
-   */
+  /// Calculates level set value of each vertex of element.
   void calculateElementLevelSetVal();
 
   /**
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 07ee6458..ea95dfbd 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -565,17 +565,15 @@ namespace AMDiS {
 				  const WorldVector<double> *grdUhAtQP,
 				  const WorldMatrix<double> *D2UhAtQP,
 				  double *result,
-				  double fac) const
+				  double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
-	fac * 
-	(*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]) * 
+	fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]) * 
 	uhAtQP[iq];
   }
 
-  void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints,
-				 std::vector<double> &C) const 
+  void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]);
@@ -605,8 +603,7 @@ namespace AMDiS {
     gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
   }
   
-  void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				    std::vector<double> &C) const 
+  void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]);
@@ -617,7 +614,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
@@ -643,6 +640,9 @@ namespace AMDiS {
 
       auxFeSpaces.insert(grads[i]->getFESpace());
     }   
+
+    vecsArg.resize(vecs_.size());
+    gradsArg.resize(grads_.size());
   }
 
   GeneralParametric_ZOT::GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
@@ -758,7 +758,7 @@ namespace AMDiS {
 				 const WorldVector<double> *grdUhAtQP,
 				 const WorldMatrix<double> *D2UhAtQP,
 				 double *result,
-				 double fac) const
+				 double fac) 
   {   
     if (grdUhAtQP)
       for (int iq = 0; iq < nPoints; iq++)
@@ -798,11 +798,11 @@ namespace AMDiS {
   }
   
   void FctVecAtQP_FOT::eval(int nPoints,
-			    const double              *uhAtQP,
+			    const double *uhAtQP,
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double fac) const
+			    double fac) 
   {
     if (grdUhAtQP)
       for (int iq = 0; iq < nPoints; iq++)
@@ -857,7 +857,7 @@ namespace AMDiS {
 			  const WorldVector<double> *grdUhAtQP,
 			  const WorldMatrix<double> *D2UhAtQP,
 			  double *result,
-			  double fac) const
+			  double fac) 
   {
     if (grdUhAtQP)
       for (int iq = 0; iq < nPoints; iq++) 
@@ -926,7 +926,7 @@ namespace AMDiS {
 			     const WorldVector<double> *grdUhAtQP,
 			     const WorldMatrix<double> *D2UhAtQP,
 			     double *result,
-			     double fac) const
+			     double fac) 
   {
     if (grdUhAtQP)
       for (int iq = 0; iq < nPoints; iq++)
@@ -1156,7 +1156,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double factor) const
+			    double factor) 
   {
     int dow = Global::getGeo(WORLD);
     
@@ -1683,7 +1683,7 @@ namespace AMDiS {
       lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);    
   }
 
-  void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+  void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     if (f) {
       for (int iq = 0; iq < nPoints; iq++)
@@ -1699,7 +1699,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double fac) const
+			 double fac) 
   {
     if (f) {
       for (int iq = 0; iq < nPoints; iq++)
@@ -1711,7 +1711,7 @@ namespace AMDiS {
   }
 
 
-  void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+  void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);    
@@ -1722,13 +1722,13 @@ namespace AMDiS {
 			     const WorldVector<double> *grdUhAtQP,
 			     const WorldMatrix<double> *D2UhAtQP,
 			     double *result,
-			     double fac) const
+			     double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
   }
 
-  void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+  void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);    
@@ -1739,13 +1739,13 @@ namespace AMDiS {
 			  const WorldVector<double> *grdUhAtQP,
 			  const WorldMatrix<double> *D2UhAtQP,
 			  double *result,
-			  double fac) const
+			  double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];    
   }
 
-  void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+  void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);    
@@ -1756,7 +1756,7 @@ namespace AMDiS {
 			  const WorldVector<double> *grdUhAtQP,
 			  const WorldMatrix<double> *D2UhAtQP,
 			  double *result,
-			  double fac) const
+			  double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
@@ -1764,8 +1764,7 @@ namespace AMDiS {
   }
 
 
-  void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, 
-				  std::vector<double> &C) const 
+  void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);    
@@ -1776,15 +1775,14 @@ namespace AMDiS {
 				  const WorldVector<double> *grdUhAtQP,
 				  const WorldMatrix<double> *D2UhAtQP,
 				  double *result,
-				  double fac) const
+				  double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq];    
   }
 
  
-  void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, 
-				   std::vector<double> &C) const 
+  void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);    
@@ -1795,15 +1793,14 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
 	fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq];    
   }
 
-  void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				   std::vector<double> &C) const 
+  void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
@@ -1814,7 +1811,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
@@ -1828,7 +1825,7 @@ namespace AMDiS {
 				 const WorldVector<double> *grdUhAtQP,
 				 const WorldMatrix<double> *D2UhAtQP,
 				 double *result,
-				 double fac) const
+				 double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
@@ -1837,15 +1834,13 @@ namespace AMDiS {
 	uhAtQP[iq];
   }
 
-  void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				 std::vector<double> &C) const 
+  void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)  
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
   }
 
-  void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				std::vector<double> &C) const 
+  void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) 
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -1856,14 +1851,13 @@ namespace AMDiS {
 				const WorldVector<double> *grdUhAtQP,
 				const WorldMatrix<double> *D2UhAtQP,
 				double *result,
-				double fac) const
+				double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];   
   }
 
-  void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				   std::vector<double> &C) const 
+  void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) 
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -1874,14 +1868,13 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
   }
 
-  void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints,
-				    std::vector<double> &C) const 
+  void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
@@ -1892,14 +1885,14 @@ namespace AMDiS {
 				    const WorldVector<double> *grdUhAtQP,
 				    const WorldMatrix<double> *D2UhAtQP,
 				    double *result,
-				    double fac) const
+				    double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += 
 	fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * uhAtQP[iq];  
   }
 
-  void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+  void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) 
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*f)(gradAtQPs[iq]);
@@ -1910,14 +1903,13 @@ namespace AMDiS {
 			     const WorldVector<double> *grdUhAtQP,
 			     const WorldMatrix<double> *D2UhAtQP,
 			     double *result,
-			     double fac) const
+			     double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*f)(gradAtQPs[iq]) * uhAtQP[iq];    
   }
 
-  void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, 
-			    std::vector<double> &C) const 
+  void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C)
   { 
     for (int iq = 0; iq < nPoints; iq++)
       C[iq] += (*g)(coordsAtQPs[iq]);    
@@ -1928,7 +1920,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double fac) const
+			    double fac) 
   {
     for (int iq = 0; iq < nPoints; iq++)
       result[iq] += fac * (*g)(coordsAtQPs[iq]) * uhAtQP[iq];    
@@ -1939,7 +1931,7 @@ namespace AMDiS {
 			   const WorldVector<double> *grdUhAtQP,
 			   const WorldMatrix<double> *D2UhAtQP,
 			   double *result,
-			   double factor) const
+			   double factor) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -1977,7 +1969,7 @@ namespace AMDiS {
 			const WorldVector<double> *grdUhAtQP,
 			const WorldMatrix<double> *D2UhAtQP,
 			double *result,
-			double factor) const
+			double factor) 
   {
     int dow = Global::getGeo(WORLD);
   
@@ -2006,7 +1998,7 @@ namespace AMDiS {
 				const WorldVector<double> *grdUhAtQP,
 				const WorldMatrix<double> *D2UhAtQP,
 				double *result,
-				double factor) const
+				double factor) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2050,7 +2042,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double fac) const
+			 double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2080,7 +2072,7 @@ namespace AMDiS {
 			  const WorldVector<double> *grdUhAtQP,
 			  const WorldMatrix<double> *D2UhAtQP,
 			  double *result,
-			  double fac) const
+			  double fac) 
   {
     int dow = Global::getGeo(WORLD);
     
@@ -2110,7 +2102,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double f) const
+			    double f) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2140,7 +2132,7 @@ namespace AMDiS {
 			     const WorldVector<double> *grdUhAtQP,
 			     const WorldMatrix<double> *D2UhAtQP,
 			     double *result,
-			     double fac) const
+			     double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2170,7 +2162,7 @@ namespace AMDiS {
 				    const WorldVector<double> *grdUhAtQP,
 				    const WorldMatrix<double> *D2UhAtQP,
 				    double *result,
-				    double fac) const
+				    double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2200,7 +2192,7 @@ namespace AMDiS {
 				       const WorldVector<double> *grdUhAtQP,
 				       const WorldMatrix<double> *D2UhAtQP,
 				       double *result,
-				       double factor) const
+				       double factor) 
   {
     int dow = Global::getGeo(WORLD);
     
@@ -2237,7 +2229,7 @@ namespace AMDiS {
 				  const WorldVector<double> *grdUhAtQP,
 				  const WorldMatrix<double> *D2UhAtQP,
 				  double *result,
-				  double fac) const
+				  double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2267,7 +2259,7 @@ namespace AMDiS {
 					 const WorldVector<double> *grdUhAtQP,
 					 const WorldMatrix<double> *D2UhAtQP,
 					 double *result,
-					 double factor) const
+					 double factor) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2304,7 +2296,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     int  dow = Global::getGeo(WORLD);
 
@@ -2334,7 +2326,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double fac) const
+			 double fac) 
   {
     int dow = Global::getGeo(WORLD);
   
@@ -2354,7 +2346,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double f) const
+			    double f) 
   {
     int dow = Global::getGeo(WORLD);
   
@@ -2374,7 +2366,7 @@ namespace AMDiS {
 			       const WorldVector<double> *grdUhAtQP,
 			       const WorldMatrix<double> *D2UhAtQP,
 			       double *result,
-			       double f) const
+			       double f) 
   {
     int dow = Global::getGeo(WORLD);
   
@@ -2395,7 +2387,7 @@ namespace AMDiS {
 				const WorldVector<double> *grdUhAtQP,
 				const WorldMatrix<double> *D2UhAtQP,
 				double *result,
-				double factor) const
+				double factor) 
   { 
     if (grdUhAtQP) {
       if (f) {
@@ -2415,7 +2407,7 @@ namespace AMDiS {
 			   const WorldVector<double> *grdUhAtQP,
 			   const WorldMatrix<double> *D2UhAtQP,
 			   double *result,
-			   double factor) const
+			   double factor) 
   {
     if (grdUhAtQP) {
       for (int iq = 0; iq < nPoints; iq++) {
@@ -2430,7 +2422,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double f) const
+			    double f) 
   {
     int dow = Global::getGeo(WORLD);
   
@@ -2479,7 +2471,7 @@ namespace AMDiS {
 			       const WorldVector<double> *grdUhAtQP,
 			       const WorldMatrix<double> *D2UhAtQP,
 			       double *result,
-			       double f) const
+			       double f) 
   {
     if (D2UhAtQP) {
       for (int iq = 0; iq < nPoints; iq++) {
@@ -2521,7 +2513,7 @@ namespace AMDiS {
 			    const WorldVector<double> *grdUhAtQP,
 			    const WorldMatrix<double> *D2UhAtQP,
 			    double *result,
-			    double fac) const
+			    double fac) 
   {
     if (D2UhAtQP) {
       for (int iq = 0; iq < nPoints; iq++) {
@@ -2551,8 +2543,7 @@ namespace AMDiS {
       vecsAtQPs[i] = getVectorAtQPs(vecs[i], elInfo, subAssembler, quad);
   }
  
-  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				  std::vector<double> &C) const 
+  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,  std::vector<double> &C)
   { 
     int size = static_cast<int>(vecs.size());
     std::vector<double> arg(size);
@@ -2570,7 +2561,7 @@ namespace AMDiS {
 				  const WorldVector<double> *grdUhAtQP,
 				  const WorldMatrix<double> *D2UhAtQP,
 				  double *result,
-				  double fac) const
+				  double fac) 
   {
     int size = static_cast<int>(vecs.size());
     std::vector<double> arg(size);
@@ -2603,8 +2594,7 @@ namespace AMDiS {
   }
  
 
-  void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
-				    std::vector<double> &C) const 
+  void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
   { 
     int size = static_cast<int>(vecs.size());
     std::vector<WorldVector<double>*> arg(size);
@@ -2622,7 +2612,7 @@ namespace AMDiS {
 				    const WorldVector<double> *grdUhAtQP,
 				    const WorldMatrix<double> *D2UhAtQP,
 				    double *result,
-				    double fac) const
+				    double fac) 
   {
     int size = static_cast<int>(vecs.size());
     std::vector<WorldVector<double>*> arg(size);
@@ -2635,8 +2625,7 @@ namespace AMDiS {
     }
   }
 
-  void VecDivergence_ZOT::getC(const ElInfo *, int nPoints,
-			       std::vector<double> &C) const 
+  void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
   { 
     int size = static_cast<int>(vecs.size());
 
@@ -2650,7 +2639,7 @@ namespace AMDiS {
 			       const WorldVector<double> *grdUhAtQP,
 			       const WorldMatrix<double> *D2UhAtQP,
 			       double *result,
-			       double fac) const
+			       double fac) 
   {
     int size = static_cast<int>(vecs.size());
 
@@ -2668,7 +2657,7 @@ namespace AMDiS {
 				const WorldVector<double> *grdUhAtQP,
 				const WorldMatrix<double> *D2UhAtQP,
 				double *result,
-				double fac) const
+				double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -2753,7 +2742,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double factor) const
+			 double factor) 
   {
     int dow = Global::getGeo(WORLD);
     int nVecs = static_cast<int>(vecs_.size());
@@ -2850,7 +2839,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double factor) const
+			 double factor) 
   {
     int dow = Global::getGeo(WORLD);  
     int nVecs = static_cast<int>(vecs_.size());
@@ -2892,18 +2881,15 @@ namespace AMDiS {
       gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
   }
 
-  void General_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
+  void General_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
   {
-    int nVecs = static_cast<int>(vecs_.size());
-    int nGrads = static_cast<int>(grads_.size());
-
-    std::vector<double> vecsArg(nVecs);
-    std::vector<WorldVector<double> > gradsArg(nGrads);
+    unsigned int nVecs = vecs_.size();
+    unsigned int nGrads = grads_.size();
 
     for (int iq = 0; iq < nPoints; iq++) {
-      for (int i = 0; i < nVecs; i++)
+      for (unsigned int i = 0; i < nVecs; i++)
 	vecsArg[i] = vecsAtQPs_[i][iq];
-      for (int i = 0; i < nGrads; i++)
+      for (unsigned int i = 0; i < nGrads; i++)
 	gradsArg[i] = gradsAtQPs_[i][iq];
 
       C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
@@ -2915,19 +2901,17 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double fac) const
+			 double fac) 
   {
-    int nVecs = static_cast<int>(vecs_.size());
-    int nGrads = static_cast<int>(grads_.size());
-
-    std::vector<double> vecsArg(nVecs);
-    std::vector<WorldVector<double> > gradsArg(nGrads);
+    unsigned int nVecs = vecs_.size();
+    unsigned int nGrads = grads_.size();
 
     for (int iq = 0; iq < nPoints; iq++) {
-      for (int i = 0; i < nVecs; i++)
+      for (unsigned int i = 0; i < nVecs; i++)
 	vecsArg[i] = vecsAtQPs_[i][iq];      
-      for (int i = 0; i < nGrads; i++) 
+      for (unsigned int i = 0; i < nGrads; i++) 
 	gradsArg[i] = gradsAtQPs_[i][iq];      
+
       result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq];
     }
   }
@@ -2974,7 +2958,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double factor) const
+				   double factor) 
   {
     int dow = Global::getGeo(WORLD);
     int nVecs = static_cast<int>(vecs_.size());
@@ -3071,7 +3055,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double factor) const
+				   double factor) 
   {
     int dow = Global::getGeo(WORLD);
     int nVecs = static_cast<int>(vecs_.size());
@@ -3114,8 +3098,7 @@ namespace AMDiS {
       gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
   }
 
-  void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, 
-				   std::vector<double> &C) const
+  void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
   {
     int nVecs = static_cast<int>(vecs_.size());
     int nGrads = static_cast<int>(grads_.size());
@@ -3124,12 +3107,11 @@ namespace AMDiS {
     std::vector<WorldVector<double> > gradsArg(nGrads);
 
     for (int iq = 0; iq < nPoints; iq++) {
-      for (int i = 0; i < nVecs; i++) {
-	vecsArg[i] = vecsAtQPs_[i][iq];
-      }
-      for (int i = 0; i < nGrads; i++) {
+      for (int i = 0; i < nVecs; i++)
+	vecsArg[i] = vecsAtQPs_[i][iq];      
+      for (int i = 0; i < nGrads; i++)
 	gradsArg[i] = gradsAtQPs_[i][iq];
-      }
+      
       C[iq] += (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
     }
   }
@@ -3139,7 +3121,7 @@ namespace AMDiS {
 				   const WorldVector<double> *grdUhAtQP,
 				   const WorldMatrix<double> *D2UhAtQP,
 				   double *result,
-				   double fac) const
+				   double fac) 
   {
     int nVecs = static_cast<int>(vecs_.size());
     int nGrads = static_cast<int>(grads_.size());
@@ -3169,7 +3151,7 @@ namespace AMDiS {
   }
 
   void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
-					  std::vector<double> &C) const
+					  std::vector<double> &C)
   {
     int size = static_cast<int>(vecs.size());
     std::vector<WorldVector<double>*> arg(size);
@@ -3187,7 +3169,7 @@ namespace AMDiS {
 					  const WorldVector<double> *grdUhAtQP,
 					  const WorldMatrix<double> *D2UhAtQP,
 					  double *result,
-					  double fac) const
+					  double fac) 
   {
     int size = static_cast<int>(vecs.size());
     std::vector<WorldVector<double>*> arg(size);
@@ -3229,7 +3211,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double factor) const
+			 double factor) 
   {
     if (grdUhAtQP) {
       for (int iq = 0; iq < nPoints; iq++) {
@@ -3244,7 +3226,7 @@ namespace AMDiS {
 			 const WorldVector<double> *grdUhAtQP,
 			 const WorldMatrix<double> *D2UhAtQP,
 			 double *result,
-			 double fac) const
+			 double fac) 
   {
     int dow = Global::getGeo(WORLD);
 
@@ -3306,7 +3288,7 @@ namespace AMDiS {
 			  const WorldVector<double> *grdUhAtQP,
 			  const WorldMatrix<double> *D2UhAtQP,
 			  double *result,
-			  double factor) const
+			  double factor) 
   {
     if (grdUhAtQP) {
       for (int iq = 0; iq < nPoints; iq++) {
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 93482fbf..29896928 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -110,7 +110,7 @@ namespace AMDiS {
 		      const WorldVector<double> *grdUhAtQP,
 		      const WorldMatrix<double> *D2UhAtQP,
 		      double *result,
-		      double factor) const = 0;
+		      double factor) = 0;
 
     /** \brief
      * Determines the value of a dof vector at the quadrature points of a given 
@@ -335,7 +335,7 @@ namespace AMDiS {
 		     const WorldVector<double> *,
 		     const WorldMatrix<double> *D2UhAtQP,
 		     double *result,
-		     double factor) const
+		     double factor) 
     {
       int dow = Global::getGeo(WORLD);
     
@@ -401,7 +401,7 @@ namespace AMDiS {
 	      const WorldVector<double> *,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double f) const
+	      double f) 
     {
       int dow = Global::getGeo(WORLD);
 
@@ -460,7 +460,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -514,7 +514,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -566,11 +566,11 @@ namespace AMDiS {
 
     /// Implenetation of SecondOrderTerm::eval().
     void eval(int nPoints,
-	      const double              *,
+	      const double *,
 	      const WorldVector<double> *,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const
+	      double fac)
     {
       if (D2UhAtQP)
 	for (int iq = 0; iq < nPoints; iq++)
@@ -626,7 +626,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -670,7 +670,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -718,7 +718,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -769,7 +769,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -821,7 +821,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -869,7 +869,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -919,7 +919,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -973,7 +973,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -1023,7 +1023,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
 		  std::vector<WorldVector<double> > &result) const;
@@ -1074,7 +1074,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -1123,7 +1123,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -1175,7 +1175,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -1232,7 +1232,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implenetation of SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -1295,7 +1295,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *,
 	      double *result,
-	      double factor) const 
+	      double factor)  
     {
       int dow = Global::getGeo(WORLD);
 
@@ -1419,7 +1419,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *,
 	      double *result,
-	      double factor) const
+	      double factor)
     {
       if (grdUhAtQP)
 	for (int iq = 0; iq < nPoints; iq++)
@@ -1464,7 +1464,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -1507,7 +1507,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     /// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -1546,7 +1546,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     /// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -1586,7 +1586,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     DOFVectorBase<double>* vec;
@@ -1625,7 +1625,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -1666,7 +1666,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     /// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -1710,7 +1710,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -1756,7 +1756,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -1779,22 +1779,24 @@ namespace AMDiS {
   class Vec2AndGradAtQP_FOT : public FirstOrderTerm
   {
   public:
-
     Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1, 
 			DOFVectorBase<double> *dv2,
 			TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_);
 
-    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+    void initElement(const ElInfo* elInfo, 
+		     SubAssembler* subAssembler,
 		     Quadrature *quad = NULL);
     
-    void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
+    void getLb(const ElInfo *elInfo, 
+	       int nPoints, 
+	       VectorOfFixVecs<DimVec<double> >& Lb) const;
     
     void eval(int nPoints,
-	      const double              *uhAtQP,
+	      const double *uhAtQP,
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     DOFVectorBase<double>* vec1;
@@ -1827,7 +1829,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
   protected:
     DOFVectorBase<double>* vec;
@@ -1858,7 +1860,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
   protected:
     DOFVectorBase<double>* vec1;
@@ -1891,7 +1893,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
     
   protected:
     DOFVectorBase<double>* vec1;
@@ -1931,7 +1933,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     std::vector<DOFVectorBase<double>*> vecs_; 
@@ -1975,7 +1977,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
   protected:
     std::vector<DOFVectorBase<double>*> vecs_; 
@@ -2018,7 +2020,7 @@ namespace AMDiS {
 
     /// Evaluates \f$ c \f$
     virtual void getC(const ElInfo *elInfo, int nPoints, 
-		      std::vector<double> &C) const = 0;
+		      std::vector<double> &C) = 0;
 
   };
 
@@ -2035,7 +2037,7 @@ namespace AMDiS {
     Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {}
 
     /// Implements ZeroOrderTerm::getC().
-    inline void getC(const ElInfo *, int nPoints, std::vector<double> &C) const 
+    inline void getC(const ElInfo *, int nPoints, std::vector<double> &C)  
     {
       for (int iq = 0; iq < nPoints; iq++)
 	C[iq] += factor; 
@@ -2047,7 +2049,7 @@ namespace AMDiS {
 		     const WorldVector<double> *,
 		     const WorldMatrix<double> *,
 		     double *result,
-		     double fac) const 
+		     double fac)  
     {
       for (int iq = 0; iq < nPoints; iq++)
 	result[iq] += fac * factor * uhAtQP[iq];
@@ -2083,7 +2085,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2091,7 +2093,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2124,7 +2126,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2132,7 +2134,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVectorBase to be evaluated at quadrature points.
@@ -2173,7 +2175,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2181,7 +2183,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// First DOFVector to be evaluated at quadrature points.
@@ -2219,7 +2221,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2227,7 +2229,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVectors to be evaluated at quadrature points.
@@ -2259,7 +2261,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements SecondOrderTerm::getC().
-    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2267,7 +2269,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     DOFVectorBase<double>* vec;
@@ -2302,7 +2304,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2310,7 +2312,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2345,7 +2347,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2353,7 +2355,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2387,7 +2389,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2395,7 +2397,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2432,7 +2434,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements SecondOrderTerm::getC().
-    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2440,7 +2442,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     DOFVectorBase<double>* vec;
@@ -2473,7 +2475,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2481,7 +2483,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2523,7 +2525,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implements SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -2573,7 +2575,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implements SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -2615,7 +2617,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2623,7 +2625,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -2665,7 +2667,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implements SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -2711,7 +2713,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double factor) const;
+	      double factor);
 
     /// Implements SecondOrderTerm::weakEval().
     void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -2745,7 +2747,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2753,7 +2755,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     /// DOFVector to be evaluated at quadrature points. 
@@ -2786,7 +2788,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2794,7 +2796,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     /// DOFVector to be evaluated at quadrature points. 
@@ -2832,7 +2834,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2840,7 +2842,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     /// Vector of DOFVectors to be evaluated at quadrature points. 
@@ -2865,7 +2867,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2873,7 +2875,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     /// Vector of DOFVectors to be evaluated at quadrature points. 
@@ -2901,7 +2903,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2909,7 +2911,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     /// Vector of DOFVectors to be evaluated at quadrature points. 
@@ -2934,7 +2936,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -2942,7 +2944,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     /// DOFVector to be evaluated at quadrature points.
@@ -2972,14 +2974,14 @@ namespace AMDiS {
     void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
 		     Quadrature *quad = NULL);
 
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     void eval(int nPoints,
 	      const double *uhAtQP,
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     DOFVectorBase<double> *vec1, *vec2;
@@ -3002,14 +3004,14 @@ namespace AMDiS {
     void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
 		     Quadrature *quad = NULL);
 
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     void eval(int nPoints,
 	      const double *uhAtQP,
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected: 
     DOFVectorBase<double>* vec1; 
@@ -3039,7 +3041,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -3047,7 +3049,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     std::vector<DOFVectorBase<double>*> vecs_; 
@@ -3064,6 +3066,10 @@ namespace AMDiS {
     std::vector<double*> vecsAtQPs_;
 
     std::vector<WorldVector<double>*> gradsAtQPs_;
+
+    std::vector<double> vecsArg;
+    
+    std::vector<WorldVector<double> > gradsArg;
   };
 
   class GeneralParametric_ZOT : public ZeroOrderTerm
@@ -3083,7 +3089,7 @@ namespace AMDiS {
 		     Quadrature *quad = NULL);
 
     /// Implements ZeroOrderTerm::getC().
-    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+    void getC(const ElInfo *, int nPoints, std::vector<double> &C);
 
     /// Implements ZeroOrderTerm::eval().
     void eval(int nPoints,
@@ -3091,7 +3097,7 @@ namespace AMDiS {
 	      const WorldVector<double> *grdUhAtQP,
 	      const WorldMatrix<double> *D2UhAtQP,
 	      double *result,
-	      double fac) const;
+	      double fac);
 
   protected:
     std::vector<DOFVectorBase<double>*> vecs_; 
@@ -3407,7 +3413,7 @@ namespace AMDiS {
     }
 
     /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
-    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
+    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c)
     {
       int myRank = omp_get_thread_num();
 
-- 
GitLab