diff --git a/src/navier_stokes.cc b/src/navier_stokes.cc
index b291234d304fd19815190acdede9d9fbc2dd4f0a..0fef4a308c069a52e3b4cf5b2e7176fa56f1dbef 100644
--- a/src/navier_stokes.cc
+++ b/src/navier_stokes.cc
@@ -50,17 +50,17 @@ int main(int argc, char** argv)
     For<0,DOW>::loop([&](auto const _i) 
     {    
 	// <1/tau * u_i, v_i>
-	auto opTime = Op::create(), opTimeOld = Op::create();
+	Op *opTime = new Op, *opTimeOld = new Op;
 	opTime->addZOT( density );
         opTimeOld->addZOT( valueOf(prob.getSolution(_i), density) );
-        prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
-        prob.addVectorOperator(opTimeOld,  _i, probInstat.getInvTau());
+        prob.addMatrixOperator(*opTime, _i, _i, probInstat.getInvTau());
+        prob.addVectorOperator(*opTimeOld,  _i, probInstat.getInvTau());
         
 #if 0
         // <(u * nabla)u_i^old, v_i>
         For<0, DOW>::loop([&](auto const _j)
         {
-          auto opNonlin = Op::create();
+          auto opNonlin = new Op;
           opNonlin->addZOT( derivativeOf(prob.getSolution(_i), _j, density) );
           prob.addMatrixOperator(opNonlin, _i, _j);
         });
@@ -68,25 +68,25 @@ int main(int argc, char** argv)
         // <(u^old * nabla)u_i, v_i>
         For<0, DOW>::loop([&](auto const _j)
         {
-          auto opNonlin = Op::create();
+          Op *opNonlin = new Op;
           opNonlin->addFOT( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
-          prob.addMatrixOperator(opNonlin, _i, _i);
+          prob.addMatrixOperator(*opNonlin, _i, _i);
         });
 #endif      
         // <viscosity*grad(u_i), grad(v_i)>
-        auto opL = Op::create();
+        Op *opL = new Op;
 	opL->addSOT( viscosity ); 
-	prob.addMatrixOperator(opL, _i, _i);
+	prob.addMatrixOperator(*opL, _i, _i);
 	
 	// <p, d_i(v_i)>
-	auto opP = Op::create();
+	Op *opP = new Op;
 	opP->addFOT( 1.0, _i, GRD_PSI ); 
-	prob.addMatrixOperator(opP, _i, DOW); 
+	prob.addMatrixOperator(*opP, _i, DOW); 
     
 	// <d_i(u_i), q>
-	auto opDiv = Op::create(); 
+	Op *opDiv = new Op; 
 	opDiv->addFOT( 1.0, _i, GRD_PHI );
-	prob.addMatrixOperator(opDiv, DOW, _i);
+	prob.addMatrixOperator(*opDiv, DOW, _i);
     });
     
     // define boundary regions