[cig-commits] commit: Compute dvy_dtt using general method. This also fixes a bug where the

[Mercurial]

changeset:   99:cf2d88c6eb97
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Mar 20 13:21:02 2012 -0700
files:       compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx
description:
Compute dvy_dtt using general method.  This also fixes a bug where the
valid vector was not being reset for each direction.


diff -r 1794f9c80966 -r cf2d88c6eb97 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx	Tue Mar 20 13:14:25 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx	Tue Mar 20 13:21:02 2012 -0700
@@ -18,16 +18,11 @@ void compute_dv_dtt(const double zx[Nx+1
                     const FTensor::Tensor1<double,2> &pos,
                     const FTensor::Tensor1<double,2> &norm,
                     const FTensor::Tensor1<double,2> &tangent,
+                    const int &ii, const int &jj,
                     FTensor::Tensor1<double,2> &ddv,
                     FTensor::Tensor1<double,2> &ddv_dv)
 {
-  /* vx */
-
   const int max_r=4;
-  FTensor::Tensor2_symmetric<std::vector<Valid>,2> valid;
-  valid(0,0).resize((2*max_r+1)*(2*max_r+1));
-  valid(0,1).resize((2*max_r+1)*(2*max_r+1));
-  valid(1,1).resize((2*max_r+1)*(2*max_r+1));
 
   const FTensor::Index<'a',2> a;
   const FTensor::Index<'b',2> b;
@@ -40,6 +35,11 @@ void compute_dv_dtt(const double zx[Nx+1
 
   for(int d=0;d<2;++d)
     {
+      FTensor::Tensor2_symmetric<std::vector<Valid>,2> valid;
+      valid(0,0).resize((2*max_r+1)*(2*max_r+1));
+      valid(0,1).resize((2*max_r+1)*(2*max_r+1));
+      valid(1,1).resize((2*max_r+1)*(2*max_r+1));
+
       int max_x(Nx-d), max_y(Ny+d-1);
       double dx(d*h/2), dy((1-d)*h/2);
 
@@ -65,12 +65,12 @@ void compute_dv_dtt(const double zx[Nx+1
               {
                 bool v;
                 if(i==0)
-                  v=sign(disty[i][j])==sign(disty[i+1][j]);
+                  v=(sign(disty[i][j])==sign(disty[i+1][j]));
                 else if(i==max_x)
-                  v=sign(disty[i][j])==sign(disty[i-1][j]);
+                  v=(sign(disty[i][j])==sign(disty[i-1][j]));
                 else
-                  v=sign(disty[i][j])==sign(disty[i+1][j])
-                    && sign(disty[i][j])==sign(disty[i-1][j]);
+                  v=(sign(disty[i][j])==sign(disty[i+1][j])
+                     && sign(disty[i][j])==sign(disty[i-1][j]));
                 
                 valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
                   Valid(sign(disty[i][j]),v,i,j,diff(a)*diff(a));
@@ -83,6 +83,7 @@ void compute_dv_dtt(const double zx[Nx+1
           {
             FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
             diff(a)=p(a)-pos(a);
+
             if(d==0)
               {
                 bool v;
@@ -165,7 +166,6 @@ void compute_dv_dtt(const double zx[Nx+1
                         else
                           compute_2nd_derivs(pm,1,v,d0,d1,Nx,zy,log_etay,
                                              ddv_pm,eta_pm);
-
                       }
                   }
               }
diff -r 1794f9c80966 -r cf2d88c6eb97 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx	Tue Mar 20 13:14:25 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx	Tue Mar 20 13:21:02 2012 -0700
@@ -124,6 +124,7 @@ void compute_dv_dtt(const double zx[Nx+1
                     const FTensor::Tensor1<double,2> &pos,
                     const FTensor::Tensor1<double,2> &norm,
                     const FTensor::Tensor1<double,2> &tangent,
+                    const int &i, const int &j,
                     FTensor::Tensor1<double,2> &ddv,
                     FTensor::Tensor1<double,2> &ddv_dv);
 
@@ -168,6 +169,8 @@ void compute_v_on_interface(const double
 
   if(xyz==0)
     {
+      /* Need to compute norm correctly when the interface is
+         horizontal and vertical */
       FTensor::Tensor1<double,2> norm, tangent;
       norm(0)=((disty[iy+1][jy] - disty[iy][jy])*(h-dx_j)
                + (disty[iy+2][jy] - disty[iy+1][jy])*dx_j)/(h*h);
@@ -183,7 +186,7 @@ void compute_v_on_interface(const double
       
       FTensor::Tensor1<double,2> ddv_dv_temp;
       compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
-                     norm,tangent,ddv,ddv_dv_temp);
+                     norm,tangent,i,j,ddv,ddv_dv_temp);
 
       // FTensor::Tensor1<double,2> pos_p, pos_pp, pos_m, pos_mm;
       // pos_p(a)=pos(a)+length*norm(a);
@@ -223,7 +226,23 @@ void compute_v_on_interface(const double
     }
   else
     {
-      ddv(1)=compute_dv_yy(zy,log_etay,dx_j,iy,j);
+      /* Need to compute norm correctly when the interface is
+         horizontal and vertical */
+      FTensor::Tensor1<double,2> norm, tangent;
+      norm(0)=((distx[ix+1][jx] - distx[ix][jx])*(h-dx_i)
+               + (distx[ix+2][jx] - distx[ix+1][jx])*dx_i)/(h*h);
+      norm(1)=((distx[ix][jx+1] - distx[ix][jx])*(h-dx_i)
+               + (distx[ix+1][jx+1] - distx[ix+1][jx])*dx_i)/(h*h);
+
+      norm(a)/=sqrt(norm(b)*norm(b));
+
+      tangent(0)=-norm(1);
+      tangent(1)=norm(0);
+
+      FTensor::Tensor1<double,2> ddv_temp, ddv_dv_temp;
+      compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
+                     norm,tangent,i,j,ddv,ddv_dv_temp);
+
       dv(0)=compute_dv_y(zx,log_etax,dzx_yx_jump(ddv(1)),dx_i,ix,j-1);
       v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dx_j,iy,j);