From 44e4efd11027c2d9da76657f7af05a3a39da1f58 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 23 Mar 2009 12:21:46 +0000
Subject: [PATCH] * hope, it works, no more bugs at all

---
 AMDiS/src/CoarseningManager.cc    |  9 ++--
 AMDiS/src/ComponentTraverseInfo.h | 22 ++++++++-
 AMDiS/src/DOFMatrix.cc            | 10 ++--
 AMDiS/src/DOFVector.cc            |  6 +++
 AMDiS/src/DOFVector.h             |  2 +-
 AMDiS/src/DOFVector.hh            |  1 +
 AMDiS/src/Operator.h              | 27 +++--------
 AMDiS/src/ProblemVec.cc           | 57 ++++++++++++----------
 AMDiS/src/ProblemVec.h            |  6 ++-
 AMDiS/src/Quadrature.h            | 74 +++++++++++++---------------
 AMDiS/src/ResidualEstimator.cc    | 80 +++++++++++++++++--------------
 AMDiS/src/SolutionDataStorage.h   |  4 +-
 AMDiS/src/SubAssembler.cc         |  4 +-
 13 files changed, 162 insertions(+), 140 deletions(-)

diff --git a/AMDiS/src/CoarseningManager.cc b/AMDiS/src/CoarseningManager.cc
index 77eb26db..3a102a2b 100644
--- a/AMDiS/src/CoarseningManager.cc
+++ b/AMDiS/src/CoarseningManager.cc
@@ -107,11 +107,10 @@ namespace AMDiS {
     do {
       doMore = false;
       el_info = stack->traverseFirst(mesh, -1, flag);
-      while (el_info)
-	{
-	  coarsenFunction(el_info);
-	  el_info = stack->traverseNext(el_info);
-	}
+      while (el_info) {
+	coarsenFunction(el_info);
+	el_info = stack->traverseNext(el_info);
+      }
     } while (doMore);
 
     DELETE stack;
diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h
index b210b7f4..0656bdb9 100644
--- a/AMDiS/src/ComponentTraverseInfo.h
+++ b/AMDiS/src/ComponentTraverseInfo.h
@@ -70,7 +70,11 @@ namespace AMDiS {
     }
     
     const FiniteElemSpace *getAuxFESpace(int i) {
-      return ((i <= static_cast<int>(auxFESpaces.size())) ? auxFESpaces[i] : NULL);
+      return ((i < static_cast<int>(auxFESpaces.size())) ? auxFESpaces[i] : NULL);
+    }
+
+    void setAuxFESpace(const FiniteElemSpace* fe, int pos) {
+      auxFESpaces[pos] = fe;
     }
 
     int getStatus() {
@@ -164,6 +168,22 @@ namespace AMDiS {
       return vectorComponents[row].getAuxFESpace(0);
     }
 
+    void fakeFESpace(const FiniteElemSpace *feOld, 
+		     const FiniteElemSpace *feNew)
+    {
+      for (int i = 0; i < nComponents; i++) {
+	for (int j = 0; j < nComponents; j++) {
+	  if (matrixComponents[i][j].getAuxFESpace(0) == feOld) {
+	    matrixComponents[i][j].setAuxFESpace(feNew, 0);
+	  }
+	}
+
+	if (vectorComponents[i].getAuxFESpace(0) == feOld) {
+	  vectorComponents[i].setAuxFESpace(feNew, 0);
+	}
+      }
+    }
+
   protected:
     int nComponents;
 
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 0ad62856..95878057 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -441,11 +441,11 @@ namespace AMDiS {
     std::vector<Operator*>::iterator it = operators.begin();
     std::vector<double*>::iterator factorIt = operatorFactor.begin();
     for (; it != operators.end(); ++it, ++factorIt) {
-	if ((*it)->getNeedDualTraverse() == false) {
-	    (*it)->getElementMatrix(elInfo, 
-				    elementMatrix, 
-				    *factorIt ? **factorIt : 1.0);
-	}
+      if ((*it)->getNeedDualTraverse() == false) {
+	  (*it)->getElementMatrix(elInfo, 
+				  elementMatrix, 
+				  *factorIt ? **factorIt : 1.0);
+      }
     }        
 
     addElementMatrix(factor, *elementMatrix, bound);   
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index b2da7a69..d6742ffe 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -16,6 +16,12 @@ namespace AMDiS {
       (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
       break;
     case COARSE_INTERPOL:
+      if (feSpace == NULL || !feSpace) {
+	ERROR_EXIT("ERR1\n");
+      }
+      if (feSpace->getBasisFcts() == NULL || !(feSpace->getBasisFcts())) {
+	ERROR_EXIT("ERR2\n");
+      }
       (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
       break;
     default:
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 9b337807..9fede244 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -498,7 +498,7 @@ namespace AMDiS {
      */
     inline double H1Norm(Quadrature* q = NULL) const {
       return sqrt(H1NormSquare());
-    }; 
+    }
 
     /** \brief
      * Calculates square of H1 norm of this DOFVector
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 9bf1874b..bf86d9e0 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -744,6 +744,7 @@ namespace AMDiS {
   DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
   {
     feSpace = rhs.feSpace;
+    DOFVectorBase<T>::feSpace = rhs.feSpace;
     this->nBasFcts = rhs.nBasFcts;
     vec = rhs.vec;
     this->elementVector = new ElementVector(this->nBasFcts);
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 998a7b6c..8b1c1f38 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -1575,20 +1575,14 @@ namespace AMDiS {
 		AbstractFunction<double, double> *af,
 		WorldVector<double> *wv);
 
-    /** \brief
-     * Implementation of \ref OperatorTerm::initElement().
-     */
+    /// Implementation of \ref OperatorTerm::initElement().
     void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
 		     Quadrature *quad = NULL);
 
-    /** \brief
-     * Implements FirstOrderTerm::getLb().
-     */
+    /// Implements FirstOrderTerm::getLb().
     void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
 
-    /** \brief
-     * Implements FirstOrderTerm::eval().
-     */
+    /// Implements FirstOrderTerm::eval().
     void eval(int nPoints,
 	      const double *uhAtQP,
 	      const WorldVector<double> *grdUhAtQP,
@@ -1597,23 +1591,16 @@ namespace AMDiS {
 	      double factor) const;
 
   protected:
-    /** \brief
-     * DOFVector to be evaluated at quadrature points.
-     */
+    /// DOFVector to be evaluated at quadrature points.
     DOFVectorBase<double>* vec;
 
-    /** \brief
-     * v at quadrature points.
-     */
+    /// v at quadrature points.
     double *vecAtQPs;
 
-    /** \brief
-     * Function f.
-     */
+    /// Function f.
     AbstractFunction<double, double> *f;
 
-    /**
-     */
+    ///
     WorldVector<double> *b;
   };
 
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 0a532940..7887a19b 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -729,7 +729,6 @@ namespace AMDiS {
 
     for (int i = 0; i < nComponents; i++) {
       for (int j = 0; j < nComponents; j++) {
-
 	// Only if this variable is true, the current matrix will be assembled.	
 	bool assembleMatrix = true;
 	// The DOFMatrix which should be assembled (or not, if assembleMatrix
@@ -780,18 +779,18 @@ namespace AMDiS {
 	  } else if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
 
 	    // Row fe space and col fe space are both equal, but right hand side has at
-	    // least one another aux fe space. In this case, do not assemble the rhs,
-	    // this will be done afterwards.
+	    // least one another aux fe space. 
 
 	    assembleOnOneMesh(componentSpaces[i],
 			      assembleFlag,
 			      assembleMatrix ? matrix : NULL,
 			      ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
 
-	    assembleOnDifMeshes2(componentSpaces[1], componentSpaces[0],
+	    assembleOnDifMeshes2(componentSpaces[i], 
+				 traverseInfo.getAuxFESpace(i, j),
 				 assembleFlag,
 				 NULL,
-				 rhs->getDOFVector(1));
+				 ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
 
 	  } else {
 	    ERROR_EXIT("Possible? If yes, not yet implemented!\n");
@@ -803,12 +802,12 @@ namespace AMDiS {
 			    assembleFlag,
 			    assembleMatrix ? matrix : NULL,
 			    ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
-
+	  
 	  assembleOnDifMeshes2(componentSpaces[i],
 			       traverseInfo.getAuxFESpace(i, j),
 			       assembleFlag,
 			       assembleMatrix ? matrix : NULL,
-			       NULL);
+			       ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
 
 	} else if (traverseInfo.getStatus(i, j) ==  SingleComponentInfo::DIF_SPACES_NO_AUX ||
 		   traverseInfo.getStatus(i, j) ==  SingleComponentInfo::DIF_SPACES_WITH_AUX) {
@@ -888,7 +887,8 @@ namespace AMDiS {
   void ProblemVec::addMatrixOperator(Operator *op, 
 				     int i, int j,
 				     double *factor,
-				     double *estFactor)
+				     double *estFactor,
+				     bool fake)
   {
     FUNCNAME("ProblemVec::addMatrixOperator()");
    
@@ -900,34 +900,41 @@ namespace AMDiS {
 	setBoundaryConditionMap((*systemMatrix)[i][i]->getBoundaryManager()->
 				getBoundaryConditionMap());
     }    
-    (*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
 
-    traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces()); 
+    (*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
 
-    for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
-      if ((op->getAuxFESpaces())[k] != componentSpaces[i] ||
-	  (op->getAuxFESpaces())[k] != componentSpaces[j]) {
-	op->setNeedDualTraverse(true);
-	break;
-      }
-    }    
+    if (!fake) {    
+      traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces()); 
+      
+      for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
+	if ((op->getAuxFESpaces())[k]->getMesh() != componentSpaces[i]->getMesh() ||
+	    (op->getAuxFESpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) {
+	  op->setNeedDualTraverse(true);
+	  break;
+	}
+      }    
+    } 
   }
 
   void ProblemVec::addVectorOperator(Operator *op, int i,
 				     double *factor,
-				     double *estFactor)
+				     double *estFactor,
+				     bool fake)
   {
     FUNCNAME("ProblemVec::addVectorOperator()");
 
     rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
-    traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces()); 
 
-    for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
-      if ((op->getAuxFESpaces())[j] != componentSpaces[i]) {
-	op->setNeedDualTraverse(true);
-	break;
-      }
-    }    
+    if (!fake) {
+      traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces()); 
+      
+      for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
+	if ((op->getAuxFESpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) {
+	  op->setNeedDualTraverse(true);
+	  break;
+	}
+      }    
+    }
   }
 
   void ProblemVec::addDirichletBC(BoundaryType type, int system,
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index 2d09bab7..37eb1985 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -216,12 +216,14 @@ namespace AMDiS {
     /// Adds an Operator to \ref A.
     void addMatrixOperator(Operator *op, int i, int j,
 			   double *factor = NULL,
-			   double *estFactor = NULL);
+			   double *estFactor = NULL,
+			   bool fake = false);
 
     /// Adds an Operator to \ref rhs.
     void addVectorOperator(Operator *op, int i,
 			   double *factor = NULL,
-			   double *estFactor = NULL);
+			   double *estFactor = NULL,
+			   bool fake = false);
 
     /// Adds dirichlet boundary conditions.
     virtual void addDirichletBC(BoundaryType type, int system,
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index 0f965185..0b34eed7 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -84,19 +84,13 @@ namespace AMDiS {
     {}
 
   public:
-    /** \brief
-     * Copy constructor
-     */
+    /// Copy constructor
     Quadrature(const Quadrature&);
 
-    /** \brief
-     * Destructor
-     */
+    /// Destructor
     ~Quadrature();
 
-    /** \brief
-     * Returns a Quadrature for dimension dim exact for degree degree.
-     */
+    /// Returns a Quadrature for dimension dim exact for degree degree.
     static Quadrature *provideQuadrature(int dim, int degree);
 
     /** \brief
@@ -116,42 +110,40 @@ namespace AMDiS {
      */
     inline const std::string& getName() { 
       return name; 
-    };
+    }
 
-    /** \brief
-     * Returns \ref n_points
-     */
+    /// Returns \ref n_points
     inline int getNumPoints() const {
       return n_points;
-    };
+    }
 
     /** \brief
      * Returns \ref w[p]
      */
     inline double getWeight(int p) const {
       return w[p];
-    };
+    }
 
     /** \brief
      * Returns \ref w.
      */
     inline double* getWeight() const { 
       return w; 
-    };
+    }
 
     /** \brief
      * Returns \ref dim
      */
     inline int getDim() const { 
       return dim; 
-    };
+    }
 
     /** \brief
      * Returns \ref degree
      */
     inline int getDegree() const { 
       return degree; 
-    };
+    }
 
     /** \brief
      * Returns a pointer to a vector storing the values of a doubled valued 
@@ -184,7 +176,7 @@ namespace AMDiS {
      */
     inline double getLambda(int a, int b) const {
       return (lambda ? (*lambda)[a][b] : 0.0);
-    };
+    }
 
     /** \brief
      * Returns \ref lambda[a] which is a DimVec<double> containing the 
@@ -192,14 +184,14 @@ namespace AMDiS {
      */
     inline const DimVec<double>& getLambda(int a) const {
       return (*lambda)[a]; 
-    };
+    }
 
     /** \brief
      * Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
      */
     VectorOfFixVecs<DimVec<double> > *getLambda() const { 
       return lambda; 
-    };
+    }
 
 
   public:
@@ -224,9 +216,7 @@ namespace AMDiS {
      */
     int dim;
 
-    /** \brief
-     * Number of quadrature points
-     */
+    /// Number of quadrature points
     int n_points;
 
     /** \brief
@@ -378,7 +368,7 @@ namespace AMDiS {
 	D2Phi(NULL), 
 	quadrature(quad), 
 	basisFunctions(basFcts) 
-    {};
+    {}
 
     /** \brief
      * Copy constructor
@@ -424,21 +414,21 @@ namespace AMDiS {
 
       ERROR_EXIT("invalid flag\n");
       return false;
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature
      */
     inline const Quadrature* getQuadrature() const { 
       return quadrature; 
-    };
+    }
 
     /** \brief
      * Returns \ref max_points
      */
     inline int getMaxQuadPoints() { 
       return max_points; 
-    };
+    }
 
     /** \brief
      * Returns (*\ref D2Phi)[q][i][j][m]
@@ -455,63 +445,63 @@ namespace AMDiS {
      */
     inline const double getGradient(int q, int i ,int j) const {
       return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
-    };
+    }
 
     /** \brief
      * Returns (*\ref grdPhi)[q]
      */
     inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const {
       return (grdPhi) ? &((*grdPhi)[q]) : NULL;
-    };
+    }
 
     /** \brief
      * Returns \ref phi[q][i]
      */
     inline const double getPhi(int q, int i) const {
       return (phi) ? phi[q][i] : 0.0;
-    };
+    }
 
     /** \brief
      * Returns \ref phi[q]
      */
     inline const double *getPhi(int q) const {
       return (phi) ? phi[q] : NULL;
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->integrateStdSimplex(f)
      */
     inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) {
       return quadrature->integrateStdSimplex(f); 
-    };
+    }
   
     /** \brief
      * Returns \ref quadrature ->getNumPoints()
      */
     inline int getNumPoints() const { 
       return quadrature->getNumPoints();
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->getWeight(p)
      */
     inline double getWeight(int p) const {
       return quadrature->getWeight(p);
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->getDim()
      */
     inline int getDim() const { 
       return quadrature->getDim(); 
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->getDegree()
      */
     inline int getDegree() const { 
       return quadrature->getDegree();
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->grdFAtQp(f, vec)
@@ -522,7 +512,7 @@ namespace AMDiS {
 	      WorldVector<double>* vec) const 
     { 
       return quadrature->grdFAtQp(f, vec);
-    };
+    }
   
     /** \brief
      * Returns \ref quadrature ->fAtQp(f, vec)
@@ -531,28 +521,28 @@ namespace AMDiS {
 			       DimVec<double> >& f, double *vec) const 
     {
       return quadrature->fAtQp(f, vec);    
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->getLambda(a,b)
      */
     inline double getLambda(int a,int b) const {
       return quadrature->getLambda(a,b);
-    };
+    }
 
     /** \brief
      * Returns \ref quadrature ->getLambda(a)
      */
     inline const DimVec<double>& getLambda(int a) const { 
       return quadrature->getLambda(a); 
-    };
+    }
 
     /** \brief
      * Returns \ref basisFunctions
      */
     inline BasisFunction* getBasisFunctions() const { 
       return basisFunctions; 
-    };
+    }
 
   protected:
     /** \brief
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index fb6e44d1..f820c4a9 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -178,13 +178,14 @@ namespace AMDiS {
 
 #if 1
   void ResidualEstimator::estimateElement(ElInfo *elInfo)
-  {
+  {    
     FUNCNAME("ResidualEstimator::estimateElement()");
 
     TEST_EXIT_DBG(nSystems > 0)("no system set\n");
 
     double val = 0.0;
     std::vector<Operator*>::iterator it;
+    std::vector<double*>::iterator itfac;
     Element *el = elInfo->getElement();
     double det = elInfo->getDet();
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
@@ -201,17 +202,20 @@ namespace AMDiS {
 	continue;
 
       // init assemblers
-      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
+      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
+	   itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
 	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
-	   ++it) {
-	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+	   ++it, ++itfac) {
+	if (**itfac != 0.0) {
+	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+	}
       }
 
       if (C0) {
 	for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
 	     it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd(); 
 	     ++it) {
-	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	
 	}
       }
 
@@ -238,20 +242,23 @@ namespace AMDiS {
       }
            
       if (C0) {  
-	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(); 
+	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
+	     itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
 	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
-	     ++it) {
-	  if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
-	    uhQP = GET_MEMORY(double, nPoints);
-	    uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
-	  }
-	  if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
-	    grdUh_qp = new WorldVector<double>[nPoints];
-	    uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
-	  }
-	  if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
-	    D2uhqp = new WorldMatrix<double>[nPoints];
-	    uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
+	     ++it, ++itfac) {
+	  if (**itfac != 0.0) {
+	    if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
+	      uhQP = GET_MEMORY(double, nPoints);
+	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
+	    }
+	    if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
+	      grdUh_qp = new WorldVector<double>[nPoints];
+	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
+	    }
+	    if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
+	      D2uhqp = new WorldMatrix<double>[nPoints];
+	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
+	    }
 	  }
 	}
 	
@@ -396,24 +403,27 @@ namespace AMDiS {
 		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
 		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
 		 ++it, ++fac) {
-	      for (int iq = 0; iq < nPointsSurface_; iq++) {
-		localJump_[iq].set(0.0);
-	      }
-		  
-	      (*it)->weakEvalSecondOrder(nPointsSurface_,
-					 grdUhEl_.getValArray(),
-					 localJump_.getValArray());
-	      double factor = *fac ? **fac : 1.0;
-	      if (factor != 1.0) {
+
+	      if (**fac != 0.0) {
+		for (int iq = 0; iq < nPointsSurface_; iq++) {
+		  localJump_[iq].set(0.0);
+		}
+		
+		(*it)->weakEvalSecondOrder(nPointsSurface_,
+					   grdUhEl_.getValArray(),
+					   localJump_.getValArray());
+		double factor = *fac ? **fac : 1.0;
+		if (factor != 1.0) {
+		  for (int i = 0; i < nPointsSurface_; i++) {
+		    localJump_[i] *= factor;
+		  }
+		}
+		
 		for (int i = 0; i < nPointsSurface_; i++) {
-		  localJump_[i] *= factor;
+		  jump_[i] += localJump_[i];
 		}
-	      }
-		  
-	      for (int i = 0; i < nPointsSurface_; i++) {
-		jump_[i] += localJump_[i];
-	      }
-	    }				     
+	      }		
+	    }
 	  }
 	      
 	  val = 0.0;
@@ -454,7 +464,7 @@ namespace AMDiS {
     est_sum += est_el;
     est_max = max(est_max, est_el);
 
-    elInfo->getElement()->setMark(0);   
+    elInfo->getElement()->setMark(0);  
   }
 
 
diff --git a/AMDiS/src/SolutionDataStorage.h b/AMDiS/src/SolutionDataStorage.h
index 8fab0084..83ebb55f 100644
--- a/AMDiS/src/SolutionDataStorage.h
+++ b/AMDiS/src/SolutionDataStorage.h
@@ -64,8 +64,8 @@ namespace AMDiS {
     ~SolutionData() 
     {
       if (!serialized) {
-	DELETE solution;
-	deleteFeSpace(feSpace);
+	//DELETE solution;
+	//deleteFeSpace(feSpace);
       }
     }
 
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index 0a40ba8f..f8dadcc0 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -156,11 +156,11 @@ namespace AMDiS {
 
     TEST_EXIT(elInfo->getMesh() == vec->getFESpace()->getMesh())("Vector and fe space do not fit together!\n");
 
+    Quadrature *localQuad = quad ? quad : quadrature;
+
     if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) 
       return valuesAtQPs[vec]->values.getValArray();
 
-    Quadrature *localQuad = quad ? quad : quadrature;
-
     if (!valuesAtQPs[vec]) {
       valuesAtQPs[vec] = new ValuesAtQPs;
     }
-- 
GitLab