[cig-commits] commit: Remove more second derivative stuff.
Mercurial
hg at geodynamics.org
Tue May 8 05:50:55 PDT 2012
changeset: 156:3949c83f0526
user: Walter Landry <wlandry at caltech.edu>
date: Sun May 06 19:00:59 2012 -0700
files: compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
description:
Remove more second derivative stuff.
diff -r 980e18575427 -r 3949c83f0526 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Sun May 06 18:49:35 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Sun May 06 19:00:59 2012 -0700
@@ -45,11 +45,6 @@ void compute_dv_dtt(const double zx[Nx+1
/* First find which points can give valid derivatives */
for(int d=0;d<2;++d)
{
- 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));
@@ -66,20 +61,12 @@ void compute_dv_dtt(const double zx[Nx+1
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), p_d(i*h+dx+h/2,j*h+dy),
- dd_diff, d_diff, diff;
- dd_diff(a)=p(a)-pos(a);
+ d_diff, diff;
+ diff(a)=p(a)-pos(a);
d_diff(a)=p_d(a)-pos(a);
if(d==0)
{
- dd_diff(a)-=length*norm(a)*sign(distx[i][j]);
- /* ddv */
- if(i>min_x && i<max_x)
- 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,dd_diff);
/* dv */
if(i<max_x)
{
@@ -90,28 +77,12 @@ void compute_dv_dtt(const double zx[Nx+1
i,j,d_diff);
}
/* v */
- diff(a)=dd_diff(a)-0.5*length*norm(a)*sign(distx[i][j]);
+ diff(a)-=1.5*length*norm(a)*sign(distx[i][j]);
valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),true,i,j,diff);
}
else
{
- dd_diff(a)-=length*norm(a)*sign(disty[i][j]);
- /* ddv */
- if(i>min_x)
- {
- bool v_dd;
- if(i==0)
- v_dd=(sign(disty[i][j])==sign(disty[i+1][j]));
- else if(i==max_x)
- v_dd=(sign(disty[i][j])==sign(disty[i-1][j]));
- else
- v_dd=(sign(disty[i][j])==sign(disty[i+1][j])
- && sign(disty[i][j])==sign(disty[i-1][j]));
-
- 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);
- }
/* dv */
bool v_d;
/* TODO: This bothers me a little. What if the
@@ -128,7 +99,7 @@ void compute_dv_dtt(const double zx[Nx+1
Valid(sign(disty[i][j]),v_d,i,j,d_diff);
/* v */
- diff(a)=dd_diff(a)-0.5*length*norm(a)*sign(disty[i][j]);
+ diff(a)-=1.5*length*norm(a)*sign(disty[i][j]);
if(i>min_x)
valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(disty[i][j]),true,i,j,diff);
@@ -141,28 +112,11 @@ void compute_dv_dtt(const double zx[Nx+1
++j)
{
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;
d_diff(a)=p_d(a)-pos(a);
if(d==0)
{
- /* ddv */
- dd_diff(a)-=length*norm(a)*sign(distx[i][j]);
- if(i>min_y)
- {
- bool v_dd;
- if(j==0)
- v_dd=(sign(distx[i][j])==sign(distx[i][j+1]));
- else if(j==max_y)
- v_dd=(sign(distx[i][j])==sign(distx[i][j-1]));
- else
- v_dd=(sign(distx[i][j])==sign(distx[i][j+1])
- && sign(distx[i][j])==sign(distx[i][j-1]));
- dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v_dd,i,j,dd_diff);
- }
-
/* dv */
bool v_d;
if(j==min_y || j==max_y)
@@ -176,15 +130,6 @@ void compute_dv_dtt(const double zx[Nx+1
}
else
{
- /* ddv */
- dd_diff(a)-=length*norm(a)*sign(disty[i][j]);
- if(j>min_y && 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);
-
/* dv */
if(j<max_y)
{
@@ -196,51 +141,6 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
}
-
- /* d/dxdy */
- for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,Nx+d-1);++i)
- for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,Ny-d);++j)
- {
- /* Switch dx and dy since the derivative of u_x ends up on
- the u_y grid and vice versa */
- FTensor::Tensor1<double,2> p_dd(i*h+dy,j*h+dx), dd_diff;
- dd_diff(a)=p_dd(a)-pos(a);
- bool v_dd;
- if(d==0)
- {
- if(j==0)
- v_dd=(sign(distx[i+1][j])==sign(distx[i][j]));
- else if(j==Ny)
- v_dd=(sign(distx[i+1][j-1])==sign(distx[i][j-1]));
- else
- v_dd=(sign(distx[i+1][j])==sign(distx[i][j])
- && sign(distx[i+1][j])==sign(distx[i+1][j-1])
- && sign(distx[i+1][j])==sign(distx[i][j-1]));
-
- dd_diff(a)-=length*norm(a)*sign(distx[i][j]);
- dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v_dd,i,j,dd_diff);
- }
- else
- {
- if(i==0)
- v_dd=(sign(disty[i][j+1])==sign(disty[i][j]));
- else if(i==Nx)
- v_dd=(sign(disty[i-1][j+1])==sign(disty[i-1][j]));
- else
- v_dd=sign(disty[i][j+1])==sign(disty[i][j])
- && sign(disty[i][j+1])==sign(disty[i-1][j+1])
- && sign(disty[i][j+1])==sign(disty[i-1][j]);
-
- dd_diff(a)-=length*norm(a)*sign(disty[i][j]);
- dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
- }
- }
-
- sort(dd_valid(0,0).begin(), dd_valid(0,0).end());
- sort(dd_valid(0,1).begin(), dd_valid(0,1).end());
- sort(dd_valid(1,1).begin(), dd_valid(1,1).end());
sort(d_valid(0).begin(), d_valid(0).end());
sort(d_valid(1).begin(), d_valid(1).end());
@@ -297,11 +197,19 @@ void compute_dv_dtt(const double zx[Nx+1
if(values[pm].size()==6)
{
if(d==0)
- compute_values(0,values[pm],Nx+1,zx,log_etax,
- norm,pm==0 ? -length : length,v_pm[pm]);
+ {
+ eta_pm[pm]=exp(log_etax[values[pm].begin()->i]
+ [values[pm].begin()->j]);
+ compute_values(0,values[pm],Nx+1,zx,log_etax,
+ norm,pm==0 ? -length : length,v_pm[pm]);
+ }
else
- compute_values(1,values[pm],Nx,zy,log_etay,
- norm,pm==0 ? -length : length,v_pm[pm]);
+ {
+ eta_pm[pm]=exp(log_etay[values[pm].begin()->i]
+ [values[pm].begin()->j]);
+ compute_values(1,values[pm],Nx,zy,log_etay,
+ norm,pm==0 ? -length : length,v_pm[pm]);
+ }
}
break;
}
@@ -347,31 +255,6 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
}
- /* Second derivative ddv */
- for(int d1=d0;d1<2;++d1)
- {
- double distances[]={0,0};
- int num[]={0,0};
- for(auto &V: dd_valid(d0,d1))
- {
- for(int pm=0;pm<2;++pm)
- {
- if(V.sign==(pm==0 ? -1 : 1) && V.valid
- && (!is_set(pm,d0,d1)
- || about_equal(V.r,distances[pm])))
- {
- distances[pm]=V.r;
- ++num[pm];
-
- is_set(pm,d0,d1)=true;
- if(d==0)
- eta_pm[pm]=exp(log_etax[V.i][V.j]);
- else
- eta_pm[pm]=exp(log_etay[V.i][V.j]);
- }
- }
- }
- }
}
}
More information about the CIG-COMMITS
mailing list