[cig-commits] commit: Remove (i, j) from the call to compute_v_on_interface since it is not
Mercurial
hg at geodynamics.org
Fri Mar 23 09:50:47 PDT 2012
changeset: 113:644a357572e3
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 23 09:50:30 2012 -0700
files: Interface.hxx Interface/Interface.cxx compute_Cp.cxx 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 (i,j) from the call to compute_v_on_interface since it is not
used. Added the beginning of support for corner cutting. Also added
notes about handling curved interfaces.
diff -r f80876a0a634 -r 644a357572e3 Interface.hxx
--- a/Interface.hxx Thu Mar 22 15:16:26 2012 -0700
+++ b/Interface.hxx Fri Mar 23 09:50:30 2012 -0700
@@ -8,7 +8,8 @@ class Interface
class Interface
{
public:
- FTensor::Tensor1<double,2> grid_pos, pos[2], corner_pos[2][2][2];
+ FTensor::Tensor1<double,2> grid_pos, pos[2], corner_pos[2][2];
+ int corner_dir;
bool intersect_dd[2], intersect_xy[2], intersect_corner[2][2];
bool intersect_dp;
@@ -41,13 +42,28 @@ public:
{
double dh[2]={sign(n0-0.5)*h/2, sign(n1-0.5)*h/2};
double delta;
- delta=dist[i+n0][j+(n1+1)%2]-dist[i+n0][j+n1];
- corner_pos[n0][n1][0](0)=grid_pos(0)+dh[0]-dist[i+n0][j+n1]*h/delta;
- corner_pos[n0][n1][0](1)=grid_pos(1)+dh[1];
- delta=dist[i][j+1]-dist[i][j];
- corner_pos[n0][n1][1](0)=grid_pos(0)+dh[0];
- corner_pos[n0][n1][1](1)=grid_pos(1)+dh[1]-dist[i+n0][j+n1]*h/delta;
+ /* TODO fix this for higher order curved surfaces */
+
+ /* Correct in the direction that is closest to the
+ corner. This should use the interface that is more
+ perpendicular to the derivative direction. */
+ if(std::fabs(dist[i+(n0+1)%2][j+n1])>std::fabs(dist[i+n0][j+(n1+1)%2]))
+ {
+ corner_dir=0;
+ delta=dist[i+n0][j+(n1+1)%2]-dist[i+n0][j+n1];
+ corner_pos[n0][n1](0)=grid_pos(0)+dh[0]-dist[i+n0][j+n1]*h/delta;
+ corner_pos[n0][n1](1)=grid_pos(1)+dh[1];
+ }
+ else
+ {
+ corner_dir=1;
+ delta=dist[i+n0][j+(n1+1)%2]-dist[i+n0][j+n1];
+ corner_pos[n0][n1](0)=grid_pos(0)+dh[0];
+ corner_pos[n0][n1](1)=grid_pos(1)+dh[1]-dist[i+n0][j+n1]*h/delta;
+ }
+
+
}
}
}
diff -r f80876a0a634 -r 644a357572e3 Interface/Interface.cxx
--- a/Interface/Interface.cxx Thu Mar 22 15:16:26 2012 -0700
+++ b/Interface/Interface.cxx Fri Mar 23 09:50:30 2012 -0700
@@ -36,6 +36,8 @@ Interface::Interface(const int &d, const
intersect_dp=(sign(disty[i][j] + disty[i][j+1])
!=sign(disty[i-1][j] + disty[i-1][j+1]));
+ /* TODO This will have an error in the position of the interface
+ of O(h^2). Need to correctt for curvature. */
double delta;
if(intersect_dd[0])
{
diff -r f80876a0a634 -r 644a357572e3 compute_Cp.cxx
--- a/compute_Cp.cxx Thu Mar 22 15:16:26 2012 -0700
+++ b/compute_Cp.cxx Fri Mar 23 09:50:30 2012 -0700
@@ -19,6 +19,7 @@ void compute_Cp(const double zx[Nx+1][Ny
FTensor::Tensor1<double,2> pos((i+0.5)*h,(j+0.5)*h);
+ /* TODO fix for curved interfaces */
FTensor::Tensor1<double,2> interface_pos[2];
if(intersects[0])
{
@@ -38,7 +39,7 @@ void compute_Cp(const double zx[Nx+1][Ny
FTensor::Tensor1<int,2> dir(0,0);
dir(dd)=1;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface_pos[dd],i,j,dd,v,dv);
+ interface_pos[dd],dd,v,dv);
double dx=(pos(dd)>interface_pos[dd](dd))
? (interface_pos[dd](dd)-(pos(dd)-h/2))
diff -r f80876a0a634 -r 644a357572e3 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Thu Mar 22 15:16:26 2012 -0700
+++ b/compute_Cxyz.cxx Fri Mar 23 09:50:30 2012 -0700
@@ -37,7 +37,7 @@ void compute_Cxyz(const double zx[Nx+1][
xyz(d)=1;
dir(dd)=1;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface.pos[dd],i,j,d,
+ interface.pos[dd],d,
v,dv,ddv);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
@@ -86,7 +86,6 @@ void compute_Cxyz(const double zx[Nx+1][
C+=(p_jump + dx*dp_jump(a)*dir(a))/h;
}
- /* TODO handle corner cutting */
if(interface.intersect_xy[dd])
{
FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
@@ -100,5 +99,19 @@ void compute_Cxyz(const double zx[Nx+1][
+ dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
}
}
+ // /* Handle corner cutting */
+ // for(int ee=0;ee<2;++ee)
+ // {
+ // if(interface.intersect_corner[dd][ee])
+ // {
+ // FTensor::Tensor1<double,2> v, dv, ddv;
+ // compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ // interface.corner_pos[dd][ee],d,
+ // v,dv,ddv);
+
+ // double dh[2]={sign(dd-0.5)*h/2, sign(ee-0.5)*h/2};
+
+ // }
+ // }
}
}
diff -r f80876a0a634 -r 644a357572e3 compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Thu Mar 22 15:16:26 2012 -0700
+++ b/compute_v_on_interface.hxx Fri Mar 23 09:50:30 2012 -0700
@@ -10,7 +10,6 @@ void compute_v_on_interface(const double
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
const FTensor::Tensor1<double,2> &pos,
- const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv);
@@ -22,7 +21,6 @@ void compute_v_on_interface(const double
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
const FTensor::Tensor1<double,2> &pos,
- const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
diff -r f80876a0a634 -r 644a357572e3 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Thu Mar 22 15:16:26 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Fri Mar 23 09:50:30 2012 -0700
@@ -20,7 +20,6 @@ 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> &dv,
FTensor::Tensor1<double,2> &v)
diff -r f80876a0a634 -r 644a357572e3 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Thu Mar 22 15:16:26 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Fri Mar 23 09:50:30 2012 -0700
@@ -16,7 +16,6 @@ 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> &dv,
FTensor::Tensor1<double,2> &v);
@@ -28,13 +27,12 @@ void compute_v_on_interface(const double
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
const FTensor::Tensor1<double,2> &pos,
- const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &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);
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,pos,xyz,v,dv,ddv);
}
void compute_v_on_interface(const double zx[Nx+1][Ny],
@@ -44,7 +42,6 @@ void compute_v_on_interface(const double
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
const FTensor::Tensor1<double,2> &pos,
- const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
@@ -53,39 +50,31 @@ void compute_v_on_interface(const double
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
int ix(pos(0)/h), iy(pos(0)/h - 0.5), jx(pos(1)/h - 0.5), jy(pos(1)/h);
- double dx_i((pos(0)-ix*h)/h), dx_j((pos(0)-iy*h)/h-0.5),
- dy_i((pos(1)-jx*h)/h-0.5), dy_j((pos(1)-jy*h)/h);
+ FTensor::Tensor2<double,2,2> dx((pos(0)-ix*h)/h,(pos(0)-iy*h)/h-0.5,
+ (pos(1)-jx*h)/h-0.5,(pos(1)-jy*h)/h);
+
+ /* TODO: fix for curved interfaces */
+ FTensor::Tensor1<double,2> norm, tangent;
if(xyz==0)
{
- 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);
- norm(1)=((disty[iy][jy+1] - disty[iy][jy])*(h-dx_j)
- + (disty[iy+1][jy+1] - disty[iy+1][jy])*dx_j)/(h*h);
-
- norm(a)/=sqrt(norm(b)*norm(b));
-
- tangent(0)=-norm(1);
- tangent(1)=norm(0);
-
- compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
- i,j,ddv,dv,v);
+ norm(0)=((disty[iy+1][jy] - disty[iy][jy])*(1-dx(0,xyz))
+ + (disty[iy+1][jy+1] - disty[iy][jy+1])*dx(0,xyz))/h;
+ norm(1)=((disty[iy][jy+1] - disty[iy][jy])*(1-dx(1,xyz))
+ + (disty[iy+1][jy+1] - disty[iy+1][jy])*dx(1,xyz))/h;
}
else
{
- 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(0)=((distx[ix+1][jx] - distx[ix][jx])*(1-dx(0,xyz))
+ + (distx[ix+1][jx+1] - distx[ix][jx+1])*dx(0,xyz))/h;
+ norm(1)=((distx[ix][jx+1] - distx[ix][jx])*(1-dx(1,xyz))
+ + (distx[ix+1][jx+1] - distx[ix+1][jx])*dx(1,xyz))/h;
+ }
+ norm(a)/=sqrt(norm(b)*norm(b));
- norm(a)/=sqrt(norm(b)*norm(b));
+ tangent(0)=-norm(1);
+ tangent(1)=norm(0);
- tangent(0)=-norm(1);
- tangent(1)=norm(0);
-
- compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
- i,j,ddv,dv,v);
- }
+ compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
+ ddv,dv,v);
}
More information about the CIG-COMMITS
mailing list