[cig-commits] commit: Compute corner cutting correction correctly.

[Mercurial]

changeset:   147:551aef7421a2
user:        Walter Landry <wlandry at caltech.edu>
date:        Mon Apr 02 20:06:10 2012 -0700
files:       compute_coefficients/compute_Cxyz.cxx
description:
Compute corner cutting correction correctly.


diff -r 0b2d53a5d42f -r 551aef7421a2 compute_coefficients/compute_Cxyz.cxx
--- a/compute_coefficients/compute_Cxyz.cxx	Mon Apr 02 20:04:00 2012 -0700
+++ b/compute_coefficients/compute_Cxyz.cxx	Mon Apr 02 20:06:10 2012 -0700
@@ -38,6 +38,7 @@ void compute_Cxyz(const double zx[Nx+1][
   const Interface interface(d,i,j,distx,disty);
   if(!interface.intersects())
     return;
+
   FTensor::Tensor1<double,2> pos;
   const FTensor::Index<'a',2> a;
   const FTensor::Index<'b',2> b;
@@ -80,7 +81,18 @@ void compute_Cxyz(const double zx[Nx+1][
                   eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
                   break;
                 case 1:
-                  eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
+                  if(j==0)
+                    {
+                      eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j]);
+                    }
+                  else if(j==Ny-1)
+                    {
+                      eta_jump*=exp(log_etax[i][j])-exp(log_etax[i][j-1]);
+                    }
+                  else
+                    {
+                      eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
+                    }
                   break;
                 default:
                   abort();
@@ -91,7 +103,18 @@ void compute_Cxyz(const double zx[Nx+1][
               switch(dd)
                 {
                 case 0:
-                  eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
+                  if(i==0)
+                    {
+                      eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i][j]);
+                    }
+                  else if(i==Nx-1)
+                    {
+                      eta_jump*=exp(log_etay[i][j])-exp(log_etay[i-1][j]);
+                    }
+                  else
+                    {
+                      eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
+                    }
                   break;
                 case 1:
                   eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
@@ -180,14 +203,16 @@ void compute_Cxyz(const double zx[Nx+1][
                 }
               compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
 
-              FTensor::Tensor1<double,2> dx;
-              for(int aa=0;aa<2;++aa)
-                dx(aa)=pos(aa) - interface.corner_pos[dd][ee](aa)
-                  - (h/2)*sign(pos(aa) - interface.corner_pos[dd][ee](aa));
+              FTensor::Tensor1<double,2> dx,yx_pos;
 
-              C+=sign(dx(0))*sign(dx(1))*zyx(a)
-                *(z_jump(a) + dz_jump(a,b)*dx(b)
-                  + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+              yx_pos(0)=pos(0)+(dd-0.5)*h;
+              yx_pos(1)=pos(1)+(ee-0.5)*h;
+
+              dx(a)=yx_pos(a)-interface.corner_pos[dd][ee](a);
+
+              C+=sign(dd-0.5)*sign(ee-0.5)*
+                zyx(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
+                        + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
             }
     }
 }