[cig-commits] commit: Remove most of the second derivative stuff.

Mercurial hg at geodynamics.org
Tue May 8 05:50:54 PDT 2012


changeset:   155:980e18575427
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun May 06 18:49:35 2012 -0700
files:       compute_coefficients/compute_Cp.cxx compute_coefficients/compute_Cxyz.cxx compute_coefficients/compute_jumps.cxx compute_coefficients/compute_v_on_interface.hxx compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx
description:
Remove most of the second derivative stuff.


diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_Cp.cxx
--- a/compute_coefficients/compute_Cp.cxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_Cp.cxx	Sun May 06 18:49:35 2012 -0700
@@ -8,7 +8,6 @@ void compute_jumps(const double &eta_jum
 void compute_jumps(const double &eta_jump,
                    const FTensor::Tensor1<double,2> &v,
                    const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
                    const FTensor::Tensor2<double,2,2> &nt_to_xy,
                    FTensor::Tensor1<double,2> &z_jump,
                    FTensor::Tensor2<double,2,2> &dz_jump);
@@ -60,14 +59,14 @@ void compute_Cp(const double zx[Nx+1][Ny
     for(int pm=0;pm<2;++pm)
       if(intersects[dd][pm])
         {
-          FTensor::Tensor1<double,2> v, dv, ddv;
+          FTensor::Tensor1<double,2> v, dv;
           FTensor::Tensor1<int,2> dir(0,0);
           FTensor::Tensor2<double,2,2> nt_to_xy;
           const int sgn(sign(pos(dd)-interface_pos[dd][pm](dd)));
           dir(dd)=sgn;
           compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
                                  dist_cell,dist_edge,
-                                 interface_pos[dd][pm],dd,v,dv,ddv,nt_to_xy);
+                                 interface_pos[dd][pm],dd,v,dv,nt_to_xy);
 
           FTensor::Tensor1<double,2> z_jump;
           FTensor::Tensor2<double,2,2> dz_jump;
@@ -85,7 +84,7 @@ void compute_Cp(const double zx[Nx+1][Ny
               abort();
               break;
             }
-          compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
+          compute_jumps(eta_jump,v,dv,nt_to_xy,z_jump,dz_jump);
 
           double dx=std::fabs(pos(dd)-interface_pos[dd][pm](dd) - (h/2)*sgn);
 
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_Cxyz.cxx
--- a/compute_coefficients/compute_Cxyz.cxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_Cxyz.cxx	Sun May 06 18:49:35 2012 -0700
@@ -7,22 +7,17 @@ void compute_jumps(const double &eta_jum
 void compute_jumps(const double &eta_jump,
                    const FTensor::Tensor1<double,2> &v,
                    const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
                    const FTensor::Tensor2<double,2,2> &nt_to_xy,
                    FTensor::Tensor1<double,2> &z_jump,
-                   FTensor::Tensor2<double,2,2> &dz_jump,
-                   FTensor::Tensor3_christof<double,2,2> &ddz_jump);
+                   FTensor::Tensor2<double,2,2> &dz_jump);
 
 void compute_jumps(const double &eta_jump,
                    const FTensor::Tensor1<double,2> &v,
                    const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
                    const FTensor::Tensor2<double,2,2> &nt_to_xy,
                    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);
+                   double &p_jump);
 
 void compute_Cxyz(const double zx[Nx+1][Ny],
                   const double zy[Nx][Ny+1],
@@ -63,73 +58,46 @@ void compute_Cxyz(const double zx[Nx+1][
                  dir=direction of 2nd derivative
                  dir2, dir3=components of mixed derivative we are correcting.
                  So for zy,xy with boundary in x, dir2=y, dir3=y. */
-              FTensor::Tensor1<double,2> v, dv, ddv;
+              FTensor::Tensor1<double,2> v, dv;
               FTensor::Tensor2<double,2,2> nt_to_xy;
               compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
                                      dist_cell,dist_edge,
                                      interface.dd_pos[dd][pm],d,
-                                     v,dv,ddv,nt_to_xy);
+                                     v,dv,nt_to_xy);
 
               int sgn(sign(pos(dd) - interface.dd_pos[dd][pm](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;
 
-
-              compute_jumps(interface.eta_jump_dd[dd][pm],v,dv,ddv,nt_to_xy,
-                            z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+              compute_jumps(interface.eta_jump_dd[dd][pm],v,dv,nt_to_xy,
+                            z_jump,dz_jump,p_jump);
 
               FTensor::Tensor1<double,2> dx(0,0);
               dx(dd)=pos(dd) - interface.dd_pos[dd][pm](dd) - h*sgn;
               double dz_dd_correction=
-                xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
-                        + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+                xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b))/(h*h);
 
               const int dz_dd_factor(d==dd ? 2 : 1);
 
               C+=dz_dd_factor*dz_dd_correction;
-
-              // if(d==0 && dd==1 && i==9 && j==2)
-              //   std::cout << __func__ << " "
-              //             << d << " "
-              //             // << pm << " "
-              //             << i << " "
-              //             << j << " "
-              //             // << dx(0) << " "
-              //             // << dx(1) << " "
-              //             // << interface.dd_pos[dd][pm](0) << " "
-              //             // << interface.dd_pos[dd][pm](1) << " "
-              //             // << dv(0) << " "
-              //             // << dv(1) << " "
-              //             // << ddv(0) << " "
-              //             // << ddv(1) << " "
-              //             << interface.eta_jump_dd[dd][pm] << " "
-              //             // << xyz(a)*(z_jump(a))/(h*h) << " "
-              //             // << xyz(a)*(dz_jump(a,b)*dx(b))/(h*h) << " "
-              //             << C << " "
-              //             << "\n";
-
             }
           if(interface.intersect_sides[dd][pm])
             {
-              FTensor::Tensor1<double,2> v, dv, ddv;
+              FTensor::Tensor1<double,2> v, dv;
               FTensor::Tensor2<double,2,2> nt_to_xy;
               compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
                                      dist_cell,dist_edge,
                                      interface.sides_pos[dd][pm],d,
-                                     v,dv,ddv,nt_to_xy);
+                                     v,dv,nt_to_xy);
 
               int sgn(sign(pos(dd) - interface.sides_pos[dd][pm](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;
 
-              compute_jumps(interface.eta_jump_sides[dd][pm],v,dv,ddv,nt_to_xy,
-                            z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+              compute_jumps(interface.eta_jump_sides[dd][pm],v,dv,nt_to_xy,
+                            z_jump,dz_jump,p_jump);
 
               try_corners=false;
               FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
@@ -139,33 +107,30 @@ void compute_Cxyz(const double zx[Nx+1][
               FTensor::Tensor1<double,2> dx(0,0);
               dx(dd)=pos(dd) - interface.sides_pos[dd][pm](dd) - (h/2)*sgn;
 
-              C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b)
-                       + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
+              C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b))/h;
             }
         }
       if(interface.intersect_dp[pm])
         {
-          FTensor::Tensor1<double,2> v, dv, ddv;
+          FTensor::Tensor1<double,2> v, dv;
           FTensor::Tensor2<double,2,2> nt_to_xy;
           compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
                                  dist_cell,dist_edge,
                                  interface.dp_pos[pm],d,
-                                 v,dv,ddv,nt_to_xy);
+                                 v,dv,nt_to_xy);
 
           int sgn(sign(pos(d) - interface.dp_pos[pm](d)));
           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;
 
-          compute_jumps(interface.eta_jump_dp[pm],v,dv,ddv,nt_to_xy,
-                        z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+          compute_jumps(interface.eta_jump_dp[pm],v,dv,nt_to_xy,
+                        z_jump,dz_jump,p_jump);
 
           FTensor::Tensor1<double,2> dx(0,0);
           dx(d)=pos(d) - interface.dp_pos[pm](d) - (h/2)*sgn;
 
-          C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
+          C+=sgn*p_jump/h;
         }
     }
 
@@ -188,19 +153,18 @@ void compute_Cxyz(const double zx[Nx+1][
         for(int ee=0;ee<2;++ee)
           if(interface.intersect_corner[dd][ee])
             {
-              FTensor::Tensor1<double,2> v, dv, ddv;
+              FTensor::Tensor1<double,2> v, dv;
               FTensor::Tensor2<double,2,2> nt_to_xy;
               compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
                                      dist_cell,dist_edge,
                                      interface.corner_pos[dd][ee],d,
-                                     v,dv,ddv,nt_to_xy);
+                                     v,dv,nt_to_xy);
 
               FTensor::Tensor1<double,2> z_jump;
               FTensor::Tensor2<double,2,2> dz_jump;
-              FTensor::Tensor3_christof<double,2,2> ddz_jump;
 
-              compute_jumps(interface.eta_jump_corner[dd][ee],v,dv,ddv,nt_to_xy,
-                            z_jump,dz_jump,ddz_jump);
+              compute_jumps(interface.eta_jump_corner[dd][ee],v,dv,nt_to_xy,
+                            z_jump,dz_jump);
 
               FTensor::Tensor1<double,2> dx,yx_pos;
 
@@ -210,8 +174,7 @@ void compute_Cxyz(const double zx[Nx+1][
               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);
+                zyx(a)*(z_jump(a) + dz_jump(a,b)*dx(b))/(h*h);
             }
     }
 }
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_jumps.cxx
--- a/compute_coefficients/compute_jumps.cxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_jumps.cxx	Sun May 06 18:49:35 2012 -0700
@@ -4,13 +4,10 @@ void compute_jumps(const double &eta_jum
 void compute_jumps(const double &eta_jump,
                    const FTensor::Tensor1<double,2> &v,
                    const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
                    const FTensor::Tensor2<double,2,2> &nt_to_xy,
                    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)
+                   double &p_jump)
 {
   const FTensor::Index<'a',2> a;
   const FTensor::Index<'b',2> b;
@@ -31,51 +28,15 @@ void compute_jumps(const double &eta_jum
   dz_jump(a,b)=dz_jump_nt(c,d)*nt_to_xy(a,c)*nt_to_xy(b,d);
 
   p_jump=-2*eta_jump*dv(t);
-
-  /* TODO: dp_jump and ddz_jump both need to be corrected for
-     curvature and density jumps. */
-
-  FTensor::Tensor1<double,2> dp_jump_nt;
-  dp_jump_nt(n)=2*eta_jump*ddv(n);
-  dp_jump_nt(t)=-2*eta_jump*ddv(t);
-  dp_jump(a)=nt_to_xy(a,b)*dp_jump_nt(b);
-
-  FTensor::Tensor3_christof<double,2,2> ddz_jump_nt;
-  ddz_jump_nt(a,t,t)=eta_jump*ddv(a);
-  ddz_jump_nt(a,n,n)=-ddz_jump_nt(a,t,t) + dp_jump_nt(a);
-  ddz_jump_nt(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
-
-  /* It would be nice if there were a way to do this directly instead
-     of iterating over all of the indices. */
-  ddz_jump(a,n,n)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(n,f);
-  ddz_jump(a,n,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(t,f);
-  ddz_jump(a,t,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(t,e)*nt_to_xy(t,f);
 }
 
 void compute_jumps(const double &eta_jump,
                    const FTensor::Tensor1<double,2> &v,
                    const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
-                   const FTensor::Tensor2<double,2,2> &nt_to_xy,
-                   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;
-  compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
-                p_jump,dp_jump);
-}
-
-void compute_jumps(const double &eta_jump,
-                   const FTensor::Tensor1<double,2> &v,
-                   const FTensor::Tensor1<double,2> &dv,
-                   const FTensor::Tensor1<double,2> &ddv, 
                    const FTensor::Tensor2<double,2,2> &nt_to_xy,
                    FTensor::Tensor1<double,2> &z_jump,
                    FTensor::Tensor2<double,2,2> &dz_jump)
 {
-  FTensor::Tensor3_christof<double,2,2> ddz_jump;
-  compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
+  double p_jump;
+  compute_jumps(eta_jump,v,dv,nt_to_xy,z_jump,dz_jump,p_jump);
 }
-
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_v_on_interface.hxx
--- a/compute_coefficients/compute_v_on_interface.hxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface.hxx	Sun May 06 18:49:35 2012 -0700
@@ -15,7 +15,6 @@ void compute_v_on_interface(const double
                             const int &xyz,
                             FTensor::Tensor1<double,2> &v,
                             FTensor::Tensor1<double,2> &dv,
-                            FTensor::Tensor1<double,2> &ddv,
                             FTensor::Tensor2<double,2,2> &nt_to_xy);
 
 #endif
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx	Fri Apr 13 16:36:09 2012 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-#ifndef GAMR_COMPUTE_2ND_DERIVS_HXX
-#define GAMR_COMPUTE_2ND_DERIVS_HXX
-
-#include "Valid.hxx"
-
-template<int ny>
-void compute_2nd_derivs(const int &d, const Valid &v,
-                        const int &d0, const int &d1,
-                        const int &nx,
-                        const double z[][ny],
-                        const double log_eta[][ny],
-                        FTensor::Tensor3_christof<double,2,2> &ddv_pm,
-                        double &eta_pm)
-{
-  /* Mixed derivative */
-  if(d0!=d1)
-    {
-      if(d==0)
-        {
-          if(v.j==0 || v.j==ny)
-            ddv_pm(d,d0,d1)=0;
-          else
-            ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
-              - vel(z,log_eta,v.i,v.j)
-              - vel(z,log_eta,v.i+1,v.j-1)
-              + vel(z,log_eta,v.i,v.j-1);
-        }
-      else
-        {
-          if(v.i==0 || v.i==nx)
-            ddv_pm(d,d0,d1)=0;
-          else
-            ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
-              - vel(z,log_eta,v.i,v.j)
-              - vel(z,log_eta,v.i-1,v.j+1)
-              + vel(z,log_eta,v.i-1,v.j);
-        }
-    }
-  else
-    {
-      switch(d0)
-        {
-        case 0:
-          if(v.i==0)
-            {
-              ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
-                - vel(z,log_eta,v.i,v.j);
-            }
-          else if(v.i==nx-1)
-            {
-              ddv_pm(d,d0,d1)=
-                - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
-            }
-          else
-            {
-              ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
-                - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
-            }
-          break;
-        case 1:
-          if(v.j==0)
-            {
-              ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
-                - vel(z,log_eta,v.i,v.j);
-            }
-          else if(v.j==ny-1)
-            {
-              ddv_pm(d,d0,d1)=
-                - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
-            }
-          else
-            {
-              ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
-                - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
-            }
-          break;
-        default:
-          abort();
-        }
-    }
-  eta_pm=exp(log_eta[v.i][v.j]);
-}
-
-#endif
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx	Sun May 06 18:49:35 2012 -0700
@@ -4,7 +4,6 @@
 #include <vector>
 #include "../sign.hxx"
 #include "vel.hxx"
-#include "compute_2nd_derivs.hxx"
 #include "compute_1st_derivs.hxx"
 #include "compute_values.hxx"
 #include <tuple>
@@ -25,7 +24,6 @@ void compute_dv_dtt(const double zx[Nx+1
                     const FTensor::Tensor1<double,2> &norm,
                     const FTensor::Tensor1<double,2> &tangent,
                     const FTensor::Tensor2<double,2,2> &nt_to_xy,
-                    FTensor::Tensor1<double,2> &ddv,
                     FTensor::Tensor1<double,2> &dv,
                     FTensor::Tensor1<double,2> &v)
 {
@@ -34,10 +32,6 @@ void compute_dv_dtt(const double zx[Nx+1
   const FTensor::Index<'a',2> a;
   const FTensor::Index<'b',2> b;
   const FTensor::Index<'c',2> c;
-  FTensor::Tensor3_christof<double,2,2> ddv_pm[2];
-  ddv_pm[0](a,b,c)=0;
-  ddv_pm[1](a,b,c)=0;
-
   FTensor::Tensor2<double,2,2> dv_pm[2];
   dv_pm[0](a,b)=0;
   dv_pm[1](a,b)=0;
@@ -371,31 +365,15 @@ void compute_dv_dtt(const double zx[Nx+1
 
                           is_set(pm,d0,d1)=true;
                           if(d==0)
-                            compute_2nd_derivs(0,V,d0,d1,Nx+1,zx,log_etax,
-                                               ddv_pm[pm],eta_pm[pm]);
+                            eta_pm[pm]=exp(log_etax[V.i][V.j]);
                           else
-                            compute_2nd_derivs(1,V,d0,d1,Nx,zy,log_etay,
-                                               ddv_pm[pm],eta_pm[pm]);
+                            eta_pm[pm]=exp(log_etay[V.i][V.j]);
                         }
                     }
                 }
-              for(int pm=0;pm<2;++pm)
-                ddv_pm[pm](d,d0,d1)/=num[pm];
             }
         }
     }
-
-  /* ddv */
-  FTensor::Tensor1<double,2> ddv_xy(0,0);
-  double temp(0);
-  for(int d=0;d<2;++d)
-    {
-      temp+=eta_pm[d]*eta_pm[d];
-      ddv_xy(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
-    }
-  ddv_xy(a)/=temp*h*h;
-
-  ddv(a)=nt_to_xy(b,a)*ddv_xy(b);
 
   /* dv */
   FTensor::Tensor1<double,2> dv_xy(0,0);
@@ -404,13 +382,6 @@ void compute_dv_dtt(const double zx[Nx+1
 
   dv(a)=nt_to_xy(b,a)*dv_xy(b);
 
-  FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
-  double eta_jump(eta_pm[1]-eta_pm[0]);
-  // /* TODO: fix for curved surfaces */
-  FTensor::Tensor1<double,2> ddz_nt_jump;
-  ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
-
-  dv(a)-=h*length*ddz_nt_jump(a);
   dv(a)/=(eta_pm[0]+eta_pm[1])*h;
 
   /* v */
@@ -420,7 +391,8 @@ void compute_dv_dtt(const double zx[Nx+1
 
   v(a)=nt_to_xy(b,a)*v_xy(b);
 
-  /* TODO: fix for curved surfaces */
+  FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+  double eta_jump(eta_pm[1]-eta_pm[0]);
   FTensor::Tensor1<double,2> dz_n_jump;
   dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
 
diff -r 3e5274e3300c -r 980e18575427 compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx	Fri Apr 13 16:36:09 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx	Sun May 06 18:49:35 2012 -0700
@@ -23,7 +23,6 @@ void compute_dv_dtt(const double zx[Nx+1
                     const FTensor::Tensor1<double,2> &norm,
                     const FTensor::Tensor1<double,2> &tangent,
                     const FTensor::Tensor2<double,2,2> &nt_to_xy,
-                    FTensor::Tensor1<double,2> &ddv,
                     FTensor::Tensor1<double,2> &dv,
                     FTensor::Tensor1<double,2> &v);
 
@@ -39,7 +38,6 @@ void compute_v_on_interface(const double
                             const int &xyz,
                             FTensor::Tensor1<double,2> &v,
                             FTensor::Tensor1<double,2> &dv,
-                            FTensor::Tensor1<double,2> &ddv,
                             FTensor::Tensor2<double,2,2> &nt_to_xy)
 {
   const FTensor::Index<'a',2> a;
@@ -56,5 +54,5 @@ void compute_v_on_interface(const double
   nt_to_xy(a,t)=tangent(a);
 
   compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,nt_to_xy,
-                 ddv,dv,v);
+                 dv,v);
 }



More information about the CIG-COMMITS mailing list