[cig-commits] commit: First efforts to generalize the interface code. Uses FTensor

Mercurial hg at geodynamics.org
Fri Mar 9 13:30:07 PST 2012


changeset:   75:7b5ea67a2e69
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Fri Mar 09 11:16:05 2012 -0800
files:       compute_Cx.cxx wscript
description:
First efforts to generalize the interface code.  Uses FTensor


diff -r 114d616b8587 -r 7b5ea67a2e69 compute_Cx.cxx
--- a/compute_Cx.cxx	Mon Mar 05 17:40:26 2012 -0800
+++ b/compute_Cx.cxx	Fri Mar 09 11:16:05 2012 -0800
@@ -2,6 +2,7 @@
 #include "Model.hxx"
 #include <iostream>
 #include "compute_v_on_interface.hxx"
+#include "FTensor.hpp"
 
 void compute_Cx(const Model &model, const double zx[Nx+1][Ny],
                 const double zy[Nx][Ny+1],
@@ -26,19 +27,34 @@ void compute_Cx(const Model &model, cons
                                  vx,vy,dvx_y,dvy_y,
                                  dvx_yy,dvy_yy);
 
+          const int d(0);
+          FTensor::Tensor1<double,2> v(vx,0), dv(0,dvy_y), ddv(dvx_yy,0), norm(1,0),
+            tangent(0,1);
+          FTensor::Tensor2_antisymmetric<double,2> eps(1);
+          const FTensor::Index<'a',2> a;
+
           double dx=(x>middle) ? (middle-(x-h))
             : (middle-(x+h));
-          double zx_jump=eta_jump*vx;
-          double dzx_x_jump=eta_jump*dvy_y;
-          double p_jump=-2*dvy_y*eta_jump;
-          double dp_x_jump=2*dvx_yy*eta_jump;
-          double dzx_xx_jump=
-            -dvx_yy*eta_jump + dp_x_jump;
-          double dzx_xx_correction=
-            -(zx_jump + dx*dzx_x_jump
-              + dx*dx*dzx_xx_jump/2)/(h*h);
-          const int dzx_xx_sign(x>middle ? -1 : 1);
-          Cx=dzx_xx_sign*2*dzx_xx_correction;
+          double z_jump=eta_jump*v(d);
+
+          double dz_n_jump=eta_jump*(dv(0)*norm(a) - dv(1)*tangent(a))*eps(d,a);
+          double dz_t_jump=eta_jump*(dv(0)*norm(d) + dv(1)*tangent(d));
+          double dz_x_jump=dz_n_jump*norm(d) + dz_t_jump*tangent(d);
+
+          double p_jump=-2*eta_jump*dv(1);
+
+          double dp_t_jump=-2*eta_jump*ddv(1);
+          double dp_n_jump=2*eta_jump*ddv(0);
+          double dp_x_jump=dp_n_jump*norm(d) + dp_t_jump*tangent(d);
+
+          double dz_tt=eta_jump*(ddv(0)*norm(d) + ddv(1)*tangent(d));
+          double dz_xx_jump=-dz_tt + dp_x_jump;
+
+          double dz_xx_correction=
+            -(z_jump - dx*dz_x_jump
+              + dx*dx*dz_xx_jump/2)/(h*h);
+          const int dz_xx_sign(x>middle ? -1 : 1);
+          Cx=dz_xx_sign*2*dz_xx_correction;
 
           double eta=exp(log_etax[i][j]);
           double delta=(h-std::fabs(dx))/h;
@@ -56,7 +72,7 @@ void compute_Cx(const Model &model, cons
           double dzx_xx_dv=-eta_jump*dvx_yy_dv + dp_x_dv;
           double dp_dv=-2*eta_jump*dvy_y_dv;
 
-          dCx_dvx=-dzx_xx_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
+          dCx_dvx=-dz_xx_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
                                   + dx*dx*dzx_xx_dv/2)/(h*h);
 
           if(x+h/2>middle && x-h/2<middle)
diff -r 114d616b8587 -r 7b5ea67a2e69 wscript
--- a/wscript	Mon Mar 05 17:40:26 2012 -0800
+++ b/wscript	Fri Mar 09 11:16:05 2012 -0800
@@ -1,6 +1,3 @@
-# top = '.'
-# out = 'build'
-
 def options(opt):
     opt.load('compiler_cxx')
 
@@ -18,6 +15,7 @@ def build(bld):
 
         target       = 'Gamr_disc',
         # cxxflags      = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG','-Wall'],
-        cxxflags      = ['-std=c++0x','-O3','-Wall'],
+        cxxflags      = ['-std=c++0x','-O3','-Wall','-Drestrict='],
         lib = ['boost_filesystem', 'boost_system', 'gsl', 'gslcblas'],
+        includes = ['/home/boo/programs/FTensor-bzr'],
         )



More information about the CIG-COMMITS mailing list