[cig-commits] commit: Convert the jump calculations for dp/dx and ddzy/dxdy for Cx to FTensor

[Mercurial]

changeset:   78:e14e5fab7ba8
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Sat Mar 10 20:11:14 2012 -0800
files:       compute_Cx.cxx
description:
Convert the jump calculations for dp/dx and ddzy/dxdy for Cx to FTensor


diff -r e7346f4e5961 -r e14e5fab7ba8 compute_Cx.cxx
--- a/compute_Cx.cxx	Sat Mar 10 15:11:33 2012 -0800
+++ b/compute_Cx.cxx	Sat Mar 10 20:11:14 2012 -0800
@@ -27,8 +27,12 @@ void compute_Cx(const Model &model, cons
                                  vx,vy,dvx_y,dvy_y,
                                  dvx_yy,dvy_yy);
 
+          /* xyz= the component we are correcting, Cx, Cy, or Cz.
+             dir=direction of derivative
+             dir2, dir3=components of mixed derivative we are correcting.
+             So for zy,xy with boundary in x, dir2=y, dir3=y. */
           FTensor::Tensor1<double,2> v(vx,0), dv(0,dvy_y), ddv(dvx_yy,0),
-            xyz(1,0), dir(1,0);
+            xyz(1,0), dir(1,0), dir2(0,1), dir3(0,1), dir4(1,0);
           FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
           const FTensor::Index<'a',2> a;
           const FTensor::Index<'b',2> b;
@@ -49,9 +53,8 @@ void compute_Cx(const Model &model, cons
           double p_jump=-2*eta_jump*dv(1);
 
           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);
+          dp_jump(t)=-2*eta_jump*ddv(t);
+          dp_jump(n)=2*eta_jump*ddv(n);
 
           FTensor::Tensor3_christof<double,2,2> ddz_jump;
 
@@ -66,6 +69,13 @@ void compute_Cx(const Model &model, cons
               + 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;
@@ -96,9 +106,12 @@ void compute_Cx(const Model &model, cons
               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-=eta_jump*(dvy_y - dx*dvx_yy)/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;