BaseProblems updated and some Warnings removed

extensions/base_problems/NavierStokesPhase_TaylorHood.cc

@@ -147,20 +147,14 @@ void NavierStokesPhase_TaylorHood::addLaplaceTerm(int i)
       for (unsigned j = 0; j < dow; ++j) {
         Operator *opLaplaceUi1 = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
         opLaplaceUi1->addTerm(new MatrixIJPhase_SOT(phase, j, i, viscosity));
-        prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta);
-
-        opLaplaceUi1->setUhOld(prob->getSolution()->getDOFVector(j));
-        prob->addVectorOperator(*opLaplaceUi1, i, &minusTheta1);
+        prob->addMatrixOperator(*opLaplaceUi1, i, j);
       }
     }
     
     /// < alpha*grad(u'_i) , grad(psi) >
     Operator *opLaplaceUi = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
       opLaplaceUi->addTerm(new Phase_SOT(phase, viscosity));
-    prob->addMatrixOperator(*opLaplaceUi, i, i, &theta);
-
-    opLaplaceUi->setUhOld(prob->getSolution()->getDOFVector(i));
-    prob->addVectorOperator(*opLaplaceUi, i, &minusTheta1);
+    prob->addMatrixOperator(*opLaplaceUi, i, i);
 }
 
 

extensions/base_problems/NavierStokesPhase_TaylorHood.h

@@ -84,7 +84,7 @@ public: // methods
   virtual void fillBoundaryConditions();
   virtual void addLaplaceTerm(int i);
 
-private: // variables
+protected: // variables
   
   double dbc_factor; // beta/epsilon^alpha
   double beta;

extensions/base_problems/NavierStokes_TH_MultiPhase.cc

@@ -18,8 +18,10 @@ NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_)
   Initfile::get(name + "->density1", density1); // density of fluid 1
   Initfile::get(name + "->density2", density2); // density of fluid 2
 
-  if (viscosity1 <= 0.0)
+  if (viscosity1 <= 0.0 && viscosity > 0.0)
     viscosity1 = viscosity;
+  
+  TEST_EXIT(viscosity1 > 0.0)("Positive viscosity is necessary!\n");
 };
 
 
@@ -71,14 +73,14 @@ void NavierStokes_TH_MultiPhase::fillOperators()
     /// < (1/tau)*rho*u'_i , psi >
     Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
     if (density1==density2)
-      opTime->addTerm(new Simple_ZOT());  
+      opTime->addTerm(new Simple_ZOT(density1));  
     else
       opTime->addTerm(new VecAtQP_ZOT(densityPhase, NULL));
     prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
     /// < (1/tau)*rho*u_i^old , psi >
     Operator *opTimeOld = new Operator(getFeSpace(i), getFeSpace(i));
     if (density1==density2)
-      opTimeOld->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(i), NULL));  
+      opTimeOld->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(i), new AMDiS::Factor(density1)));  
     else
       opTimeOld->addTerm(new Vec2AtQP_ZOT(densityPhase, prob->getSolution()->getDOFVector(i), NULL));
     prob->addVectorOperator(*opTimeOld, i, getInvTau(), getInvTau());
@@ -86,7 +88,7 @@ void NavierStokes_TH_MultiPhase::fillOperators()
     /// < u^old*grad(u_i^old) , psi >
     Operator *opUGradU0 = new Operator(getFeSpace(i),getFeSpace(i));
     if (density1==density2)
-        opUGradU0->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PHI);
+        opUGradU0->addTerm(new WorldVec_FOT(vel, -density1), GRD_PHI);
     else
         opUGradU0->addTerm(new WorldVecPhase_FOT(densityPhase, vel, -1.0), GRD_PHI);
     opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
@@ -145,7 +147,7 @@ void NavierStokes_TH_MultiPhase::addLaplaceTerm(int i)
 { FUNCNAME("NavierStokes_TH_MultiPhase::addLaplaceTerm()");
 
   /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
-  if (viscosity1!=viscosity2) {
+  if (viscosity1 != viscosity2) {
     for (unsigned j = 0; j < dow; ++j) {
       Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
       opLaplaceUi1->addTerm(new VecAtQP_IJ_SOT(  viscosityPhase, NULL, j, i));
@@ -157,8 +159,13 @@ void NavierStokes_TH_MultiPhase::addLaplaceTerm(int i)
   
   /// < alpha*grad(u'_i) , grad(psi) >
   Operator *opLaplaceUi = new Operator(getFeSpace(i), getFeSpace(i));
-  opLaplaceUi->addTerm(new VecAtQP_SOT(viscosityPhase, NULL));
+  if (viscosity1 != viscosity2) {
+    opLaplaceUi->addTerm(new VecAtQP_SOT(viscosityPhase, NULL));
+  } else {
+    opLaplaceUi->addTerm(new Simple_SOT(viscosity1));
+  }
   opLaplaceUi->setUhOld(getSolution()->getDOFVector(i));
   prob->addMatrixOperator(*opLaplaceUi, i, i, &theta, &theta);  
   prob->addVectorOperator(*opLaplaceUi, i, &minusTheta1, &minusTheta1);
+  
 };