some modifications in the base_problems

extensions/Helpers.h

@@ -83,7 +83,7 @@ namespace Helpers {
   inline int signum(int val) { return val > 0 ? 1 : (val < 0 ? -1 : 0); }
   inline int signum(bool val) { return val ? 1 : -1; }
   inline int toInt(bool val) { return val ? 1 : 0; }
-  inline bool isNumber(double val) { return !boost::math::isnan(val) && !boost::math::isinf(val); }
+
   inline bool toBool(double val, double tol = DBL_TOL) { return val > tol || val < -tol ? true : false; }
   inline bool isPositiv(double val, double tol = DBL_TOL) { return val > -tol; }
   inline bool isStrictPositiv(double val, double tol = DBL_TOL) { return val > tol; }

extensions/Refinement.h

@@ -283,6 +283,6 @@ protected:
 };
 
 #include "Refinement_Level.h"
-#include "Refinement_MeshSize.h"
+// #include "Refinement_MeshSize.h"
 
 #endif // EXTENSIONS_REFINEMENT_H

extensions/Refinement_Level.h

@@ -217,8 +217,8 @@ public:
         center += coords[i];
       }
       center *= 1.0 / static_cast<double>(coords.size());
-
-      int refineLevel = (*refineFct)(std::make_pair<WorldVector<double>, double>(center, meanValue));
+      std::pair<WorldVector<double>, double> swap(center, meanValue);
+      int refineLevel = (*refineFct)(swap);
       int oldLevel = elInfo->getLevel();
       elInfo->getElement()->setMark( calcMark(refineLevel, oldLevel) );
 

extensions/Refinement_MeshSize.h

@@ -166,7 +166,8 @@ public:
       }
       center *= 1.0 / static_cast<double>(coords.size());
 
-      double refineH = (*refineFct)(std::make_pair<WorldVector<double>, double>(center, meanValue));
+      std::pair<WorldVector<double>, double> swap(center, meanValue);
+      int refineH = (*refineFct)(swap);
       double oldH = calcMeshSize(elInfo);
       elInfo->getElement()->setMark( calcMark(refineH, oldH) );
 

extensions/Views.h

@@ -433,21 +433,21 @@ T evalAtPoint(const DOFVector<T> &obj,
 
 // DOFView can be accessed by own method
 template<typename T >
-inline T evalAtPoint(const DOFView<T> &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
+T evalAtPoint(const DOFView<T> &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
 {
   return obj.evalAtPoint(p, elInfo);
 }
 
 // intrinsic types can return values directly
 template<typename T>
-inline T evalAtPoint(T &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
+typename boost::enable_if< boost::is_pod<T>, T >::type evalAtPoint(const T &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
 {
   return obj;
 }
 
 // AbstractFunctions with argument-type = WorldVector can be accessed directly
 template<typename T>
-inline T evalAtPoint(AbstractFunction<T, WorldVector<double> > &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
+T evalAtPoint(const AbstractFunction<T, WorldVector<double> > &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
 {
   return obj(p);
 }

extensions/base_problems/BaseProblem.h

@@ -24,6 +24,12 @@
 #include "VtuReader.h"
 #include "pugixml.hpp"
 
+#if __cplusplus > 199711L
+  #define _OVERRIDE_ override
+#else
+  #define _OVERRIDE_
+#endif
+
 using namespace AMDiS;
 
 const Flag MESH_ADOPTED = 1<<2;
@@ -40,6 +46,14 @@ public:
   {
     if (prob)
       delete prob;
+    
+    for (size_t i = 0; i < operatorTermList.size(); i++)
+      delete operatorTermList[i];
+    operatorTermList.clear();
+    
+    for (size_t i = 0; i < operatorList.size(); i++)
+      delete operatorList[i];
+    operatorList.clear();
   }
   
   /// Initialisation of the problem.
@@ -51,12 +65,16 @@ public:
   /// is called in \ref initTimeInteface after feSpace and mesh are initialized
   virtual void initData() {};
   
+  /// Method is called at the end of \ref initTimeInteface
+  virtual void finalizeData() {};
+  
   /// calls \ref initData, \ref fillOperators and \ref fillBoundaryConditions in this ordering
   virtual void initTimeInterface()
   { 
     initData();
     fillOperators();
     fillBoundaryConditions();
+    finalizeData();
   }
   
   /// read solution DOFVectors from .arh, .dat or .vtu files
@@ -204,6 +222,10 @@ protected:
   int oldMeshChangeIdx;
 
   double oldTimestep;
+  
+  /// lists of objects to delete in the destructor
+  std::vector<OperatorTerm*> operatorTermList;
+  std::vector<Operator*> operatorList;
 
 private:
   std::vector<std::pair<int,int> > fixedMatrixTimestep;

extensions/base_problems/BaseProblem_RB.h

@@ -42,10 +42,11 @@ public:
   
   Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) 
   {
-    prob->oneIteration(adaptInfo, toDo);
+    Flag flag = prob->oneIteration(adaptInfo, toDo);
     
     oldTimestep = adaptInfo->getTimestep();
     oldMeshChangeIdx = getMesh()->getChangeIndex();
+    return flag;
   }
   
   Flag stageIteration(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix, bool asmVector)

extensions/base_problems/CahnHilliard.cc

@@ -34,12 +34,12 @@ CahnHilliard::CahnHilliard(const std::string &name_) :
   minusEpsSqr(-0.01)
 {
   // parameters for CH
-  Parameters::get(name + "->use mobility", useMobility); // mobility
-  Parameters::get(name + "->gamma", gamma); // mobility
-  Parameters::get(name + "->epsilon", eps); // interface width
+  Parameters::get(name_ + "->use mobility", useMobility); // mobility
+  Parameters::get(name_ + "->gamma", gamma); // mobility
+  Parameters::get(name_ + "->epsilon", eps); // interface width
   
   // type of double well: 0= [0,1], 1= [-1,1]
-  Parameters::get(name + "->double-well type", doubleWell); 
+  Parameters::get(name_ + "->double-well type", doubleWell); 
 
   // transformation of the parameters
   minusEps = -eps;
@@ -50,10 +50,6 @@ CahnHilliard::CahnHilliard(const std::string &name_) :
 }
 
 
-CahnHilliard::~CahnHilliard() 
-{ FUNCNAME("CahnHilliard::~CahnHilliard()"); }
-
-
 void CahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo) 
 { FUNCNAME("CahnHilliard::solveInitialProblem()");
 
@@ -123,21 +119,21 @@ void CahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo)
 
 void CahnHilliard::fillOperators()
 { FUNCNAME("CahnHilliard::fillOperators()");
-  MSG("CahnHilliard::fillOperators()\n");
   
-  int degree = prob->getFeSpace(0)->getBasisFcts()->getDegree();
+  const FiniteElemSpace* feSpace = prob->getFeSpace();
+  int degree = feSpace->getBasisFcts()->getDegree();
 
   // c
-  Operator *opChMnew = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+  Operator *opChMnew = new Operator(feSpace, feSpace);
   opChMnew->addZeroOrderTerm(new Simple_ZOT);
-  Operator *opChMold = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+  Operator *opChMold = new Operator(feSpace, feSpace);
   opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0)));
   // -nabla*(grad(c))
-  Operator *opChL = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  Operator *opChL = new Operator(feSpace, feSpace);
   opChL->addSecondOrderTerm(new Simple_SOT);
   
   // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
-  Operator *opChLM = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+  Operator *opChLM = new Operator(feSpace, feSpace);
   if (useMobility) {
     if (doubleWell == 0)
       opChLM->addSecondOrderTerm(new VecAtQP_SOT(
@@ -151,7 +147,7 @@ void CahnHilliard::fillOperators()
     opChLM->addSecondOrderTerm(new Simple_SOT(gamma));
   
   // -2*c_old^3 + 3/2*c_old^2
-  Operator *opChMPowExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  Operator *opChMPowExpl = new Operator(feSpace, feSpace);
   opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
     prob->getSolution()->getDOFVector(0),
     new AMDiS::Pow<3>(-2.0, 3*degree)));
@@ -162,7 +158,7 @@ void CahnHilliard::fillOperators()
   }
 
   // -3*c_old^2 * c
-  Operator *opChMPowImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  Operator *opChMPowImpl = new Operator(feSpace, feSpace);
   opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
     prob->getSolution()->getDOFVector(0),
     new AMDiS::Pow<2>(-3.0, 2*degree)));
@@ -183,7 +179,6 @@ void CahnHilliard::fillOperators()
   // . . . vectorOperators . . . . . . . . . . . . . . .
   prob->addVectorOperator(*opChMPowExpl,0);            /// < -2*phi*c_old^3 , psi >
 
-//   setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
   
   // dt(c) = laplace(mu) - u*grad(c)
   // -----------------------------------
@@ -192,5 +187,15 @@ void CahnHilliard::fillOperators()
   // . . . vectorOperators . . . . . . . . . . . . . . .
   prob->addVectorOperator(*opChMold,1, getInvTau());   /// < phi*c^old/tau , psi >
   
-//   setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
 }
+
+
+void CahnHilliard::finalizeData()
+{ FUNCNAME("CahnHilliard::finalizeData()");
+  
+  setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
+  setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
+  if (!useMobility)
+    setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
+}
+

extensions/base_problems/CahnHilliard.h

@@ -20,30 +20,30 @@
 
 #include "AMDiS.h"
 #include "BaseProblem.h"
-#include "ExtendedProblemStat.h"
 #include "chns.h"
 
 using namespace AMDiS;
-
-class CahnHilliard : public BaseProblem<ExtendedProblemStat>
+    
+class CahnHilliard : public BaseProblem<ProblemStat>
 {
 public: // definition of types
 
-  typedef BaseProblem<ExtendedProblemStat> super;
+  typedef BaseProblem<ProblemStat> super;
 
 public: // public methods
 
   CahnHilliard(const std::string &name_);
-  ~CahnHilliard();
+  ~CahnHilliard() {};
 
-  virtual void initData() {};
-  virtual void solveInitialProblem(AdaptInfo *adaptInfo);
+  void solveInitialProblem(AdaptInfo *adaptInfo) _OVERRIDE_;
 
   double getEpsilon() { return eps; }
   int getDoubleWellType() { return doubleWell; }
 
-  virtual void fillOperators();
-  virtual void fillBoundaryConditions() {}
+  void fillOperators() _OVERRIDE_;
+  void fillBoundaryConditions() _OVERRIDE_ {}
+  
+  void finalizeData() _OVERRIDE_;
 
 protected: // protected variables
 
@@ -61,5 +61,7 @@ protected: // protected variables
   double epsSqr;
   double minusEpsSqr;
 };
+  
+
 
 #endif // CAHN_HILLIARD_H

extensions/base_problems/NavierStokes_TH_MultiPhase.cc

@@ -45,7 +45,6 @@ NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_,
 
 NavierStokes_TH_MultiPhase::~NavierStokes_TH_MultiPhase()
 {
-  if (multiPhase) { delete multiPhase; multiPhase = NULL; }
   if (densityPhase) { delete densityPhase; densityPhase = NULL; }
   if (viscosityPhase) { delete viscosityPhase; viscosityPhase = NULL; }
 }
@@ -193,4 +192,3 @@ void NavierStokes_TH_MultiPhase::addLaplaceTerm(int i)
   }
   prob->addMatrixOperator(*opLaplaceUi, i, i);
 }
-

extensions/base_problems/NavierStokes_TH_MultiPhase.h

@@ -22,13 +22,13 @@
 using namespace AMDiS;
 
 /** \ingroup NavierStokes_TaylorHood
- * \brief
- * Navier-Stokes multi-phase problem, using Taylor Hood elements
- * Standard implementation for two phases, but can be extended to
- * many more, by defining the multiPhase variable and overloading
- * the initTimestep, where densityPhase and viscosityPhase are defined
- * depending on the phases
- */
+* \brief
+* Navier-Stokes multi-phase problem, using Taylor Hood elements
+* Standard implementation for two phases, but can be extended to
+* many more, by defining the multiPhase variable and overloading
+* the initTimestep, where densityPhase and viscosityPhase are defined
+* depending on the phases
+*/
 class NavierStokes_TH_MultiPhase : public NavierStokes_TaylorHood
 {
 public: // typedefs
@@ -43,25 +43,25 @@ public: // methods
   ~NavierStokes_TH_MultiPhase();
 
   /** definition of constant multiPhase and const
-   * viscosity/density DOFVectors,
-   * initTimeInterface of super
-   **/
-  virtual void initData();
+  * viscosity/density DOFVectors,
+  * initTimeInterface of super
+  **/
+  void initData() _OVERRIDE_;
 
   /** initTimestep of super and
-   * initialization of densityPhase and viscosityPhase
-   **/
-  virtual void initTimestep(AdaptInfo *adaptInfo);
+  * initialization of densityPhase and viscosityPhase
+  **/
+  void initTimestep(AdaptInfo *adaptInfo) _OVERRIDE_;
 
   /** pointer to the multiPhase DOFVector, that
-   * indicates the different phases. Will be initialized
-   * in \ref initTimeInterface with const 1
-   **/
+  * indicates the different phases. Will be initialized
+  * in \ref initTimeInterface with const 1
+  **/
   void setMultiPhase(DOFVector<double>* p) { multiPhase=p; }
 
-  virtual void fillOperators();
+  void fillOperators() _OVERRIDE_;
 
-  virtual void addLaplaceTerm(int i);
+  void addLaplaceTerm(int i) _OVERRIDE_;
   
   // get-methods
   

extensions/base_problems/NavierStokes_TH_MultiPhase_RB.h

@@ -47,12 +47,12 @@ public: // methods
    * viscosity/density DOFVectors,
    * initTimeInterface of super
    **/
-  virtual void initData();
+  void initData() override;
 
   /** initTimestep of super and
    * initialization of densityPhase and viscosityPhase
    **/
-  virtual void initTimestep(AdaptInfo *adaptInfo);
+  void initTimestep(AdaptInfo *adaptInfo) override;
 
   /** pointer to the multiPhase DOFVector, that
    * indicates the different phases. Will be initialized
@@ -60,8 +60,8 @@ public: // methods
    **/
   DOFVector<double> *getMultiPhase() { return multiPhase; }
 
-  virtual void fillOperators();
-  virtual void addLaplaceTerm(int i);
+  void fillOperators() override;
+  void addLaplaceTerm(int i) override;
   
   // get-methods
   

extensions/base_problems/NavierStokes_TaylorHood.cc

@@ -33,33 +33,33 @@ NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_, bool
   fileWriter(NULL)
 {  
   // fluid viscosity
-  Initfile::get(name + "->viscosity", viscosity);
+  Initfile::get(self::name + "->viscosity", viscosity);
   // fluid density
-  Initfile::get(name + "->density", density);
+  Initfile::get(self::name + "->density", density);
 
   // type of laplace operator: 
   //   0... div(nu*grad(u)), 
   //   1... div(0.5*nu*(grad(u)+grad(u)^T))
-  Initfile::get(name + "->laplace operator", laplaceType); 
+  Initfile::get(self::name + "->laplace operator", laplaceType); 
 
   // type of non-linear term: 
   //   0... u^old*grad(u_i^old), 
   //   1... u'*grad(u_i^old), 
   //   2... u^old*grad(u'_i)
   //   3... (1)+(2)-(0)
-  Initfile::get(name + "->non-linear term", nonLinTerm); 
+  Initfile::get(self::name + "->non-linear term", nonLinTerm); 
 
   // theta-scheme parameter
-  Initfile::get(name + "->theta", theta); 
+  Initfile::get(self::name + "->theta", theta); 
   theta1 = 1.0-theta;
   minusTheta1 = -theta1;
 
   // volume force
   force.set(0.0);
-  Initfile::get(name + "->force", force);
+  Initfile::get(self::name + "->force", force);
 
-  oldSolution.resize(dow);
-  for (size_t i = 0; i < dow; i++)
+  oldSolution.resize(self::dow);
+  for (size_t i = 0; i < self::dow; i++)
     oldSolution[i] = NULL;
 }
 
@@ -72,7 +72,7 @@ NavierStokes_TaylorHood::~NavierStokes_TaylorHood()
     velocity = NULL;
   }
   
-  for (size_t i = 0; i < dow; i++) {
+  for (size_t i = 0; i < self::dow; i++) {
     if (oldSolution[i] != NULL)
       delete oldSolution[i];
     oldSolution[i] = NULL;
@@ -89,12 +89,12 @@ void NavierStokes_TaylorHood::initData()
 { FUNCNAME("NavierStokes_TaylorHood::initTimeInterface()");
 
   if (velocity == NULL)
-    velocity = new DOFVector<WorldVector<double> >(getFeSpace(0), "velocity");
+    velocity = new DOFVector<WorldVector<double> >(self::getFeSpace(0), "velocity");
   
-  for (size_t i = 0; i < dow; i++)
-    oldSolution[i] = new DOFVector<double>(getFeSpace(i), "old(v_"+Helpers::toString(i)+")");
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i] = new DOFVector<double>(self::getFeSpace(i), "old(v_"+Helpers::toString(i)+")");
   
-  fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
+  fileWriter = new FileVectorWriter(name + "->velocity->output", self::getFeSpace()->getMesh(), velocity);
 
   super::initData();
 }
@@ -104,8 +104,8 @@ void NavierStokes_TaylorHood::transferInitialSolution(AdaptInfo *adaptInfo)
 { FUNCNAME("NavierStokes_TaylorHood::transferInitialSolution()");
 
   calcVelocity();
-  for (size_t i = 0; i < dow; i++)
-    oldSolution[i]->copy(*prob->getSolution()->getDOFVector(i));
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i]->copy(*self::prob->getSolution()->getDOFVector(i));
   
   super::transferInitialSolution(adaptInfo);
 }
@@ -116,48 +116,51 @@ void NavierStokes_TaylorHood::fillOperators()
 
   WorldVector<DOFVector<double>* > vel;
   for (size_t k = 0; k < dow; k++)
-    vel[k] = prob->getSolution()->getDOFVector(k);
+    vel[k] = self::prob->getSolution()->getDOFVector(k);
+  
+  const FiniteElemSpace* feSpace_u = self::getFeSpace(0);
+  const FiniteElemSpace* feSpace_p = self::getFeSpace(self::dow);
   
   // fill operators for prob
-  for (size_t i = 0; i < dow; ++i) {
+  for (size_t i = 0; i < self::dow; ++i) {
     /// < (1/tau)*u'_i , psi >
-    Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
+    Operator *opTime = new Operator(feSpace_u, feSpace_u);
     opTime->addTerm(new Simple_ZOT(density));
-    prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
+    self::prob->addMatrixOperator(*opTime, i, i, self::getInvTau(), self::getInvTau());
     
     /// < (1/tau)*u_i^old , psi >
-    Operator *opTimeOld = new Operator(getFeSpace(i), getFeSpace(i));
+    Operator *opTimeOld = new Operator(feSpace_u, feSpace_u);
     opTimeOld->addTerm(new Phase_ZOT(getOldSolution(i), density));
-    prob->addVectorOperator(*opTimeOld, i, getInvTau(), getInvTau());
+    self::prob->addVectorOperator(*opTimeOld, i, self::getInvTau(), self::getInvTau());
 
     /// < u^old*grad(u_i^old) , psi >
     if (nonLinTerm == 0 || nonLinTerm == 3) {
-      Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
+      Operator *opUGradU0 = new Operator(feSpace_u, feSpace_u);
       opUGradU0->addTerm(new WorldVec_FOT(vel, density), GRD_PHI);
       opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
       if (nonLinTerm == 0)
-        prob->addVectorOperator(*opUGradU0, i);
+        self::prob->addVectorOperator(*opUGradU0, i);
       else
-        prob->addVectorOperator(*opUGradU0, i, &minus1, &minus1);
+        self::prob->addVectorOperator(*opUGradU0, i, &minus1, &minus1);
     }
     
     /// < u*grad(u_i^old) , psi >
     if (nonLinTerm == 1 || nonLinTerm == 3) {
-      for (size_t k = 0; k < dow; ++k) {
-        Operator *opUGradU1 = new Operator(getFeSpace(i), getFeSpace(k));
+      for (size_t k = 0; k < self::dow; ++k) {
+        Operator *opUGradU1 = new Operator(feSpace_u, feSpace_u);
         opUGradU1->addTerm(new PartialDerivative_ZOT(
-	  prob->getSolution()->getDOFVector(i), k, density));
+	  self::prob->getSolution()->getDOFVector(i), k, density));
         prob->addMatrixOperator(*opUGradU1, i, k);
       }
     }
 
     /// < u^old*grad(u_i) , psi >
     if (nonLinTerm == 2 || nonLinTerm == 3) {
-      for (size_t k = 0; k < dow; ++k) {
-        Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
+      for (size_t k = 0; k < self::dow; ++k) {
+        Operator *opUGradU2 = new Operator(feSpace_u, feSpace_u);
         opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
 	  vel[k], k, density), GRD_PHI);
-        prob->addMatrixOperator(*opUGradU2, i, i);
+        self::prob->addMatrixOperator(*opUGradU2, i, i);
       }
     }
 
@@ -165,57 +168,59 @@ void NavierStokes_TaylorHood::fillOperators()
     addLaplaceTerm(i);
   
     /// < p , d_i(psi) >
-    Operator *opGradP = new Operator(getFeSpace(i), getFeSpace(dow));
+    Operator *opGradP = new Operator(feSpace_u, feSpace_p);
     opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
-    prob->addMatrixOperator(*opGradP, i, dow);
+    self::prob->addMatrixOperator(*opGradP, i, self::dow);
   
     /// external force, i.e. gravitational force
     if (std::abs(force[i]) > DBL_TOL) {
-      Operator *opForce = new Operator(getFeSpace(i), getFeSpace(i));
+      Operator *opForce = new Operator(feSpace_u, feSpace_u);
       opForce->addTerm(new Simple_ZOT(density*force[i]));
-      prob->addVectorOperator(*opForce, i);
+      self::prob->addVectorOperator(*opForce, i);
     }
   }
 
   /// div(u) = 0
   for (size_t i = 0; i < dow; ++i) {
     /// < d_i(u'_i) , psi >
-    Operator *opDivU = new Operator(getFeSpace(dow), getFeSpace(i));
+    Operator *opDivU = new Operator(feSpace_p, feSpace_u);
     opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
-    prob->addMatrixOperator(*opDivU, dow, i);
+    self::prob->addMatrixOperator(*opDivU, self::dow, i);
   }
 
-  Operator *opNull = new Operator(getFeSpace(dow), getFeSpace(dow));
+  Operator *opNull = new Operator(feSpace_p, feSpace_p);
   opNull->addTerm(new Simple_ZOT(0.0));
-  prob->addMatrixOperator(*opNull, dow, dow);
+  self::prob->addMatrixOperator(*opNull, self::dow, self::dow);
 }
 
 
 void NavierStokes_TaylorHood::addLaplaceTerm(int i)
 { FUNCNAME("NavierStokes_TaylorHood::addLaplaceTerm()");
 
+  const FiniteElemSpace* feSpace_u = self::getFeSpace(0);
+  
   /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
   if (laplaceType == 1) {
     for (size_t j = 0; j < dow; ++j) {
-      Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
+      Operator *opLaplaceUi1 = new Operator(feSpace_u, feSpace_u);
       opLaplaceUi1->addTerm(new MatrixIJ_SOT(j, i, viscosity));
-      prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta, &theta);
+      self::prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta, &theta);
 
       if (std::abs(minusTheta1) > DBL_TOL) {