[cig-commits] commit: switch jump computation out to compute_jumps

[Mercurial]

changeset:   115:06aaacbd9afb
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Mar 25 19:14:56 2012 -0700
files:       compute_Cxyz.cxx compute_jumps.cxx wscript
description:
switch jump computation out to compute_jumps


diff -r 2690d006b243 -r 06aaacbd9afb compute_Cxyz.cxx
--- a/compute_Cxyz.cxx	Sun Mar 25 16:28:14 2012 -0700
+++ b/compute_Cxyz.cxx	Sun Mar 25 19:14:56 2012 -0700
@@ -3,6 +3,15 @@
 #include "compute_v_on_interface.hxx"
 #include "FTensor.hpp"
 #include "Interface.hxx"
+
+void compute_jumps(const FTensor::Tensor1<double,2> &v,
+                   const FTensor::Tensor1<double,2> &dv,
+                   const FTensor::Tensor1<double,2> &ddv, 
+                   FTensor::Tensor1<double,2> &z_jump,
+                   FTensor::Tensor2<double,2,2> &dz_jump,
+                   FTensor::Tensor3_christof<double,2,2> &ddz_jump,
+                   double &p_jump,
+                   FTensor::Tensor1<double,2> &dp_jump);
 
 void compute_Cxyz(const double zx[Nx+1][Ny],
                   const double zy[Nx][Ny+1],
@@ -41,38 +50,22 @@ void compute_Cxyz(const double zx[Nx+1][
                                  interface.pos[dd],d,
                                  v,dv,ddv);
 
-          FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
-  
-          double dx=(pos(dd)>interface.pos[dd](dd))
-            ? (interface.pos[dd](dd)-(pos(dd)-h))
-            : (interface.pos[dd](dd)-(pos(dd)+h));
-          double z_jump=eta_jump*v(a)*xyz(a);
+          
+          FTensor::Tensor1<double,2> dx(0,0);
+          dx(dd)=interface.pos[dd](dd) - pos(dd)
+            - h*sign(interface.pos[dd](dd)-pos(dd));
+          FTensor::Tensor1<double,2> z_jump;
+          FTensor::Tensor2<double,2,2> dz_jump;
+          FTensor::Tensor3_christof<double,2,2> ddz_jump;
+          double p_jump;
+          FTensor::Tensor1<double,2> dp_jump;
 
-          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);
-
-          /* TODO: dp_jump and ddz_jump both need to be corrected for
-             curvature and density jumps. */
-
-          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);
-          
-          double ddz_dd_jump=ddz_jump(a,b,c)*xyz(a)*dir(b)*dir(c);
+          compute_jumps(v,dv,ddv,z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
 
           double dz_dd_correction=
-            -(z_jump - dx*dz_d_jump + dx*dx*ddz_dd_jump/2)/(h*h);
+            -xyz(a)*(z_jump(a) - dz_jump(a,b)*dx(b)
+                     + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+
           const int dz_dd_factor((pos(dd)>interface.pos[dd](dd) ? -1 : 1)
                                  *(d==dd ? 2 : 1));
 
@@ -80,11 +73,11 @@ void compute_Cxyz(const double zx[Nx+1][
 
           if(d==dd && interface.intersect_dp)
             {
-              dx=(pos(dd)>interface.pos[dd](dd))
+              dx(dd)=(pos(dd)>interface.pos[dd](dd))
                 ? ((pos(dd)-h/2)-interface.pos[dd](dd))
                 : ((pos(dd)+h/2)-interface.pos[dd](dd));
 
-              C+=(p_jump + dx*dp_jump(a)*dir(a))/h;
+              C+=(p_jump + dp_jump(a)*dx(a))/h;
             }
 
           if(interface.intersect_sides[dd])
@@ -93,12 +86,12 @@ void compute_Cxyz(const double zx[Nx+1][
               FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
               dir2((d+1)%2)=1;
               dir3((dd+1)%2)=1;
-              dx=(pos(dd)>interface.pos[dd](dd))
+              dx(dd)=(pos(dd)>interface.pos[dd](dd))
                 ? ((pos(dd)-h/2)-interface.pos[dd](dd))
                 : ((pos(dd)+h/2)-interface.pos[dd](dd));
 
               C+=-(dz_jump(a,b)*dir2(a)*dir3(b)
-                   + dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
+                   + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
             }
         }
     }
diff -r 2690d006b243 -r 06aaacbd9afb compute_jumps.cxx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_jumps.cxx	Sun Mar 25 19:14:56 2012 -0700
@@ -0,0 +1,36 @@
+#include "FTensor.hpp"
+#include "constants.hxx"
+
+void compute_jumps(const FTensor::Tensor1<double,2> &v,
+                   const FTensor::Tensor1<double,2> &dv,
+                   const FTensor::Tensor1<double,2> &ddv, 
+                   FTensor::Tensor1<double,2> &z_jump,
+                   FTensor::Tensor2<double,2,2> &dz_jump,
+                   FTensor::Tensor3_christof<double,2,2> &ddz_jump,
+                   double &p_jump,
+                   FTensor::Tensor1<double,2> &dp_jump)
+{
+  const FTensor::Index<'a',2> a;
+  const FTensor::Index<'b',2> b;
+
+  z_jump(a)=eta_jump*v(a);
+  
+  FTensor::Number<0> n;
+  FTensor::Number<1> t;
+  FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+  
+  dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
+  dz_jump(a,t)=eta_jump*dv(a);
+
+  p_jump=-2*eta_jump*dv(1);
+
+  /* TODO: dp_jump and ddz_jump both need to be corrected for
+     curvature and density jumps. */
+
+  dp_jump(n)=2*eta_jump*ddv(n);
+  dp_jump(t)=-2*eta_jump*ddv(t);
+
+  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);
+}
diff -r 2690d006b243 -r 06aaacbd9afb wscript
--- a/wscript	Sun Mar 25 16:28:14 2012 -0700
+++ b/wscript	Sun Mar 25 19:14:56 2012 -0700
@@ -13,6 +13,7 @@ def build(bld):
                         'compute_Cxyz.cxx',
                         'compute_Cp.cxx',
                         'compute_dC_dv.cxx',
+                        'compute_jumps.cxx',
                         'Interface/Interface.cxx'],
 
         target       = 'Gamr_disc',