[cig-commits] commit: Remove all of the logic for computing dC_dv within compute_v_on_interface.
Mercurial
hg at geodynamics.org
Thu Mar 22 11:44:59 PDT 2012
changeset: 109:8a00c28704c1
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Thu Mar 22 11:44:33 2012 -0700
files: compute_Cxyz.cxx compute_v_on_interface.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx
description:
Remove all of the logic for computing dC_dv within compute_v_on_interface.
diff -r 115c2a1443ff -r 8a00c28704c1 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Thu Mar 22 11:36:26 2012 -0700
+++ b/compute_Cxyz.cxx Thu Mar 22 11:44:33 2012 -0700
@@ -14,7 +14,6 @@ void compute_Cxyz(const double zx[Nx+1][
double &C)
{
C=0;
- double dC_dz=0;
const Interface interface(d,i,j,distx,disty);
if(!interface.intersects())
@@ -25,7 +24,6 @@ void compute_Cxyz(const double zx[Nx+1][
const FTensor::Index<'c',2> c;
pos(a)=interface.grid_pos(a);
- double dC_dv(0);
for(int dd=0;dd<2;++dd)
{
if(interface.intersect_dd[dd])
@@ -34,13 +32,13 @@ 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, v_dv, dv_dv, ddv_dv;
+ FTensor::Tensor1<double,2> v, dv, ddv;
FTensor::Tensor1<int,2> xyz(0,0), dir(0,0);
xyz(d)=1;
dir(dd)=1;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
interface.pos[dd],i,j,d,
- v,dv,ddv,v_dv,dv_dv,ddv_dv);
+ v,dv,ddv);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
@@ -76,24 +74,6 @@ void compute_Cxyz(const double zx[Nx+1][
C+=dz_dd_factor*dz_dd_correction;
- FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
- -2*eta_jump*ddv_dv(t));
- double p_dv=-2*eta_jump*dv_dv(t);
-
- FTensor::Tensor2<double,2,2> dz_dv;
- dz_dv(a,n)=-eta_jump*dv_dv(b)*ddel(a,b);
- dz_dv(a,t)=eta_jump*dv_dv(a);
-
- FTensor::Tensor3_christof<double,2,2> ddz_dv;
- ddz_dv(a,t,t)=eta_jump*ddv_dv(a);
- ddz_dv(a,n,n)=-ddz_dv(a,t,t) + dp_dv(a);
- ddz_dv(a,n,t)=-eta_jump*ddv_dv(b)*ddel(a,b);
-
- dC_dv-=dz_dd_factor
- *(eta_jump*v_dv(a)*xyz(a)
- - dx*dz_dv(a,b)*xyz(a)*dir(b)
- + dx*dx*ddz_dv(a,b,c)*xyz(a)*dir(b)*dir(c)/2)/(h*h);
-
if(d==dd && interface.intersect_dp)
{
dx=(pos(dd)>interface.pos[dd](dd))
@@ -101,7 +81,6 @@ void compute_Cxyz(const double zx[Nx+1][
: ((pos(dd)+h/2)-interface.pos[dd](dd));
C+=(p_jump + dx*dp_jump(a)*dir(a))/h;
- dC_dv+=(p_dv + dx*dp_dv(a)*dir(a))/h;
}
/* TODO handle corner cutting */
@@ -116,22 +95,7 @@ void compute_Cxyz(const double zx[Nx+1][
C+=-(dz_jump(a,b)*dir2(a)*dir3(b)
+ dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
-
- dC_dv+=-(dz_dv(a,b)*dir2(a)*dir3(b)
- + dx*ddz_dv(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
}
}
}
-
- switch(d)
- {
- case 0:
- dC_dz+=dC_dv*exp(-log_etax[i][j]);
- break;
- case 1:
- dC_dz+=dC_dv*exp(-log_etay[i][j]);
- break;
- default:
- abort();
- }
}
diff -r 115c2a1443ff -r 8a00c28704c1 compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Thu Mar 22 11:36:26 2012 -0700
+++ b/compute_v_on_interface.hxx Thu Mar 22 11:44:33 2012 -0700
@@ -26,9 +26,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::Tensor1<double,2> &v_dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &ddv_dv);
+ FTensor::Tensor1<double,2> &ddv);
#endif
diff -r 115c2a1443ff -r 8a00c28704c1 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Thu Mar 22 11:36:26 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Thu Mar 22 11:44:33 2012 -0700
@@ -22,11 +22,8 @@ void compute_dv_dtt(const double zx[Nx+1
const FTensor::Tensor1<double,2> &tangent,
const int &ii, const int &jj,
FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &ddv_dv,
FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &v_dv)
+ FTensor::Tensor1<double,2> &v)
{
const int max_r=4;
diff -r 115c2a1443ff -r 8a00c28704c1 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Thu Mar 22 11:36:26 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Thu Mar 22 11:44:33 2012 -0700
@@ -18,11 +18,8 @@ void compute_dv_dtt(const double zx[Nx+1
const FTensor::Tensor1<double,2> &tangent,
const int &i, const int &j,
FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &ddv_dv,
FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &v_dv);
+ FTensor::Tensor1<double,2> &v);
void compute_v_on_interface(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -36,9 +33,8 @@ void compute_v_on_interface(const double
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv)
{
- FTensor::Tensor1<double,2> ddv, v_dv, dv_dv, ddv_dv;
- compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,pos,i,j,xyz,v,dv,ddv,
- v_dv,dv_dv,ddv_dv);
+ FTensor::Tensor1<double,2> ddv;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,pos,i,j,xyz,v,dv,ddv);
}
void compute_v_on_interface(const double zx[Nx+1][Ny],
@@ -52,10 +48,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> &v_dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &ddv_dv)
+ FTensor::Tensor1<double,2> &ddv)
{
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
@@ -78,33 +71,8 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- FTensor::Tensor1<double,2> ddv_dv_temp,dv_dv_temp,v_temp,v_dv_temp;
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
- i,j,ddv,ddv_dv_temp,dv,dv_dv_temp,v,v_dv_temp);
-
- double eta=exp(log_etax[i][j]);
- double x_i(i*h);
- double delta_x=(x_i>middle) ? (middle-(x_i-h))
- : (middle-(x_i+h));
- double delta=(h-std::fabs(delta_x))/h;
- double dvx_yy_dv=-2*eta*eta/(h*h*(min_eta*min_eta + max_eta*max_eta));
- if(j==0 || j==Ny-1)
- dvx_yy_dv/=2;
-
- double dzy_yx_dv=-eta_jump*dvx_yy_dv;
- double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
- double dzx_x_dv=eta_jump*dvy_y_dv;
- double dvp_dv=delta/2 + delta*delta/2;
- double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
- double dvx_dv=(eta*(4*dvp_dv - dvpp_dv) + 2*h*dzx_x_dv)
- /(3*(min_eta+max_eta));
-
- ddv_dv(0)=dvx_yy_dv;
- ddv_dv(1)=0;
- dv_dv(0)=0;
- dv_dv(1)=dvy_y_dv;
- v_dv(0)=dvx_dv;
- v_dv(1)=0;
+ i,j,ddv,dv,v);
}
else
{
@@ -121,30 +89,7 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- FTensor::Tensor1<double,2> ddv_dv_temp,dv_dv_temp,v_temp,v_dv_temp;
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
- i,j,ddv,ddv_dv_temp,dv,dv_dv_temp,v,v_dv_temp);
-
- double eta=exp(log_etay[i][j]);
- double x_i((i+0.5)*h);
- double delta_x=(x_i>middle) ? (middle-(x_i-h))
- : (middle-(x_i+h));
- double delta=(h-std::fabs(delta_x))/h;
- double dvy_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
- double dzx_yx_dv=-eta_jump*dvy_yy_dv;
- double dvx_y_dv=(-h/(min_eta+max_eta))*dzx_yx_dv;
- double dzy_x_dv=-eta_jump*dvx_y_dv;
- double dvp_dv=delta/2 + delta*delta/2;
- double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
- double dvy_dv=
- (eta*(4*dvp_dv - dvpp_dv) - 2*h*dzy_x_dv)/(3*(min_eta+max_eta));
-
- ddv_dv(0)=0;
- ddv_dv(1)=dvy_yy_dv;
- dv_dv(0)=dvx_y_dv;
- dv_dv(1)=0;
- v_dv(0)=0;
- v_dv(1)=dvy_dv;
-
+ i,j,ddv,dv,v);
}
}
More information about the CIG-COMMITS
mailing list