[cig-commits] commit: Change dCx to use FTensor

[Mercurial]

changeset:   79:168fc9952868
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Mar 11 01:47:57 2012 -0800
files:       compute_Cx.cxx
description:
Change dCx to use FTensor


diff -r e14e5fab7ba8 -r 168fc9952868 compute_Cx.cxx
--- a/compute_Cx.cxx	Sat Mar 10 20:11:14 2012 -0800
+++ b/compute_Cx.cxx	Sun Mar 11 01:47:57 2012 -0800
@@ -52,12 +52,10 @@ void compute_Cx(const Model &model, cons
 
           double p_jump=-2*eta_jump*dv(1);
 
-          FTensor::Tensor1<double,2> dp_jump;
-          dp_jump(t)=-2*eta_jump*ddv(t);
-          dp_jump(n)=2*eta_jump*ddv(n);
+          FTensor::Tensor1<double,2> dp_jump(2*eta_jump*ddv(n),
+                                             -2*eta_jump*ddv(t));
 
           FTensor::Tensor3_christof<double,2,2> ddz_jump;
-
           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);
@@ -71,21 +69,11 @@ void compute_Cx(const Model &model, cons
           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);
@@ -94,27 +82,42 @@ void compute_Cx(const Model &model, cons
           double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
           double dvx_dv=(eta*(4*dvp_dv - dvpp_dv)  + 2*h*dzx_x_dv)
             /(3*(min_eta+max_eta));
-          double dp_x_dv=2*eta_jump*dvx_yy_dv;
-          double dzx_xx_dv=-eta_jump*dvx_yy_dv + dp_x_dv;
-          double dp_dv=-2*eta_jump*dvy_y_dv;
 
-          dCx_dvx=-dz_dd_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
-                                  + dx*dx*dzx_xx_dv/2)/(h*h);
+
+          FTensor::Tensor1<double,2> ddv_dv(dvx_yy_dv,0);
+          FTensor::Tensor1<double,2> dv_dv(0,dvy_y_dv);
+          FTensor::Tensor1<double,2> v_dv(dvx_dv,0);
+
+          FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
+                                           -2*eta_jump*ddv_dv(t));
+          double p_dv=-2*eta_jump*dv_dv(t);
+
+          FTensor::Tensor2<double,2,2> dz_dv;
+          dz_dv(a,n)=-eta_jump*dv_dv(b)*ddel(a,b);
+          dz_dv(a,t)=eta_jump*dv_dv(a);
+
+          FTensor::Tensor3_christof<double,2,2> ddz_dv;
+          ddz_dv(a,t,t)=eta_jump*ddv_dv(a);
+          ddz_dv(a,n,n)=-ddz_dv(a,t,t) + dp_dv(a);
+          ddz_dv(a,n,t)=-eta_jump*ddv_dv(b)*ddel(a,b);
+
+          dCx_dvx=-dz_dd_sign*2*(eta_jump*v_dv(a)*xyz(a)
+                                 - dx*dz_dv(a,b)*xyz(a)*dir(b)
+                                 + dx*dx*ddz_dv(a,b,c)*xyz(a)*dir(b)*dir(c)/2)
+            /(h*h);
 
           if(x+h/2>middle && x-h/2<middle)
             {
               dx=(x>middle) ? ((x-h/2)-middle)
                 : ((x+h/2)-middle);
 
-              double dp_d_jump=dp_jump(a)*dir(a);
-              Cx+=(p_jump + dx*dp_d_jump)/h;
+              Cx+=(p_jump + dx*dp_jump(a)*dir(a))/h
+                - (dz_jump(a,b)*dir2(a)*dir3(b)
+                   + dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
 
-              double dz2_3_jump=dz_jump(a,b)*dir2(a)*dir3(b);
-              double ddz2_34_jump=ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir4(c);
-              Cx-=(dz2_3_jump + dx*ddz2_34_jump)/h;
-
-              dCx_dvx+= -(eta_jump*dvy_y_dv + dx*dzy_yx_dv)/h
-                + (dp_dv + dx*dp_x_dv)/h;
+              dCx_dvx+=(p_dv + dx*dp_dv(a)*dir(a))/h
+                - (dz_dv(a,b)*dir2(a)*dir3(b)
+                   + dx*ddz_dv(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
             }
           dCx_dzx=dCx_dvx/eta;
         }