[cig-commits] commit: Factor out computing the norm into compute_norm.
Mercurial
hg at geodynamics.org
Fri Apr 13 16:41:20 PDT 2012
changeset: 154:3e5274e3300c
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Apr 13 16:36:09 2012 -0700
files: compute_coefficients/compute_Cp.cxx compute_coefficients/compute_Cxyz.cxx compute_coefficients/compute_coefficients.cxx compute_coefficients/compute_v_on_interface.hxx compute_coefficients/compute_v_on_interface/compute_norm.cxx compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx compute_coefficients/simplified_Rp.cxx wscript
description:
Factor out computing the norm into compute_norm.
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_Cp.cxx
--- a/compute_coefficients/compute_Cp.cxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/compute_Cp.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -20,6 +20,7 @@ void compute_Cp(const double zx[Nx+1][Ny
const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const int &i, const int &j,
double &Cp)
{
@@ -65,6 +66,7 @@ void compute_Cp(const double zx[Nx+1][Ny
const int sgn(sign(pos(dd)-interface_pos[dd][pm](dd)));
dir(dd)=sgn;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ dist_cell,dist_edge,
interface_pos[dd][pm],dd,v,dv,ddv,nt_to_xy);
FTensor::Tensor1<double,2> z_jump;
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_Cxyz.cxx
--- a/compute_coefficients/compute_Cxyz.cxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/compute_Cxyz.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -66,6 +66,7 @@ void compute_Cxyz(const double zx[Nx+1][
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,
+ dist_cell,dist_edge,
interface.dd_pos[dd][pm],d,
v,dv,ddv,nt_to_xy);
@@ -116,6 +117,7 @@ void compute_Cxyz(const double zx[Nx+1][
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,
+ dist_cell,dist_edge,
interface.sides_pos[dd][pm],d,
v,dv,ddv,nt_to_xy);
@@ -146,6 +148,7 @@ void compute_Cxyz(const double zx[Nx+1][
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,
+ dist_cell,dist_edge,
interface.dp_pos[pm],d,
v,dv,ddv,nt_to_xy);
@@ -188,6 +191,7 @@ void compute_Cxyz(const double zx[Nx+1][
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,
+ dist_cell,dist_edge,
interface.corner_pos[dd][ee],d,
v,dv,ddv,nt_to_xy);
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_coefficients.cxx
--- a/compute_coefficients/compute_coefficients.cxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/compute_coefficients.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -48,7 +48,8 @@ double simplified_Rp(const int &i, const
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
- const double dist_cell[Nx][Ny]);
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1]);
void compute_coefficients(const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
@@ -171,7 +172,7 @@ void compute_coefficients(const double l
{
zx[di][dj]=1;
double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
- distx,disty,dist_cell);
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfp[0][i][j].push_back(Coefficient(di,dj,R));
zx[di][dj]=0;
@@ -183,7 +184,7 @@ void compute_coefficients(const double l
{
zy[di][dj]=1;
double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
- distx,disty,dist_cell);
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfp[1][i][j].push_back(Coefficient(di,dj,R));
zy[di][dj]=0;
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_v_on_interface.hxx
--- a/compute_coefficients/compute_v_on_interface.hxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface.hxx Fri Apr 13 16:36:09 2012 -0700
@@ -9,6 +9,8 @@ void compute_v_on_interface(const double
const double log_etay[Nx][Ny+1],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const FTensor::Tensor1<double,2> &pos,
const int &xyz,
FTensor::Tensor1<double,2> &v,
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_v_on_interface/compute_norm.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_norm.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -0,0 +1,30 @@
+#include "../constants.hxx"
+#include "FTensor.hpp"
+
+FTensor::Tensor1<double,2> compute_norm(const FTensor::Tensor1<double,2> &pos,
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
+{
+ 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);
+
+ /* TODO: fix for curved interfaces, especially when the interface
+ point is not along a grid coordinate. */
+ FTensor::Tensor1<double,2> norm;
+ if(iy<Ny)
+ {
+ norm(0)=(disty[iy+1][jy] - disty[iy][jy])/h;
+ norm(1)=(disty[iy][jy+1] - disty[iy][jy])/h;
+ }
+ else
+ {
+ norm(0)=(distx[ix+1][jx] - distx[ix][jx])/h;
+ norm(1)=(distx[ix][jx+1] - distx[ix][jx])/h;
+ }
+ norm(a)/=sqrt(norm(b)*norm(b));
+
+ return norm;
+}
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -6,6 +6,12 @@
#include "FTensor.hpp"
#include "../compute_v_on_interface.hxx"
#include "vel.hxx"
+
+FTensor::Tensor1<double,2> compute_norm(const FTensor::Tensor1<double,2> &pos,
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1]);
void compute_dv_dtt(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -27,6 +33,8 @@ void compute_v_on_interface(const double
const double log_etay[Nx][Ny+1],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const FTensor::Tensor1<double,2> &pos,
const int &xyz,
FTensor::Tensor1<double,2> &v,
@@ -35,24 +43,10 @@ void compute_v_on_interface(const double
FTensor::Tensor2<double,2,2> &nt_to_xy)
{
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);
- /* TODO: fix for curved interfaces, especially when the interface
- point is not along a grid coordinate. */
- FTensor::Tensor1<double,2> norm, tangent;
- if(xyz==0)
- {
- norm(0)=(disty[iy+1][jy] - disty[iy][jy])/h;
- norm(1)=(disty[iy][jy+1] - disty[iy][jy])/h;
- }
- else
- {
- norm(0)=(distx[ix+1][jx] - distx[ix][jx])/h;
- norm(1)=(distx[ix][jx+1] - distx[ix][jx])/h;
- }
- norm(a)/=sqrt(norm(b)*norm(b));
+ FTensor::Tensor1<double,2> norm(compute_norm(pos,distx,disty,dist_cell,dist_edge));
+ FTensor::Tensor1<double,2> tangent;
tangent(0)=-norm(1);
tangent(1)=norm(0);
diff -r dd661cbf8cc4 -r 3e5274e3300c compute_coefficients/simplified_Rp.cxx
--- a/compute_coefficients/simplified_Rp.cxx Fri Apr 13 16:34:40 2012 -0700
+++ b/compute_coefficients/simplified_Rp.cxx Fri Apr 13 16:36:09 2012 -0700
@@ -10,6 +10,7 @@ extern void compute_Cp(const double zx[N
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const int &i, const int &j,
double &Cp);
@@ -23,12 +24,13 @@ double simplified_Rp(const int &i, const
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
- const double dist_cell[Nx][Ny])
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
{
double dzx_x=(zx[i+1][j] - zx[i][j])/h;
double dzy_y=(zy[i][j+1] - zy[i][j])/h;
double Cp;
- compute_Cp(zx,zy,log_etax,log_etay,eta_cell,distx,disty,dist_cell,i,j,Cp);
+ compute_Cp(zx,zy,log_etax,log_etay,eta_cell,distx,disty,dist_cell,dist_edge,i,j,Cp);
return dzx_x + dzy_y + Cp;
}
diff -r dd661cbf8cc4 -r 3e5274e3300c wscript
--- a/wscript Fri Apr 13 16:34:40 2012 -0700
+++ b/wscript Fri Apr 13 16:36:09 2012 -0700
@@ -15,6 +15,7 @@ def build(bld):
'compute_coefficients/simplified_Ry.cxx',
'compute_coefficients/simplified_Rp.cxx',
'compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx',
+ 'compute_coefficients/compute_v_on_interface/compute_norm.cxx',
'compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx',
'compute_coefficients/compute_jumps.cxx',
'compute_coefficients/Interface/Interface.cxx'],
More information about the CIG-COMMITS
mailing list