Still not working

src/solutioninstance.h

@@ -21,7 +21,7 @@ namespace td{
     vd rho_,m_,rhoE_;
     d time=0;
   public:
-    SolutionInstance(int gp,d rho=1.2);
+    SolutionInstance(int gp,d rho=0);
     ~SolutionInstance(){}
     d getTime() const {return time;}
     void setTime(d t) {time=t;}

src/tube.cpp

@@ -9,7 +9,7 @@ d Tube::pleft(d t) {
   d pleft = std::sin(gc.getomg() * t);
   // cout << "t:" << t << endl;
   // cout << "omg:" << gc.getomg() << endl;
-  cout << "pleft:"<< pleft<<endl;
+  // cout << "pleft:"<< pleft<<endl;
   return pleft;
 }
   Tube::Tube(d L, int gp) throw(int) : L(L), gp(gp), sol(SolutionInstance(gp)) {
@@ -49,31 +49,31 @@ void TubeLF::Integrate(d dt) {
   vd Mflux=sol.Mflux();
   vd Eflux=sol.Eflux();
 
-  // std::cout << "poL" << sol.p()(0) << std::endl;
   {  // Left boundary
     d la = dt/dx;
     rho(0) = oldrho(0);
     rho(0)+=-la*(Cflux(1)-Cflux(0));
     d oldu0=oldm(0)/oldrho(0);
     // std::cout << "oldu:"<< oldu0 << std::endl;
-    d momfluxl = oldrho(0)*pow(oldu0, 2) + pleft(t);
+    d momfluxl = pow(oldm(0), 2)/oldrho(0) + pleft(t);
     m(0) = oldm(0);
     m(0)+=-la*(Mflux(1)-momfluxl);
     rhoE(0) = oldrhoE(0);
     rhoE(0)+=-la*(Eflux(1)-Eflux(0));
-
   } //  End left boundary
+
   {				// Inner nodes
     d lambda = dt/(2*dx);
-    rho.subvec(1,gp-2)=0.5*(oldrho.subvec(0,gp-3) + oldrho.subvec(2,gp-1));
-    rho.subvec(1,gp-2)+=-lambda*(Cflux.subvec(0,gp-3) -Cflux.subvec(2,gp-1));
+    rho.subvec(1,gp-2)=0.5*(oldrho.head(gp-2) + oldrho.tail(gp-2));
+    rho.subvec(1,gp-2)+=-lambda*(Cflux.head(gp-2) -Cflux.tail(gp-2));
   
-    m.subvec(1,gp-2)=0.5*(oldm.subvec(0,gp-3) + oldm.subvec(2,gp-1));
-    m.subvec(1,gp-2)+=-lambda*(Mflux.subvec(0,gp-3) -Mflux.subvec(2,gp-1));
+    m.subvec(1,gp-2)=0.5*(oldm.head(gp-2) + oldm.tail(gp-2));
+    m.subvec(1,gp-2)+=-lambda*(Mflux.head(gp-2) -Mflux.tail(gp-2));
 
-    rhoE.subvec(1,gp-2)=0.5*(oldrhoE.subvec(0,gp-3) + oldrhoE.subvec(2,gp-1));
-    rhoE.subvec(1,gp-2)+=-lambda*(Eflux.subvec(0,gp-3) -Eflux.subvec(2,gp-1));
+    rhoE.subvec(1,gp-2)=0.5*(oldrhoE.head(gp-2) + oldrhoE.tail(gp-2));
+    rhoE.subvec(1,gp-2)+=-lambda*(Eflux.head(gp-2) -Eflux.tail(gp-2));
   } // End inner nodes
+
   {				// Right boundary
     int i = gp - 1;
     d la = dt/dx;

src/tube.h

@@ -17,7 +17,7 @@ namespace td {
   typedef double d;
 
   class Tube {
-  protected:
+  public:
     d dx, L;              // Grid spacing, total length
     int gp;               // Number of gridpoints
     SolutionInstance sol; // Solutions at time instances

timedomain.i

@@ -43,6 +43,10 @@ class vd{
 %typemap(out) vd {
   $result=npy_from_vd($1);
 }
+%typemap(out) vd& {
+  const vd& res=*$1;
+  $result=npy_from_vd(res);
+}
 typedef double d;
 typedef unsigned us;
 
@@ -98,24 +102,34 @@ namespace td{
   public:
     SolutionInstance(int gp,d rho=1.2);
     ~SolutionInstance(){}
-    vd rho() const;
     vd p() const;
     vd u() const;
+    const vd& rho() const;
+    const vd& m() const;
+    const vd& rhoE() const;
+
     void setTime(d t);
     void setrho(d rho);
   };
 
   tasystem::Globalconf gc;
 
-  class Tube{
-    virtual void Integrate(d dt)=0;
+  class Tube {
   public:
-    Tube(double L,int gp) throw(int);
-    virtual ~Tube();
-    SolutionInstance& getSol();
-    void setSol(const SolutionInstance& sol);
-    void DoIntegration(d dt,int n=1);
-    d getTime();
+    d dx, L;              // Grid spacing, total length
+    int gp;               // Number of gridpoints
+    SolutionInstance sol; // Solutions at time instances
+    d t = 0;              // Current time
+    d pleft(d t);         // Compute pressure bc
+    virtual void Integrate(d dt) = 0;
+
+  public:
+    Tube(double L, int gp) throw(int);
+    virtual ~Tube() {}
+    SolutionInstance &getSol();
+    void setSol(const SolutionInstance &sol);
+    void DoIntegration(d dt, int n = 1);
+    d getTime() { return t; }
   };
   class TubeLF:public Tube{
     virtual void Integrate(d dt);