[cig-commits] commit: Add some support for computing the first derivative. Also fix a minor
Mercurial
hg at geodynamics.org
Tue Mar 20 14:53:55 PDT 2012
changeset: 101:65c1c31f67b1
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 20 14:53:35 2012 -0700
files: compute_v_on_interface/compute_dv_dtt.cxx
description:
Add some support for computing the first derivative. Also fix a minor
issue that that it now uses mixed derivatives from the bottom (x=0 or
y=0).
diff -r b98a93816ad6 -r 65c1c31f67b1 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 20 13:41:38 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 20 14:53:35 2012 -0700
@@ -35,10 +35,14 @@ 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));
+ FTensor::Tensor2_symmetric<std::vector<Valid>,2> dd_valid;
+ dd_valid(0,0).resize((2*max_r+1)*(2*max_r+1));
+ dd_valid(0,1).resize((2*max_r+1)*(2*max_r+1));
+ dd_valid(1,1).resize((2*max_r+1)*(2*max_r+1));
+
+ FTensor::Tensor1<std::vector<Valid>,2> d_valid;
+ d_valid(0).resize((2*max_r+1)*(2*max_r+1));
+ d_valid(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);
@@ -49,31 +53,52 @@ void compute_dv_dtt(const double zx[Nx+1
for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,max_x);++i)
for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,max_y);++j)
{
- FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
- diff(a)=p(a)-pos(a);
+ FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), p_d(i*h+dx+h/2,j*h+dy),
+ dd_diff, d_diff;
+ dd_diff(a)=p(a)-pos(a);
+ d_diff(a)=p_d(a)-pos(a);
+
if(d==0)
{
if(i>0 && i<max_x)
- valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),
sign(distx[i][j])==sign(distx[i+1][j])
&& sign(distx[i][j])==sign(distx[i-1][j]),
- i,j,diff(a)*diff(a));
+ i,j,dd_diff(a)*dd_diff(a));
+ if(i<max_x)
+ d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),
+ sign(distx[i][j])==sign(distx[i+1][j]),
+ i,j,d_diff(a)*d_diff(a));
+
}
else
{
- bool v;
+ bool v_dd;
if(i==0)
- v=(sign(disty[i][j])==sign(disty[i+1][j]));
+ v_dd=(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_dd=(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_dd=(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));
+ dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),v_dd,i,j,dd_diff(a)*dd_diff(a));
+
+ bool v_d;
+ /* This bothers me a little. What if the interface
+ lies between the point and the boundary? Shouldn't
+ that make the derivative calculation invalid? */
+ if(i==max_x)
+ v_d=true;
+ else
+ v_d=(sign(disty[i][j])==sign(disty[i+1][j]));
+
+ d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),v_d,i,j,d_diff(a)*d_diff(a));
}
}
@@ -81,8 +106,10 @@ void compute_dv_dtt(const double zx[Nx+1
for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,max_x);++i)
for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,max_y);++j)
{
- FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
- diff(a)=p(a)-pos(a);
+ FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), p_d(i*h+dx,j*h+dy+h/2),
+ dd_diff, d_diff;
+ dd_diff(a)=p(a)-pos(a);
+ d_diff(a)=p_d(a)-pos(a);
if(d==0)
{
@@ -94,25 +121,40 @@ void compute_dv_dtt(const double zx[Nx+1
else
v=sign(distx[i][j])==sign(distx[i][j+1])
&& sign(distx[i][j])==sign(distx[i][j-1]);
- valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,dd_diff(a)*dd_diff(a));
+
+ bool v_d;
+ if(j==max_y)
+ v_d=true;
+ else
+ v_d=sign(distx[i][j])==sign(distx[i][j+1]);
+ d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v_d,i,j,d_diff(a)*d_diff(a));
}
- else if(j>0 && j<max_y)
- valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),sign(disty[i][j])==sign(disty[i][j+1])
- && sign(disty[i][j])==sign(disty[i][j-1]),
- i,j,diff(a)*diff(a));
+ else
+ {
+ if(j>0 && j<max_y)
+ dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),sign(disty[i][j])==sign(disty[i][j+1])
+ && sign(disty[i][j])==sign(disty[i][j-1]),
+ i,j,dd_diff(a)*dd_diff(a));
+ if(j<max_y)
+ d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),sign(disty[i][j])==sign(disty[i][j+1]),
+ i,j,d_diff(a)*d_diff(a));
+ }
}
/* d/dxdy */
for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,max_x);++i)
for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,max_y);++j)
{
- FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
- diff(a)=p(a)-pos(a);
+ FTensor::Tensor1<double,2> p_dd(i*h+dx+h/2,j*h+dy+h/2), dd_diff;
+ dd_diff(a)=p_dd(a)-pos(a);
if(d==0)
{
- if(i>0 && i<max_x)
+ if(i<max_x)
{
bool v;
if(j==max_y)
@@ -122,13 +164,13 @@ void compute_dv_dtt(const double zx[Nx+1
&& sign(distx[i][j])==sign(distx[i+1][j])
&& sign(distx[i][j])==sign(distx[i+1][j+1]);
- valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,dd_diff(a)*dd_diff(a));
}
}
else
{
- if(j>0 && j<max_y)
+ if(j<max_y)
{
bool v;
if(i==max_x)
@@ -138,22 +180,24 @@ void compute_dv_dtt(const double zx[Nx+1
&& sign(distx[i][j])==sign(distx[i+1][j])
&& sign(distx[i][j])==sign(distx[i+1][j+1]);
- valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,dd_diff(a)*dd_diff(a));
}
}
}
- std::sort(valid(0,0).begin(),valid(0,0).end());
- std::sort(valid(0,1).begin(),valid(0,1).end());
- std::sort(valid(1,1).begin(),valid(1,1).end());
+ std::sort(dd_valid(0,0).begin(),dd_valid(0,0).end());
+ std::sort(dd_valid(0,1).begin(),dd_valid(0,1).end());
+ std::sort(dd_valid(1,1).begin(),dd_valid(1,1).end());
+ std::sort(d_valid(0).begin(),d_valid(0).end());
+ std::sort(d_valid(1).begin(),d_valid(1).end());
FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
false,false);
for(int d0=0;d0<2;++d0)
for(int d1=d0;d1<2;++d1)
{
- for(auto &v: valid(d0,d1))
+ for(auto &v: dd_valid(d0,d1))
{
for(int pm=0;pm<2;++pm)
{
More information about the CIG-COMMITS
mailing list