[cig-commits] commit: Fix interpolation on the boundary

[Mercurial]

changeset:   46:fdae064362c2
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Mon Feb 13 20:06:14 2012 -0800
files:       compute_v_on_interface.cxx
description:
Fix interpolation on the boundary


diff -r f2e7db30c967 -r fdae064362c2 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx	Mon Feb 13 15:39:07 2012 -0800
+++ b/compute_v_on_interface.cxx	Mon Feb 13 20:06:14 2012 -0800
@@ -13,6 +13,8 @@ namespace {
                      const double &dy,
                      double m[], double rhs[], const int sign=1)
   {
+    if(j<0 || j>Ny-1)
+      return;
     rhs[n]+=sign*z[i][j]*exp(-log_eta[i][j]);
     m[0+n*n_terms]+=sign*u0;
     m[1+n*n_terms]+=sign*dy;
@@ -86,6 +88,19 @@ void interpolate_to_surface(double z[][N
       add_to_matrix(n_terms,8,z,log_eta,i_center+di,j+1,1,-di*(dx-h),0,h,m,rhs);
       add_to_matrix(n_terms,8,z,log_eta,i_center+di,j-1,1,-di*(dx-h),0,
                     -h,m,rhs,-1);
+
+      /* At the j==0 boundary, all of the j-1 terms are not applied.
+         Instead, for n=2, apply dvx/dy=0. */
+      if(j==0)
+        {
+          m[1+2*n_terms]=1;
+          m[2+2*n_terms]=-h/2;
+        }
+      else if(j==Ny-1)
+        {
+          m[1+1*n_terms]=1;
+          m[2+1*n_terms]=h/2;
+        }
     }
   else
     {
@@ -124,7 +139,7 @@ void interpolate_to_surface(double z[][N
                     h/2,m,rhs);
     }
 
-  // if(j==128)
+  // if(j>510)
   //   for(int a=0;a<9;++a)
   //     {
   //       std::cout << a << " "
@@ -169,10 +184,21 @@ void compute_v_on_interface(double zx[Nx
       jy=j;
     }
 
+  // std::cout << "compute v "
+  //           << x << " "
+  //           << j << " "
+  //           << xyz << " "
+  //           << "\n";
+
   /* vx */
   int ix(x/h);
   double dx_x(x-ix*h);
   interpolate_to_surface(zx,log_etax,ix,jx,dx_x,xyz==0,vx,dvx_y,dvx_yy);
+
+  // std::cout 
+  //           << vx << " "
+  //           << "\n";
+
 
   /* vy */
   int iy=x/h-0.5;