[cig-commits] commit: Convert jumps from (norm, tangent) into xyz coordinates.
Mercurial
hg at geodynamics.org
Wed Mar 28 05:50:59 PDT 2012
changeset: 126:42d38edd5d9c
user: Walter Landry <wlandry at caltech.edu>
date: Wed Mar 28 04:49:04 2012 -0700
files: compute_Cp.cxx compute_Cxyz.cxx compute_jumps.cxx compute_v_on_interface.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx
description:
Convert jumps from (norm,tangent) into xyz coordinates.
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_Cp.cxx
--- a/compute_Cp.cxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_Cp.cxx Wed Mar 28 04:49:04 2012 -0700
@@ -8,6 +8,7 @@ void compute_jumps(const double &eta_jum
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);
@@ -44,9 +45,10 @@ void compute_Cp(const double zx[Nx+1][Ny
{
FTensor::Tensor1<double,2> v, dv, ddv;
FTensor::Tensor1<int,2> dir(0,0);
+ FTensor::Tensor2<double,2,2> nt_to_xy;
dir(dd)=1;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface_pos[dd],dd,v,dv,ddv);
+ interface_pos[dd],dd,v,dv,ddv,nt_to_xy);
FTensor::Tensor1<double,2> z_jump;
FTensor::Tensor2<double,2,2> dz_jump;
@@ -64,7 +66,7 @@ void compute_Cp(const double zx[Nx+1][Ny
abort();
break;
}
- compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump);
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
double dx=std::fabs(pos(dd)-interface_pos[dd](dd)
- (h/2)*sign(pos(dd)-interface_pos[dd](dd)));
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_Cxyz.cxx Wed Mar 28 04:49:04 2012 -0700
@@ -8,6 +8,7 @@ void compute_jumps(const double &eta_jum
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);
@@ -16,6 +17,7 @@ void compute_jumps(const double &eta_jum
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,
@@ -54,30 +56,31 @@ void compute_Cxyz(const double zx[Nx+1][
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::Tensor2<double,2,2> nt_to_xy;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
interface.pos[dd],d,
- v,dv,ddv);
+ v,dv,ddv,nt_to_xy);
+ int sgn(sign(pos(dd) - interface.pos[dd](dd)));
FTensor::Tensor1<double,2> dx(0,0);
- dx(dd)=pos(dd) - interface.pos[dd](dd)
- - h*sign(pos(dd) - interface.pos[dd](dd));
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - h*sgn;
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 eta_jump;
+ double eta_jump(sgn);
switch(d)
{
case 0:
switch(dd)
{
case 0:
- eta_jump=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
+ 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]);
+ eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
break;
default:
abort();
@@ -88,10 +91,10 @@ 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]);
+ 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]);
+ eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
break;
default:
abort();
@@ -103,22 +106,20 @@ void compute_Cxyz(const double zx[Nx+1][
break;
}
- compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump,
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
p_jump,dp_jump);
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);
- const int dz_dd_factor(-sign(dx(dd))*(d==dd ? 2 : 1));
-
+ const int dz_dd_factor(d==dd ? 2 : 1);
C+=dz_dd_factor*dz_dd_correction;
if(d==dd && interface.intersect_dp)
{
- dx(dd)=pos(dd) - interface.pos[dd](dd)
- - (h/2)*sign(pos(dd) - interface.pos[dd](dd));
- C+=(p_jump + dp_jump(a)*dx(a))/h;
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
+ C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
}
if(interface.intersect_sides[dd])
@@ -127,11 +128,10 @@ void compute_Cxyz(const double zx[Nx+1][
FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
dir2((d+1)%2)=1;
dir3((dd+1)%2)=1;
- dx(dd)=pos(dd) - interface.pos[dd](dd)
- - (h/2)*sign(pos(dd) - interface.pos[dd](dd));
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
- C+=-(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)
+ + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
}
}
}
@@ -151,9 +151,10 @@ void compute_Cxyz(const double zx[Nx+1][
if(interface.intersect_corner[dd][ee])
{
FTensor::Tensor1<double,2> v, dv, ddv;
+ FTensor::Tensor2<double,2,2> nt_to_xy;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
interface.corner_pos[dd][ee],d,
- v,dv,ddv);
+ v,dv,ddv,nt_to_xy);
FTensor::Tensor1<double,2> z_jump;
FTensor::Tensor2<double,2,2> dz_jump;
@@ -173,7 +174,7 @@ void compute_Cxyz(const double zx[Nx+1][
abort();
break;
}
- compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump);
+ 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)
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_jumps.cxx
--- a/compute_jumps.cxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_jumps.cxx Wed Mar 28 04:49:04 2012 -0700
@@ -5,6 +5,7 @@ void compute_jumps(const double &eta_jum
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,
@@ -13,50 +14,68 @@ void compute_jumps(const double &eta_jum
{
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ const FTensor::Index<'d',2> d;
+ const FTensor::Index<'e',2> e;
+ const FTensor::Index<'f',2> f;
- z_jump(a)=eta_jump*v(a);
+ z_jump(a)=eta_jump*v(b)*nt_to_xy(a,b);
FTensor::Number<0> n;
FTensor::Number<1> t;
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
-
- dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
- dz_jump(a,t)=eta_jump*dv(a);
+
+ FTensor::Tensor2<double,2,2> dz_jump_nt;
+ dz_jump_nt(a,n)=-eta_jump*dv(b)*ddel(a,b);
+ dz_jump_nt(a,t)=eta_jump*dv(a);
+ 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. */
- dp_jump(n)=2*eta_jump*ddv(n);
- dp_jump(t)=-2*eta_jump*ddv(t);
+ 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);
- ddz_jump(a,t,t)=eta_jump*ddv(a);
- ddz_jump(a,n,n)=-ddz_jump(a,t,t) + dp_jump(a);
- ddz_jump(a,n,t)=-eta_jump*ddv(b)*ddel(a,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(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,z_jump,dz_jump,ddz_jump,p_jump,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,z_jump,dz_jump,ddz_jump);
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
}
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_v_on_interface.hxx Wed Mar 28 04:49:04 2012 -0700
@@ -13,6 +13,7 @@ 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::Tensor1<double,2> &ddv,
+ FTensor::Tensor2<double,2,2> &nt_to_xy);
#endif
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 28 04:49:04 2012 -0700
@@ -20,6 +20,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 FTensor::Tensor2<double,2,2> &nt_to_xy,
FTensor::Tensor1<double,2> &ddv,
FTensor::Tensor1<double,2> &dv,
FTensor::Tensor1<double,2> &v)
@@ -351,25 +352,21 @@ void compute_dv_dtt(const double zx[Nx+1
/* ddv */
double temp(0);
FTensor::Tensor1<double,2> ddv_xy(0,0);
- FTensor::Number<0> n;
- FTensor::Number<1> t;
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(n)=ddv_xy(a)*norm(a);
- ddv(t)=ddv_xy(a)*tangent(a);
+
+ ddv(a)=nt_to_xy(b,a)*ddv_xy(b);
/* dv */
FTensor::Tensor1<double,2> dv_xy(0,0);
for(int d=0;d<2;++d)
dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
- dv(n)=dv_xy(a)*norm(a);
- dv(t)=dv_xy(a)*tangent(a);
-
+ dv(a)=nt_to_xy(b,a)*dv_xy(b);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
FTensor::Tensor1<double,2> ddz_nt_jump;
@@ -384,9 +381,7 @@ void compute_dv_dtt(const double zx[Nx+1
for(int d=0;d<2;++d)
v_xy(a)+=v_pm[d](a)*eta_pm[d];
- v(n)=v_xy(a)*norm(a);
- v(t)=v_xy(a)*tangent(a);
-
+ v(a)=nt_to_xy(b,a)*v_xy(b);
FTensor::Tensor1<double,2> dz_n_jump;
dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
diff -r 045f5d7030f2 -r 42d38edd5d9c compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 28 02:44:07 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 28 04:49:04 2012 -0700
@@ -16,6 +16,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 FTensor::Tensor2<double,2,2> &nt_to_xy,
FTensor::Tensor1<double,2> &ddv,
FTensor::Tensor1<double,2> &dv,
FTensor::Tensor1<double,2> &v);
@@ -30,7 +31,8 @@ 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::Tensor1<double,2> &ddv,
+ FTensor::Tensor2<double,2,2> &nt_to_xy)
{
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
@@ -54,6 +56,11 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
+ FTensor::Number<0> n;
+ FTensor::Number<1> t;
+ nt_to_xy(a,n)=norm(a);
+ 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);
}
More information about the CIG-COMMITS
mailing list