[cig-commits] commit: Convert all of the jump calculation in compute_Cx to FTensor

[Mercurial]

changeset:   77:e7346f4e5961
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Sat Mar 10 15:11:33 2012 -0800
files:       compute_Cx.cxx
description:
Convert all of the jump calculation in compute_Cx to FTensor


diff -r 4db4b87a3fdf -r e7346f4e5961 compute_Cx.cxx
--- a/compute_Cx.cxx	Sat Mar 10 15:11:15 2012 -0800
+++ b/compute_Cx.cxx	Sat Mar 10 15:11:33 2012 -0800
@@ -27,40 +27,56 @@ void compute_Cx(const Model &model, cons
                                  vx,vy,dvx_y,dvy_y,
                                  dvx_yy,dvy_yy);
 
-          const int d(0);
-          FTensor::Tensor1<double,2> v(vx,0), dv(0,dvy_y), ddv(dvx_yy,0), norm(1,0),
-            tangent(0,1);
-          FTensor::Tensor2_antisymmetric<double,2> eps(1);
+          FTensor::Tensor1<double,2> v(vx,0), dv(0,dvy_y), ddv(dvx_yy,0),
+            xyz(1,0), dir(1,0);
+          FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
           const FTensor::Index<'a',2> a;
+          const FTensor::Index<'b',2> b;
+          const FTensor::Index<'c',2> c;
 
           double dx=(x>middle) ? (middle-(x-h))
             : (middle-(x+h));
-          double z_jump=eta_jump*v(d);
+          double z_jump=eta_jump*v(a)*xyz(a);
 
-          double dz_n_jump=eta_jump*(dv(0)*norm(a) - dv(1)*tangent(a))*eps(d,a);
-          double dz_t_jump=eta_jump*(dv(0)*norm(d) + dv(1)*tangent(d));
-          double dz_x_jump=dz_n_jump*norm(d) + dz_t_jump*tangent(d);
+          FTensor::Tensor2<double,2,2> dz_jump;
+          FTensor::Number<0> n;
+          FTensor::Number<1> t;
+
+          dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
+          dz_jump(a,t)=eta_jump*dv(a);
+          double dz_d_jump=dz_jump(a,b)*xyz(a)*dir(b);
 
           double p_jump=-2*eta_jump*dv(1);
 
-          double dp_t_jump=-2*eta_jump*ddv(1);
-          double dp_n_jump=2*eta_jump*ddv(0);
-          double dp_x_jump=dp_n_jump*norm(d) + dp_t_jump*tangent(d);
+          FTensor::Tensor1<double,2> dp_jump;
+          dp_jump(1)=-2*eta_jump*ddv(1);
+          dp_jump(0)=2*eta_jump*ddv(0);
+          double dp_d_jump=dp_jump(a)*dir(a);
 
-          double dz_tt=eta_jump*(ddv(0)*norm(d) + ddv(1)*tangent(d));
-          double dz_xx_jump=-dz_tt + dp_x_jump;
+          FTensor::Tensor3_christof<double,2,2> ddz_jump;
 
-          double dz_xx_correction=
-            -(z_jump - dx*dz_x_jump
-              + dx*dx*dz_xx_jump/2)/(h*h);
-          const int dz_xx_sign(x>middle ? -1 : 1);
-          Cx=dz_xx_sign*2*dz_xx_correction;
+          ddz_jump(a,t,t)=eta_jump*ddv(a);
+          ddz_jump(a,n,n)=-ddz_jump(a,t,t) + dp_jump(a);
+          ddz_jump(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
+
+          double ddz_dd_jump=ddz_jump(a,b,c)*xyz(a)*dir(b)*dir(c);
+
+          double dz_dd_correction=
+            -(z_jump - dx*dz_d_jump
+              + dx*dx*ddz_dd_jump/2)/(h*h);
+          const int dz_dd_sign(x>middle ? -1 : 1);
+          Cx=dz_dd_sign*2*dz_dd_correction;
 
           double eta=exp(log_etax[i][j]);
           double delta=(h-std::fabs(dx))/h;
           double dvx_yy_dv=-2*eta*eta/(h*h*(min_eta*min_eta + max_eta*max_eta));
           if(j==0 || j==Ny-1)
             dvx_yy_dv/=2;
+
+
+          // double ddv_dv=dvx_yy_dv,0;
+          // double dz_nt_dv=-eta_jump*ddv_dv;
+
           double dzy_yx_dv=-eta_jump*dvx_yy_dv;
           double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
           double dzx_x_dv=eta_jump*dvy_y_dv;
@@ -72,7 +88,7 @@ void compute_Cx(const Model &model, cons
           double dzx_xx_dv=-eta_jump*dvx_yy_dv + dp_x_dv;
           double dp_dv=-2*eta_jump*dvy_y_dv;
 
-          dCx_dvx=-dz_xx_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
+          dCx_dvx=-dz_dd_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
                                   + dx*dx*dzx_xx_dv/2)/(h*h);
 
           if(x+h/2>middle && x-h/2<middle)
@@ -80,7 +96,7 @@ void compute_Cx(const Model &model, cons
               dx=(x>middle) ? ((x-h/2)-middle)
                 : ((x+h/2)-middle);
 
-              Cx+=(p_jump + dx*dp_x_jump)/h;
+              Cx+=(p_jump + dx*dp_d_jump)/h;
               Cx-=eta_jump*(dvy_y - dx*dvx_yy)/h;