[cig-commits] commit: compute_Cp uses FTensor

[Mercurial]

changeset:   87:d068549788a4
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Mon Mar 12 20:12:48 2012 -0700
files:       compute_Cp.cxx
description:
compute_Cp uses FTensor


diff -r a55f7cbcd509 -r d068549788a4 compute_Cp.cxx
--- a/compute_Cp.cxx	Mon Mar 12 18:07:24 2012 -0700
+++ b/compute_Cp.cxx	Mon Mar 12 20:12:48 2012 -0700
@@ -12,29 +12,43 @@ void compute_Cp(const Model &model, cons
                 double &Cp)
 {
   Cp=0;
-  if(model==Model::solCx)
-    {
-      const double x=(i+0.5)*h;
-      if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
-        {
-          FTensor::Tensor1<double,2> v, dv, dir(1,0);
-          compute_v_on_interface(zx,zy,log_etax,log_etay,
-                                 FTensor::Tensor1<double,2>(middle,0),
-                                 i,j,0,v,dv);
 
-          double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
+  FTensor::Tensor1<bool,2> intersects;
 
-          FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
-          const FTensor::Index<'a',2> a;
-          const FTensor::Index<'b',2> b;
-          FTensor::Number<0> n;
-          FTensor::Number<1> t;
+  FTensor::Tensor1<double,2> pos((i+0.5)*h,(j+0.5)*h);
+  intersects(0)=(pos(0)+h/2>middle && pos(0)-h/2<middle);
+  intersects(1)=false;
+  
+  FTensor::Tensor1<double,2> interface_pos[2];
+  interface_pos[0](0)=middle;
+  interface_pos[0](1)=pos(1);
+  interface_pos[1](0)=0;
+  interface_pos[1](1)=0;
 
-          FTensor::Tensor2<double,2,2> dz_jump;
-          dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
-          dz_jump(a,t)=eta_jump*dv(a);
+  for(int dd=0;dd<2;++dd)
+    if(intersects(dd))
+      {
+        FTensor::Tensor1<double,2> v, dv;
+        FTensor::Tensor1<int,2> dir(0,0);
+        dir(dd)=1;
+        compute_v_on_interface(zx,zy,log_etax,log_etay,
+                               interface_pos[dd],
+                               i,j,dd,v,dv);
 
-          Cp=(-eta_jump*v(a)*dir(a) + dx*dz_jump(a,b)*dir(a)*dir(b))/h;
-        }
-    }
+        double dx=(pos(dd)>interface_pos[dd](dd))
+          ? (interface_pos[dd](dd)-(pos(dd)-h/2))
+          : (interface_pos[dd](dd)-(pos(dd)+h/2));
+
+        FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
+        const FTensor::Index<'a',2> a;
+        const FTensor::Index<'b',2> b;
+        FTensor::Number<0> n;
+        FTensor::Number<1> t;
+
+        FTensor::Tensor2<double,2,2> dz_jump;
+        dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
+        dz_jump(a,t)=eta_jump*dv(a);
+
+        Cp=(-eta_jump*v(a)*dir(a) + dx*dz_jump(a,b)*dir(a)*dir(b))/h;
+      }
 }