[cig-commits] commit: Reorganize most things into compute_coefficients.
Mercurial
hg at geodynamics.org
Fri Mar 30 09:32:35 PDT 2012
changeset: 142:7d5b5ce74b59
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 30 09:31:59 2012 -0700
files: Interface.hxx Interface/Interface.cxx compute_Cp.cxx compute_Cxyz.cxx compute_Cxyz.hxx compute_coefficients/Interface.hxx compute_coefficients/Interface/Interface.cxx compute_coefficients/compute_Cp.cxx compute_coefficients/compute_Cxyz.cxx compute_coefficients/compute_Cxyz.hxx compute_coefficients/compute_jumps.cxx compute_coefficients/compute_v_on_interface.hxx compute_coefficients/compute_v_on_interface/Valid.hxx compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx compute_coefficients/compute_v_on_interface/compute_values.hxx compute_coefficients/compute_v_on_interface/vel.hxx compute_coefficients/sign.hxx compute_coefficients/simplified_Rx.cxx compute_coefficients/simplified_Ry.cxx compute_jumps.cxx compute_v_on_interface.hxx compute_v_on_interface/Valid.hxx compute_v_on_interface/compute_1st_derivs.hxx compute_v_on_interface/compute_2nd_derivs.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx compute_v_on_interface/compute_values.hxx compute_v_on_interface/vel.hxx sign.hxx wscript
description:
Reorganize most things into compute_coefficients.
diff -r da0a37975701 -r 7d5b5ce74b59 Interface.hxx
--- a/Interface.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,86 +0,0 @@
-#ifndef GAMR_INTERFACE_HXX
-#define GAMR_INTERFACE_HXX
-
-#include "constants.hxx"
-#include "sign.hxx"
-#include "FTensor.hpp"
-#include <iostream>
-
-class Interface
-{
-public:
- FTensor::Tensor1<double,2> grid_pos, pos[2], corner_pos[2][2],
- anticorner_pos[2];
- int corner_dir;
- bool intersect_dd[2], intersect_sides[2], intersect_corner[2][2],
- intersect_anticorner[2][2];
- bool intersect_dp;
-
- Interface(const int &d, const int &i, const int &j,
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1]);
-
- bool intersects() const
- {
- return intersect_dd[0] || intersect_dd[1]
- || intersect_sides[0] || intersect_sides[1] || intersect_dp
- || intersect_corner[0][0] || intersect_corner[0][1]
- || intersect_corner[1][0] || intersect_corner[1][1];
- }
-
- template <int ny> void compute_xy(const int &i, const int &j,
- const double ¢er_dist,
- const double dist[][ny])
- {
- intersect_sides[0]=(sign(dist[i][j])!=sign(dist[i+1][j])
- && sign(dist[i][j+1])!=sign(dist[i+1][j+1])
- && sign(dist[i][j])==sign(dist[i][j+1]));
- intersect_sides[1]=(sign(dist[i][j])!=sign(dist[i][j+1])
- && sign(dist[i+1][j])!=sign(dist[i+1][j+1])
- && sign(dist[i][j])==sign(dist[i+1][j]));
- int center(sign(center_dist));
- for(int n0=0;n0<2;++n0)
- for(int n1=0;n1<2;++n1)
- {
- intersect_corner[n0][n1]=(center!=sign(dist[i+n0][j+n1]));
- if(intersect_corner[n0][n1])
- {
- double delta(center_dist-dist[i+n0][j+n1]);
- corner_pos[n0][n1](0)=
- grid_pos(0)+sign(n0-.5)*center_dist*h/(2*delta);
- corner_pos[n0][n1](1)=
- grid_pos(1)+sign(n1-.5)*center_dist*h/(2*delta);
- }
- }
- for(int n0=0;n0<2;++n0)
- for(int n1=0;n1<2;++n1)
- {
- intersect_anticorner[n0][n1]=
- !intersect_corner[n0][n1]
- && intersect_corner[(n0+1)%2][n1]
- && intersect_corner[n0][(n1+1)%2]
- && intersect_corner[(n0+1)%2][(n1+1)%2];
- if(intersect_anticorner[n0][n1])
- {
- /* Only need to compute positions for the diagonal
- points, since the center point has already been
- calculated for interface_corner. */
- double delta=dist[i+n0][j+n1]-dist[i+(n0+1)%2][j+n1];
- anticorner_pos[0](0)=grid_pos(0)
- +sign(n0-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
- anticorner_pos[0](1)=grid_pos(1)+sign(n1-.5)*h/2;
-
- delta=dist[i+n0][j+n1]-dist[i+n0][j+(n1+1)%2];
- anticorner_pos[1](0)=grid_pos(0)+sign(n0-.5)*h/2;
- anticorner_pos[1](1)=grid_pos(1)
- +sign(n1-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
-
- std::cerr << "Anticorners not tested yet"
- << std::endl;
- abort();
- }
- }
- }
-
-};
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 Interface/Interface.cxx
--- a/Interface/Interface.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,107 +0,0 @@
-#include "../Interface.hxx"
-
-Interface::Interface(const int &d, const int &i, const int &j,
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1])
-{
- switch(d)
- {
- case 0:
- grid_pos(0)=i*h;
- grid_pos(1)=(j+0.5)*h;
- break;
- case 1:
- grid_pos(0)=(i+0.5)*h;
- grid_pos(1)=j*h;
- break;
- default:
- abort();
- }
-
- if(d==0)
- {
- intersect_dd[0]=sign(distx[i][j])!=sign(distx[i+1][j])
- || sign(distx[i][j])!=sign(distx[i-1][j]);
- intersect_dd[1]=false;
- if(j!=Ny-1)
- {
- intersect_dd[1]=intersect_dd[1]
- || sign(distx[i][j])!=sign(distx[i][j+1]);
- }
- if(j!=0)
- {
- intersect_dd[1]=intersect_dd[1]
- || sign(distx[i][j-1])!=sign(distx[i][j]);
- }
-
- 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])
- {
- pos[0](1)=grid_pos(1);
- if(sign(distx[i][j])!=sign(distx[i+1][j]))
- delta=distx[i+1][j]-distx[i][j];
- else
- delta=distx[i][j]-distx[i-1][j];
-
- pos[0](0)=grid_pos(0)-distx[i][j]*h/delta;
- }
- if(intersect_dd[1])
- {
- pos[1](0)=grid_pos(0);
- if(j!=Ny-1 && sign(distx[i][j])!=sign(distx[i][j+1]))
- delta=distx[i][j+1]-distx[i][j];
- else
- delta=distx[i][j]-distx[i][j-1];
- pos[1](1)=grid_pos(1)-distx[i][j]*h/delta;
- }
-
- compute_xy(i-1,j,distx[i][j],disty);
- }
- else
- {
- intersect_dd[0]=false;
- if(i!=Nx-1)
- {
- intersect_dd[0]=intersect_dd[0]
- || sign(disty[i][j])!=sign(disty[i+1][j]);
- }
- if(i!=0)
- {
- intersect_dd[0]=intersect_dd[0]
- || sign(disty[i][j])!=sign(disty[i-1][j]);
- }
- intersect_dd[1]=sign(disty[i][j])!=sign(disty[i][j+1])
- || sign(disty[i][j-1])!=sign(disty[i][j]);
-
- intersect_dp=(sign(distx[i][j] + distx[i+1][j])
- !=sign(distx[i][j-1] + distx[i+1][j-1]));
-
-
- double delta;
- if(intersect_dd[0])
- {
- pos[0](1)=grid_pos(1);
- if(i!=Ny-1 && sign(disty[i][j])!=sign(disty[i+1][j]))
- delta=disty[i+1][j]-disty[i][j];
- else
- delta=disty[i][j]-disty[i-1][j];
-
- pos[0](0)=grid_pos(0)-disty[i][j]*h/delta;
- }
- if(intersect_dd[1])
- {
- pos[1](0)=grid_pos(0);
- if(sign(disty[i][j])!=sign(disty[i][j+1]))
- delta=disty[i][j+1]-disty[i][j];
- else
- delta=disty[i][j]-disty[i][j-1];
- pos[1](1)=grid_pos(1)-disty[i][j]*h/delta;
- }
-
- compute_xy(i,j-1,disty[i][j],distx);
- }
-}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_Cp.cxx
--- a/compute_Cp.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-#include "constants.hxx"
-#include <iostream>
-#include "compute_v_on_interface.hxx"
-#include "FTensor.hpp"
-#include "sign.hxx"
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump);
-
-void compute_Cp(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
- const int &i, const int &j,
- double &Cp)
-{
- Cp=0;
-
- bool intersects[]={sign(distx[i][j])!=sign(distx[i+1][j]),
- sign(disty[i][j])!=sign(disty[i][j+1])};
-
- 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])
- {
- interface_pos[0](0)=i*h-distx[i][j]*h/(distx[i+1][j]-distx[i][j]);
- interface_pos[0](1)=(j+0.5)*h;
- }
- if(intersects[1])
- {
- interface_pos[1](0)=(i+0.5)*h;
- interface_pos[1](1)=j*h-disty[i][j]*h/(disty[i][j+1]-disty[i][j]);
- }
-
- for(int dd=0;dd<2;++dd)
- if(intersects[dd])
- {
- FTensor::Tensor1<double,2> v, dv, ddv;
- FTensor::Tensor1<int,2> dir(0,0);
- FTensor::Tensor2<double,2,2> nt_to_xy;
- const int sgn(sign(pos(dd)-interface_pos[dd](dd)));
- dir(dd)=sgn;
- compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface_pos[dd],dd,v,dv,ddv,nt_to_xy);
-
- FTensor::Tensor1<double,2> z_jump;
- FTensor::Tensor2<double,2,2> dz_jump;
-
- double eta_jump(sgn);
- switch(dd)
- {
- case 0:
- eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i][j]);
- break;
- case 1:
- eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j]);
- break;
- default:
- abort();
- break;
- }
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
-
- double dx=std::fabs(pos(dd)-interface_pos[dd](dd) - (h/2)*sgn);
-
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
-
- Cp-=dir(a)*(z_jump(a) - dx*dz_jump(a,b)*dir(b))/h;
- }
-}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,193 +0,0 @@
-#include "constants.hxx"
-#include <iostream>
-#include "compute_v_on_interface.hxx"
-#include "FTensor.hpp"
-#include "Interface.hxx"
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump,
- FTensor::Tensor3_christof<double,2,2> &ddz_jump);
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump,
- FTensor::Tensor3_christof<double,2,2> &ddz_jump,
- double &p_jump,
- FTensor::Tensor1<double,2> &dp_jump);
-
-void compute_Cxyz(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
- const int &d,
- const int &i, const int &j,
- double &C)
-{
- C=0;
-
- const Interface interface(d,i,j,distx,disty);
- if(!interface.intersects())
- return;
- FTensor::Tensor1<double,2> pos;
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
- const FTensor::Index<'c',2> c;
- pos(a)=interface.grid_pos(a);
-
- FTensor::Tensor1<int,2> xyz(0,0);
- xyz(d)=1;
- bool try_corners(true);
- for(int dd=0;dd<2;++dd)
- {
- if(interface.intersect_dd[dd])
- {
- /* xyz= the component we are correcting, Cx, Cy, or Cz.
- dir=direction of 2nd derivative
- dir2, dir3=components of mixed derivative we are correcting.
- So for zy,xy with boundary in x, dir2=y, dir3=y. */
- 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,
- interface.pos[dd],d,
- v,dv,ddv,nt_to_xy);
-
- int sgn(sign(pos(dd) - interface.pos[dd](dd)));
- FTensor::Tensor1<double,2> dx(0,0);
- dx(dd)=pos(dd) - interface.pos[dd](dd) - h*sgn;
- FTensor::Tensor1<double,2> z_jump;
- FTensor::Tensor2<double,2,2> dz_jump;
- FTensor::Tensor3_christof<double,2,2> ddz_jump;
- double p_jump;
- FTensor::Tensor1<double,2> dp_jump;
-
- double eta_jump(sgn);
- switch(d)
- {
- case 0:
- switch(dd)
- {
- case 0:
- eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
- break;
- case 1:
- eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
- break;
- default:
- abort();
- break;
- }
- break;
- case 1:
- switch(dd)
- {
- case 0:
- eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
- break;
- case 1:
- eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
- break;
- default:
- abort();
- break;
- }
- break;
- default:
- abort();
- break;
- }
-
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
- p_jump,dp_jump);
-
- double dz_dd_correction=
- xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
- + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
-
- const int dz_dd_factor(d==dd ? 2 : 1);
- C+=dz_dd_factor*dz_dd_correction;
-
- if(d==dd && interface.intersect_dp)
- {
- dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
- C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
- }
-
- if(interface.intersect_sides[dd])
- {
- try_corners=false;
- FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
- dir2((d+1)%2)=1;
- dir3((dd+1)%2)=1;
- dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
-
- C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b)
- + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
- }
- }
- }
- for(int dd=0;dd<2;++dd)
- for(int ee=0;ee<2;++ee)
- if(interface.intersect_anticorner[dd][ee])
- {
- try_corners=false;
- std::cerr << "Anticorners not implemented"
- << std::endl;
- abort();
- }
- /* Handle corner cutting */
- if(try_corners)
- {
- FTensor::Tensor1<double,2> zyx(0,0);
- zyx((d+1)%2)=1;
-
- for(int dd=0;dd<2;++dd)
- for(int ee=0;ee<2;++ee)
- if(interface.intersect_corner[dd][ee])
- {
- 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,
- interface.corner_pos[dd][ee],d,
- v,dv,ddv,nt_to_xy);
-
- FTensor::Tensor1<double,2> z_jump;
- FTensor::Tensor2<double,2,2> dz_jump;
- FTensor::Tensor3_christof<double,2,2> ddz_jump;
-
- double eta_jump;
- switch(d)
- {
- case 0:
- eta_jump=exp(log_etax[i][j])-exp(log_etay[i-1+dd][j+ee]);
- break;
-
- case 1:
- eta_jump=exp(log_etay[i][j])-exp(log_etax[i+dd][j-1+ee]);
- break;
- default:
- abort();
- break;
- }
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
-
- FTensor::Tensor1<double,2> dx;
- for(int aa=0;aa<2;++aa)
- dx(aa)=pos(aa) - interface.corner_pos[dd][ee](aa)
- - (h/2)*sign(pos(aa) - interface.corner_pos[dd][ee](aa));
-
- C+=sign(dx(0))*sign(dx(1))*zyx(a)
- *(z_jump(a) + dz_jump(a,b)*dx(b)
- + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
- }
- }
-}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_Cxyz.hxx
--- a/compute_Cxyz.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,15 +0,0 @@
-#ifndef GAMR_COMPUTE_CXYZ_HXX
-#define GAMR_COMPUTE_CXYZ_HXX
-
-#include "constants.hxx"
-void compute_Cxyz(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
- const int &d,
- const int &i, const int &j,
- double &C);
-
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/Interface.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/Interface.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,86 @@
+#ifndef GAMR_INTERFACE_HXX
+#define GAMR_INTERFACE_HXX
+
+#include "../constants.hxx"
+#include "sign.hxx"
+#include "FTensor.hpp"
+#include <iostream>
+
+class Interface
+{
+public:
+ FTensor::Tensor1<double,2> grid_pos, pos[2], corner_pos[2][2],
+ anticorner_pos[2];
+ int corner_dir;
+ bool intersect_dd[2], intersect_sides[2], intersect_corner[2][2],
+ intersect_anticorner[2][2];
+ bool intersect_dp;
+
+ Interface(const int &d, const int &i, const int &j,
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1]);
+
+ bool intersects() const
+ {
+ return intersect_dd[0] || intersect_dd[1]
+ || intersect_sides[0] || intersect_sides[1] || intersect_dp
+ || intersect_corner[0][0] || intersect_corner[0][1]
+ || intersect_corner[1][0] || intersect_corner[1][1];
+ }
+
+ template <int ny> void compute_xy(const int &i, const int &j,
+ const double ¢er_dist,
+ const double dist[][ny])
+ {
+ intersect_sides[0]=(sign(dist[i][j])!=sign(dist[i+1][j])
+ && sign(dist[i][j+1])!=sign(dist[i+1][j+1])
+ && sign(dist[i][j])==sign(dist[i][j+1]));
+ intersect_sides[1]=(sign(dist[i][j])!=sign(dist[i][j+1])
+ && sign(dist[i+1][j])!=sign(dist[i+1][j+1])
+ && sign(dist[i][j])==sign(dist[i+1][j]));
+ int center(sign(center_dist));
+ for(int n0=0;n0<2;++n0)
+ for(int n1=0;n1<2;++n1)
+ {
+ intersect_corner[n0][n1]=(center!=sign(dist[i+n0][j+n1]));
+ if(intersect_corner[n0][n1])
+ {
+ double delta(center_dist-dist[i+n0][j+n1]);
+ corner_pos[n0][n1](0)=
+ grid_pos(0)+sign(n0-.5)*center_dist*h/(2*delta);
+ corner_pos[n0][n1](1)=
+ grid_pos(1)+sign(n1-.5)*center_dist*h/(2*delta);
+ }
+ }
+ for(int n0=0;n0<2;++n0)
+ for(int n1=0;n1<2;++n1)
+ {
+ intersect_anticorner[n0][n1]=
+ !intersect_corner[n0][n1]
+ && intersect_corner[(n0+1)%2][n1]
+ && intersect_corner[n0][(n1+1)%2]
+ && intersect_corner[(n0+1)%2][(n1+1)%2];
+ if(intersect_anticorner[n0][n1])
+ {
+ /* Only need to compute positions for the diagonal
+ points, since the center point has already been
+ calculated for interface_corner. */
+ double delta=dist[i+n0][j+n1]-dist[i+(n0+1)%2][j+n1];
+ anticorner_pos[0](0)=grid_pos(0)
+ +sign(n0-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
+ anticorner_pos[0](1)=grid_pos(1)+sign(n1-.5)*h/2;
+
+ delta=dist[i+n0][j+n1]-dist[i+n0][j+(n1+1)%2];
+ anticorner_pos[1](0)=grid_pos(0)+sign(n0-.5)*h/2;
+ anticorner_pos[1](1)=grid_pos(1)
+ +sign(n1-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
+
+ std::cerr << "Anticorners not tested yet"
+ << std::endl;
+ abort();
+ }
+ }
+ }
+
+};
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/Interface/Interface.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/Interface/Interface.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,107 @@
+#include "../Interface.hxx"
+
+Interface::Interface(const int &d, const int &i, const int &j,
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1])
+{
+ switch(d)
+ {
+ case 0:
+ grid_pos(0)=i*h;
+ grid_pos(1)=(j+0.5)*h;
+ break;
+ case 1:
+ grid_pos(0)=(i+0.5)*h;
+ grid_pos(1)=j*h;
+ break;
+ default:
+ abort();
+ }
+
+ if(d==0)
+ {
+ intersect_dd[0]=sign(distx[i][j])!=sign(distx[i+1][j])
+ || sign(distx[i][j])!=sign(distx[i-1][j]);
+ intersect_dd[1]=false;
+ if(j!=Ny-1)
+ {
+ intersect_dd[1]=intersect_dd[1]
+ || sign(distx[i][j])!=sign(distx[i][j+1]);
+ }
+ if(j!=0)
+ {
+ intersect_dd[1]=intersect_dd[1]
+ || sign(distx[i][j-1])!=sign(distx[i][j]);
+ }
+
+ 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])
+ {
+ pos[0](1)=grid_pos(1);
+ if(sign(distx[i][j])!=sign(distx[i+1][j]))
+ delta=distx[i+1][j]-distx[i][j];
+ else
+ delta=distx[i][j]-distx[i-1][j];
+
+ pos[0](0)=grid_pos(0)-distx[i][j]*h/delta;
+ }
+ if(intersect_dd[1])
+ {
+ pos[1](0)=grid_pos(0);
+ if(j!=Ny-1 && sign(distx[i][j])!=sign(distx[i][j+1]))
+ delta=distx[i][j+1]-distx[i][j];
+ else
+ delta=distx[i][j]-distx[i][j-1];
+ pos[1](1)=grid_pos(1)-distx[i][j]*h/delta;
+ }
+
+ compute_xy(i-1,j,distx[i][j],disty);
+ }
+ else
+ {
+ intersect_dd[0]=false;
+ if(i!=Nx-1)
+ {
+ intersect_dd[0]=intersect_dd[0]
+ || sign(disty[i][j])!=sign(disty[i+1][j]);
+ }
+ if(i!=0)
+ {
+ intersect_dd[0]=intersect_dd[0]
+ || sign(disty[i][j])!=sign(disty[i-1][j]);
+ }
+ intersect_dd[1]=sign(disty[i][j])!=sign(disty[i][j+1])
+ || sign(disty[i][j-1])!=sign(disty[i][j]);
+
+ intersect_dp=(sign(distx[i][j] + distx[i+1][j])
+ !=sign(distx[i][j-1] + distx[i+1][j-1]));
+
+
+ double delta;
+ if(intersect_dd[0])
+ {
+ pos[0](1)=grid_pos(1);
+ if(i!=Ny-1 && sign(disty[i][j])!=sign(disty[i+1][j]))
+ delta=disty[i+1][j]-disty[i][j];
+ else
+ delta=disty[i][j]-disty[i-1][j];
+
+ pos[0](0)=grid_pos(0)-disty[i][j]*h/delta;
+ }
+ if(intersect_dd[1])
+ {
+ pos[1](0)=grid_pos(0);
+ if(sign(disty[i][j])!=sign(disty[i][j+1]))
+ delta=disty[i][j+1]-disty[i][j];
+ else
+ delta=disty[i][j]-disty[i][j-1];
+ pos[1](1)=grid_pos(1)-disty[i][j]*h/delta;
+ }
+
+ compute_xy(i,j-1,disty[i][j],distx);
+ }
+}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_Cp.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_Cp.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,79 @@
+#include "../constants.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+#include "FTensor.hpp"
+#include "sign.hxx"
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump);
+
+void compute_Cp(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cp)
+{
+ Cp=0;
+
+ bool intersects[]={sign(distx[i][j])!=sign(distx[i+1][j]),
+ sign(disty[i][j])!=sign(disty[i][j+1])};
+
+ 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])
+ {
+ interface_pos[0](0)=i*h-distx[i][j]*h/(distx[i+1][j]-distx[i][j]);
+ interface_pos[0](1)=(j+0.5)*h;
+ }
+ if(intersects[1])
+ {
+ interface_pos[1](0)=(i+0.5)*h;
+ interface_pos[1](1)=j*h-disty[i][j]*h/(disty[i][j+1]-disty[i][j]);
+ }
+
+ for(int dd=0;dd<2;++dd)
+ if(intersects[dd])
+ {
+ FTensor::Tensor1<double,2> v, dv, ddv;
+ FTensor::Tensor1<int,2> dir(0,0);
+ FTensor::Tensor2<double,2,2> nt_to_xy;
+ const int sgn(sign(pos(dd)-interface_pos[dd](dd)));
+ dir(dd)=sgn;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface_pos[dd],dd,v,dv,ddv,nt_to_xy);
+
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
+
+ double eta_jump(sgn);
+ switch(dd)
+ {
+ case 0:
+ eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i][j]);
+ break;
+ case 1:
+ eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
+
+ double dx=std::fabs(pos(dd)-interface_pos[dd](dd) - (h/2)*sgn);
+
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+
+ Cp-=dir(a)*(z_jump(a) - dx*dz_jump(a,b)*dir(b))/h;
+ }
+}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_Cxyz.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_Cxyz.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,193 @@
+#include "../constants.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+#include "FTensor.hpp"
+#include "Interface.hxx"
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump);
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump,
+ double &p_jump,
+ FTensor::Tensor1<double,2> &dp_jump);
+
+void compute_Cxyz(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const int &d,
+ const int &i, const int &j,
+ double &C)
+{
+ C=0;
+
+ const Interface interface(d,i,j,distx,disty);
+ if(!interface.intersects())
+ return;
+ FTensor::Tensor1<double,2> pos;
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ pos(a)=interface.grid_pos(a);
+
+ FTensor::Tensor1<int,2> xyz(0,0);
+ xyz(d)=1;
+ bool try_corners(true);
+ for(int dd=0;dd<2;++dd)
+ {
+ if(interface.intersect_dd[dd])
+ {
+ /* xyz= the component we are correcting, Cx, Cy, or Cz.
+ dir=direction of 2nd derivative
+ dir2, dir3=components of mixed derivative we are correcting.
+ So for zy,xy with boundary in x, dir2=y, dir3=y. */
+ 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,
+ interface.pos[dd],d,
+ v,dv,ddv,nt_to_xy);
+
+ int sgn(sign(pos(dd) - interface.pos[dd](dd)));
+ FTensor::Tensor1<double,2> dx(0,0);
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - h*sgn;
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ double p_jump;
+ FTensor::Tensor1<double,2> dp_jump;
+
+ double eta_jump(sgn);
+ switch(d)
+ {
+ case 0:
+ switch(dd)
+ {
+ case 0:
+ eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
+ break;
+ case 1:
+ eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ break;
+ case 1:
+ switch(dd)
+ {
+ case 0:
+ eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
+ break;
+ case 1:
+ eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ break;
+ default:
+ abort();
+ break;
+ }
+
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
+ p_jump,dp_jump);
+
+ double dz_dd_correction=
+ xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
+ + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+
+ const int dz_dd_factor(d==dd ? 2 : 1);
+ C+=dz_dd_factor*dz_dd_correction;
+
+ if(d==dd && interface.intersect_dp)
+ {
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
+ C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
+ }
+
+ if(interface.intersect_sides[dd])
+ {
+ try_corners=false;
+ FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
+ dir2((d+1)%2)=1;
+ dir3((dd+1)%2)=1;
+ dx(dd)=pos(dd) - interface.pos[dd](dd) - (h/2)*sgn;
+
+ C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b)
+ + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
+ }
+ }
+ }
+ for(int dd=0;dd<2;++dd)
+ for(int ee=0;ee<2;++ee)
+ if(interface.intersect_anticorner[dd][ee])
+ {
+ try_corners=false;
+ std::cerr << "Anticorners not implemented"
+ << std::endl;
+ abort();
+ }
+ /* Handle corner cutting */
+ if(try_corners)
+ {
+ FTensor::Tensor1<double,2> zyx(0,0);
+ zyx((d+1)%2)=1;
+
+ for(int dd=0;dd<2;++dd)
+ for(int ee=0;ee<2;++ee)
+ if(interface.intersect_corner[dd][ee])
+ {
+ 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,
+ interface.corner_pos[dd][ee],d,
+ v,dv,ddv,nt_to_xy);
+
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+
+ double eta_jump;
+ switch(d)
+ {
+ case 0:
+ eta_jump=exp(log_etax[i][j])-exp(log_etay[i-1+dd][j+ee]);
+ break;
+
+ case 1:
+ eta_jump=exp(log_etay[i][j])-exp(log_etax[i+dd][j-1+ee]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
+
+ FTensor::Tensor1<double,2> dx;
+ for(int aa=0;aa<2;++aa)
+ dx(aa)=pos(aa) - interface.corner_pos[dd][ee](aa)
+ - (h/2)*sign(pos(aa) - interface.corner_pos[dd][ee](aa));
+
+ C+=sign(dx(0))*sign(dx(1))*zyx(a)
+ *(z_jump(a) + dz_jump(a,b)*dx(b)
+ + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+ }
+ }
+}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_Cxyz.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_Cxyz.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,15 @@
+#ifndef GAMR_COMPUTE_CXYZ_HXX
+#define GAMR_COMPUTE_CXYZ_HXX
+
+#include "../constants.hxx"
+void compute_Cxyz(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const int &d,
+ const int &i, const int &j,
+ double &C);
+
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_jumps.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_jumps.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,81 @@
+#include "FTensor.hpp"
+#include "../constants.hxx"
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump,
+ double &p_jump,
+ FTensor::Tensor1<double,2> &dp_jump)
+{
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ const FTensor::Index<'d',2> d;
+ const FTensor::Index<'e',2> e;
+ const FTensor::Index<'f',2> f;
+
+ z_jump(a)=eta_jump*v(b)*nt_to_xy(a,b);
+
+ FTensor::Number<0> n;
+ FTensor::Number<1> t;
+ FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+
+ FTensor::Tensor2<double,2,2> dz_jump_nt;
+ dz_jump_nt(a,n)=-eta_jump*dv(b)*ddel(a,b);
+ dz_jump_nt(a,t)=eta_jump*dv(a);
+ dz_jump(a,b)=dz_jump_nt(c,d)*nt_to_xy(a,c)*nt_to_xy(b,d);
+
+ p_jump=-2*eta_jump*dv(t);
+
+ /* TODO: dp_jump and ddz_jump both need to be corrected for
+ curvature and density jumps. */
+
+ FTensor::Tensor1<double,2> dp_jump_nt;
+ dp_jump_nt(n)=2*eta_jump*ddv(n);
+ dp_jump_nt(t)=-2*eta_jump*ddv(t);
+ dp_jump(a)=nt_to_xy(a,b)*dp_jump_nt(b);
+
+ FTensor::Tensor3_christof<double,2,2> ddz_jump_nt;
+ ddz_jump_nt(a,t,t)=eta_jump*ddv(a);
+ ddz_jump_nt(a,n,n)=-ddz_jump_nt(a,t,t) + dp_jump_nt(a);
+ ddz_jump_nt(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
+
+ /* It would be nice if there were a way to do this directly instead
+ of iterating over all of the indices. */
+ ddz_jump(a,n,n)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(n,f);
+ ddz_jump(a,n,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(t,f);
+ ddz_jump(a,t,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(t,e)*nt_to_xy(t,f);
+}
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump)
+{
+ double p_jump;
+ FTensor::Tensor1<double,2> dp_jump;
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
+ p_jump,dp_jump);
+}
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump)
+{
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
+}
+
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,19 @@
+#ifndef GAMR_COMPUTE_V_ON_INTERFACE_HXX
+#define GAMR_COMPUTE_V_ON_INTERFACE_HXX
+
+#include "FTensor.hpp"
+
+void compute_v_on_interface(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const FTensor::Tensor1<double,2> &pos,
+ const int &xyz,
+ FTensor::Tensor1<double,2> &v,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor2<double,2,2> &nt_to_xy);
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/Valid.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/Valid.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,30 @@
+#ifndef GAMR_VALID_HXX
+#define GAMR_VALID_HXX
+
+#include <limits>
+
+class Valid
+{
+public:
+ int sign,i,j;
+ bool valid;
+ FTensor::Tensor1<double,2> diff;
+ double r;
+
+ Valid(): sign(0), valid(false), r(std::numeric_limits<double>::infinity()) {}
+ Valid(const int &Sign, const bool &val, const int &I, const int &J,
+ const FTensor::Tensor1<double,2> &Diff):
+ sign(Sign), i(I), j(J), valid(val), diff(Diff),
+ r(diff(0)*diff(0)+diff(1)*diff(1)) {}
+ bool operator<(const Valid &v) const
+ {
+ return r<v.r;
+ }
+ bool operator>(const Valid &v) const
+ {
+ return r>v.r;
+ }
+};
+
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,66 @@
+#ifndef GAMR_COMPUTE_1ST_DERIVS_HXX
+#define GAMR_COMPUTE_1ST_DERIVS_HXX
+
+#include <cassert>
+#include "Valid.hxx"
+#include <gsl/gsl_linalg.h>
+
+template<int ny>
+void compute_1st_derivs(const int &d, const std::vector<Valid> &derivs,
+ const int &d0,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ FTensor::Tensor2<double,2,2> &dv_pm)
+{
+ /* Create and invert the Vandermonde matrix */
+ assert(derivs.size()==3);
+ double m[3*3], rhs[3];
+ for(int i=0;i<3;++i)
+ {
+ switch(d0)
+ {
+ case 0:
+ if(derivs[i].i==nx-1)
+ {
+ rhs[i]=0;
+ }
+ else
+ {
+ rhs[i]=vel(z,log_eta,derivs[i].i+1,derivs[i].j)
+ - vel(z,log_eta,derivs[i].i,derivs[i].j);
+ }
+ break;
+ case 1:
+ if(derivs[i].j==ny-1)
+ {
+ rhs[i]=0;
+ }
+ else
+ {
+ rhs[i]=vel(z,log_eta,derivs[i].i,derivs[i].j+1)
+ - vel(z,log_eta,derivs[i].i,derivs[i].j);
+ }
+ break;
+ default:
+ abort();
+ }
+ m[0+3*i]=1;
+ m[1+3*i]=derivs[i].diff(0);
+ m[2+3*i]=derivs[i].diff(1);
+ }
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,3,3);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,3);
+ gsl_vector *x=gsl_vector_alloc(3);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(3);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
+
+ dv_pm(d,d0)=gsl_vector_get(x,0);
+ gsl_vector_free(x);
+ gsl_permutation_free(perm);
+}
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,74 @@
+#ifndef GAMR_COMPUTE_2ND_DERIVS_HXX
+#define GAMR_COMPUTE_2ND_DERIVS_HXX
+
+#include "Valid.hxx"
+
+template<int ny>
+void compute_2nd_derivs(const int &d, const Valid &v,
+ const int &d0, const int &d1,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ FTensor::Tensor3_christof<double,2,2> &ddv_pm,
+ double &eta_pm)
+{
+ if(d0!=d1)
+ {
+ if(v.i==nx-1 || v.j==ny-1)
+ {
+ ddv_pm(d,d0,d1)=0;
+ }
+ else
+ {
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j+1)
+ - vel(z,log_eta,v.i,v.j+1)
+ - vel(z,log_eta,v.i+1,v.j)
+ + vel(z,log_eta,v.i,v.j);
+ }
+ }
+ else
+ {
+ switch(d0)
+ {
+ case 0:
+ if(v.i==0)
+ {
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ - vel(z,log_eta,v.i,v.j);
+ }
+ else if(v.i==nx-1)
+ {
+ ddv_pm(d,d0,d1)=
+ - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
+ }
+ else
+ {
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
+ }
+ break;
+ case 1:
+ if(v.j==0)
+ {
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
+ - vel(z,log_eta,v.i,v.j);
+ }
+ else if(v.j==ny-1)
+ {
+ ddv_pm(d,d0,d1)=
+ - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
+ }
+ else
+ {
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
+ - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
+ }
+ break;
+ default:
+ abort();
+ }
+ }
+ eta_pm=exp(log_eta[v.i][v.j]);
+}
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,413 @@
+#include "../constants.hxx"
+#include <cmath>
+#include "FTensor.hpp"
+#include <vector>
+#include "../sign.hxx"
+#include "vel.hxx"
+#include "compute_2nd_derivs.hxx"
+#include "compute_1st_derivs.hxx"
+#include "compute_values.hxx"
+#include <tuple>
+#include <algorithm>
+#include <queue>
+
+void compute_dv_dtt(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const FTensor::Tensor1<double,2> &pos,
+ const FTensor::Tensor1<double,2> &norm,
+ const FTensor::Tensor1<double,2> &tangent,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &v)
+{
+ double length=h/std::max(std::abs(norm(0)),std::abs(norm(1)));
+
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ FTensor::Tensor3_christof<double,2,2> ddv_pm[2];
+ ddv_pm[0](a,b,c)=0;
+ ddv_pm[1](a,b,c)=0;
+
+ FTensor::Tensor2<double,2,2> dv_pm[2];
+ dv_pm[0](a,b)=0;
+ dv_pm[1](a,b)=0;
+
+ FTensor::Tensor1<double,2> v_pm[2];
+ v_pm[0](a)=0;
+ v_pm[1](a)=0;
+
+ double eta_pm[2];
+
+ /* 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));
+
+ std::vector<Valid> valid((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);
+
+ int starti((pos(0)-dx)/h), startj((pos(1)-dy)/h);
+
+ /* d/dx^2 */
+ 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), p_d(i*h+dx+h/2,j*h+dy),
+ dd_diff, d_diff, diff;
+ dd_diff(a)=p(a)-pos(a);
+ d_diff(a)=p_d(a)-pos(a);
+
+ if(d==0)
+ {
+ /* ddv */
+ if(i>0 && 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)
+ {
+ d_diff(a)-=length*norm(a)*sign(distx[i][j]);
+ 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);
+ }
+ /* v */
+ diff(a)=dd_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
+ {
+ /* ddv */
+ 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
+ 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_diff(a)-=length*norm(a)*sign(disty[i][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);
+
+ /* v */
+ diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(disty[i][j]);
+ valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),true,i,j,diff);
+ }
+ }
+
+ /* d/dy^2 */
+ 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), 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)
+ {
+ /* ddv */
+ bool v;
+ if(j==0)
+ v=sign(distx[i][j])==sign(distx[i][j+1]);
+ else if(j==max_y)
+ v=sign(distx[i][j])==sign(distx[i][j-1]);
+ else
+ v=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,i,j,dd_diff);
+
+ /* dv */
+ bool v_d;
+ if(j==max_y)
+ v_d=true;
+ else
+ v_d=sign(distx[i][j])==sign(distx[i][j+1]);
+
+ d_diff(a)-=length*norm(a)*sign(distx[i][j]);
+ 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);
+ }
+ else
+ {
+ /* ddv */
+ 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);
+
+ /* dv */
+ if(j<max_y)
+ {
+ d_diff(a)-=length*norm(a)*sign(disty[i][j]);
+ 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);
+ }
+ }
+ }
+
+ /* 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_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<max_x)
+ {
+ bool v;
+ if(j==max_y)
+ v=sign(distx[i][j])==sign(distx[i+1][j]);
+ else
+ v=sign(distx[i][j])==sign(distx[i][j+1])
+ && sign(distx[i][j])==sign(distx[i+1][j])
+ && sign(distx[i][j])==sign(distx[i+1][j+1]);
+
+ 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);
+ }
+ }
+ else
+ {
+ if(j<max_y)
+ {
+ bool v;
+ if(i==max_x)
+ v=sign(distx[i][j])==sign(distx[i][j+1]);
+ else
+ v=sign(distx[i][j])==sign(distx[i][j+1])
+ && sign(distx[i][j])==sign(distx[i+1][j])
+ && sign(distx[i][j])==sign(distx[i+1][j+1]);
+
+ 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);
+ }
+ }
+ }
+
+ /* Put the valid vectors into a priority queue so that we can
+ take the sorted versions out as needed. A bit messier than
+ just sorting the vector, but also significantly faster. */
+
+ FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
+ false,false);
+
+ typedef std::priority_queue<Valid, std::vector<Valid>, std::greater<Valid> > mpq;
+ mpq ddp_valid[]={{dd_valid(0,0).begin(),dd_valid(0,0).end()},
+ {dd_valid(0,1).begin(),dd_valid(0,1).end()},
+ {dd_valid(1,1).begin(),dd_valid(1,1).end()}};
+
+ mpq dp_valid[]={{d_valid(0).begin(),d_valid(0).end()},
+ {d_valid(1).begin(),d_valid(1).end()}};
+ std::priority_queue<Valid, std::vector<Valid>, std::greater<Valid> >
+ p_valid(valid.begin(),valid.end());
+
+ /* Compute everything */
+ /* v */
+
+ std::vector<Valid> values[2];
+ for(auto V=p_valid.top(); !p_valid.empty(); p_valid.pop(), V=p_valid.top())
+ {
+ if(values[0].size()==6 && values[1].size()==6)
+ break;
+
+ if(V.sign==0)
+ continue;
+ const int pm=(V.sign==-1 ? 0 : 1);
+ switch(values[pm].size())
+ {
+ case 0:
+ case 1:
+ case 2:
+ values[pm].push_back(V);
+ break;
+ case 6:
+ break;
+ default:
+ int x(0),y(0);
+ std::vector<int> x_coords={V.i};
+ std::vector<int> y_coords={V.j};
+ for(auto &value: values[pm])
+ {
+ if(V.i==value.i)
+ ++x;
+ if(V.j==value.j)
+ ++y;
+ if(values[pm].size()==5)
+ {
+ if(find(x_coords.begin(),x_coords.end(),value.i)==x_coords.end())
+ x_coords.push_back(value.i);
+ if(find(y_coords.begin(),y_coords.end(),value.j)==y_coords.end())
+ y_coords.push_back(value.j);
+ }
+ }
+ if(x<3 && y<3 && !(values[pm].size()==5 &&
+ (x_coords.size()<3 || y_coords.size()<3)))
+ values[pm].push_back(V);
+ 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]);
+ else
+ compute_values(1,values[pm],Nx,zy,log_etay,
+ norm,pm==0 ? -length : length,v_pm[pm]);
+ }
+ break;
+ }
+ }
+
+ for(int d0=0;d0<2;++d0)
+ {
+ /* First derivative dv */
+ std::vector<Valid> derivs[2];
+ for(auto V=dp_valid[d0].top(); !dp_valid[d0].empty();
+ dp_valid[d0].pop(), V=dp_valid[d0].top())
+ {
+ if(derivs[0].size()==3 && derivs[1].size()==3)
+ break;
+
+ if(V.sign==0)
+ continue;
+ if(V.valid)
+ {
+ const int pm=(V.sign==-1 ? 0 : 1);
+ switch(derivs[pm].size())
+ {
+ default:
+ break;
+ case 0:
+ case 1:
+ derivs[pm].push_back(V);
+ break;
+ case 2:
+ Valid &v0(derivs[pm][0]);
+ Valid &v1(derivs[pm][1]);
+ if(!((v0.i==v1.i && v1.i==V.i)
+ || (v0.j==v1.j && v1.j==V.j)))
+ {
+ derivs[pm].push_back(V);
+ if(d==0)
+ compute_1st_derivs(0,derivs[pm],d0,Nx+1,
+ zx,log_etax,dv_pm[pm]);
+ else
+ compute_1st_derivs(1,derivs[pm],d0,Nx,
+ zy,log_etay,dv_pm[pm]);
+ }
+ break;
+ }
+ }
+ }
+ /* Second derivative ddv */
+ for(int d1=d0;d1<2;++d1)
+ {
+ for(auto V=ddp_valid[2*d0+d1].top(); !ddp_valid[2*d0+d1].empty();
+ ddp_valid[2*d0+d1].pop(), V=ddp_valid[2*d0+d1].top())
+ {
+ if(is_set(0,d0,d1) && is_set(1,d0,d1))
+ break;
+
+ for(int pm=0;pm<2;++pm)
+ {
+ if(V.sign==(pm==0 ? -1 : 1) && V.valid && !is_set(pm,d0,d1))
+ {
+ is_set(pm,d0,d1)=true;
+ if(d==0)
+ compute_2nd_derivs(0,V,d0,d1,Nx+1,zx,log_etax,
+ ddv_pm[pm],eta_pm[pm]);
+ else
+ compute_2nd_derivs(1,V,d0,d1,Nx,zy,log_etay,
+ ddv_pm[pm],eta_pm[pm]);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* ddv */
+ double temp(0);
+ FTensor::Tensor1<double,2> ddv_xy(0,0);
+ for(int d=0;d<2;++d)
+ {
+ temp+=eta_pm[d]*eta_pm[d];
+ ddv_xy(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
+ }
+ ddv_xy(a)/=temp*h*h;
+
+ ddv(a)=nt_to_xy(b,a)*ddv_xy(b);
+
+ /* dv */
+ FTensor::Tensor1<double,2> dv_xy(0,0);
+ for(int d=0;d<2;++d)
+ dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
+
+ dv(a)=nt_to_xy(b,a)*dv_xy(b);
+
+ FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+ FTensor::Tensor1<double,2> ddz_nt_jump;
+ double eta_jump(eta_pm[1]-eta_pm[0]);
+ /* TODO: fix for curved surfaces */
+ ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
+
+ dv(a)-=h*length*ddz_nt_jump(a);
+ dv(a)/=(eta_pm[0]+eta_pm[1])*h;
+
+ /* v */
+ FTensor::Tensor1<double,2> v_xy(0,0);
+ for(int d=0;d<2;++d)
+ v_xy(a)+=v_pm[d](a)*eta_pm[d];
+
+ v(a)=nt_to_xy(b,a)*v_xy(b);
+
+ /* TODO: fix for curved surfaces */
+ FTensor::Tensor1<double,2> dz_n_jump;
+ dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
+
+ v(a)-=2*length*dz_n_jump(a)/3;
+ v(a)/=(eta_pm[0] + eta_pm[1]);
+}
+
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,66 @@
+#include <algorithm>
+#include <cmath>
+#include "../constants.hxx"
+#include <limits>
+#include <iostream>
+#include "FTensor.hpp"
+#include "../compute_v_on_interface.hxx"
+#include "vel.hxx"
+
+void compute_dv_dtt(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const FTensor::Tensor1<double,2> &pos,
+ const FTensor::Tensor1<double,2> &norm,
+ const FTensor::Tensor1<double,2> &tangent,
+ const FTensor::Tensor2<double,2,2> &nt_to_xy,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &v);
+
+void compute_v_on_interface(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const FTensor::Tensor1<double,2> &pos,
+ const int &xyz,
+ FTensor::Tensor1<double,2> &v,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &ddv,
+ 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));
+
+ tangent(0)=-norm(1);
+ tangent(1)=norm(0);
+
+ FTensor::Number<0> n;
+ FTensor::Number<1> t;
+ nt_to_xy(a,n)=norm(a);
+ nt_to_xy(a,t)=tangent(a);
+
+ compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,nt_to_xy,
+ ddv,dv,v);
+}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/compute_values.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/compute_values.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,61 @@
+#ifndef GAMR_COMPUTE_VALUES_HXX
+#define GAMR_COMPUTE_VALUES_HXX
+
+#include <cassert>
+#include "Valid.hxx"
+#include <gsl/gsl_linalg.h>
+
+template<int ny>
+void compute_values(const int &d, const std::vector<Valid> &values,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ const FTensor::Tensor1<double,2> &norm,
+ const double &signed_length,
+ FTensor::Tensor1<double,2> &v_pm)
+{
+ /* Create and invert the Vandermonde matrix */
+ assert(values.size()==6);
+ double m[6*6], rhs[6];
+ for(int i=0;i<6;++i)
+ {
+ rhs[i]=vel(z,log_eta,values[i].i,values[i].j);
+
+ m[0+6*i]=1;
+ m[1+6*i]=values[i].diff(0);
+ m[2+6*i]=values[i].diff(1);
+ m[3+6*i]=values[i].diff(0)*values[i].diff(0)/2;
+ m[4+6*i]=values[i].diff(1)*values[i].diff(1)/2;
+ m[5+6*i]=values[i].diff(0)*values[i].diff(1);
+ }
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,6,6);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,6);
+ gsl_vector *x=gsl_vector_alloc(6);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(6);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
+
+ double v=gsl_vector_get(x,0);
+
+ FTensor::Tensor1<double,2> dv_xy;
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ dv_xy(0)=gsl_vector_get(x,1);
+ dv_xy(1)=gsl_vector_get(x,2);
+ double dv=dv_xy(a)*norm(a)*signed_length/2;
+
+ FTensor::Tensor2_symmetric<double,2> ddv_xy;
+ ddv_xy(0,0)=gsl_vector_get(x,3);
+ ddv_xy(1,1)=gsl_vector_get(x,4);
+ ddv_xy(0,1)=gsl_vector_get(x,5);
+ double ddv=ddv_xy(a,b)*norm(a)*norm(b)*signed_length*signed_length/4;
+
+ v_pm(d)=(3*v - 5*dv + 1.5*ddv)/3;
+
+ gsl_vector_free(x);
+ gsl_permutation_free(perm);
+}
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/compute_v_on_interface/vel.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_v_on_interface/vel.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,75 @@
+#ifndef GAMR_VEL_HXX
+#define GAMR_VEL_HXX
+
+inline double analytic(const double x, const double y, const bool return_ux)
+{
+ const double eta_jump=max_eta-min_eta;
+
+ double ux(1.1), ux_xp(1.2), ux_xxp(1.3), ux_xyp(1.4),
+ ux_xm(1.5), ux_xxm(1.6), ux_xym(1.7),
+ uy(2.1), uy_xp(2.2), uy_xxp(2.3), uy_xyp(2.4),
+ uy_xm(2.5), uy_xxm(2.6), uy_xym(2.7);
+ double ux_y(-(max_eta*uy_xp - min_eta*uy_xm)/eta_jump),
+ ux_yy(-(max_eta*uy_xyp - min_eta*uy_xym)/eta_jump),
+ uy_y(-(max_eta*ux_xp - min_eta*ux_xm)/eta_jump),
+ uy_yy(-(max_eta*ux_xyp - min_eta*ux_xym)/eta_jump);
+
+ // std::cout << "u y "
+ // << ux << " "
+ // << ux_y << " "
+ // << ux_yy << " "
+ // << uy << " "
+ // << uy_y << " "
+ // << uy_yy << " "
+ // << "\n";
+ // exit(0);
+
+ if(return_ux)
+ {
+ if(x>0)
+ return ux + x*ux_xp + x*x*ux_xxp/2 + x*y*ux_xyp
+ + y*ux_y + y*y*ux_yy/2;
+ else
+ return ux + x*ux_xm + x*x*ux_xxm/2 + x*y*ux_xym
+ + y*ux_y + y*y*ux_yy/2;
+ }
+ else
+ {
+ if(x>0)
+ return uy + x*uy_xp + x*x*uy_xxp/2 + x*y*uy_xyp
+ + y*uy_y + y*y*uy_yy/2;
+ else
+ return uy + x*uy_xm + x*x*uy_xxm/2 + x*y*uy_xym
+ + y*uy_y + y*y*uy_yy/2;
+ }
+
+}
+
+
+template<int ny>
+double vel(const double z[][ny], const double log_eta[][ny],
+ const int i, const int j)
+{
+ // if(j>ny/8-4 && j<ny/8+4)
+ // {
+ // double x,y;
+ // if(ny%2==0)
+ // {
+ // x=i*h-middle;
+ // y=(j+0.5)*h-1.0/8;
+ // // y=j*h-1.0/8;
+ // return analytic(x,y,true);
+ // }
+ // else
+ // {
+ // x=(i+0.5)*h-middle;
+ // y=j*h-1.0/8;
+ // // y=(j-0.5)*h-1.0/8;
+ // return analytic(x,y,false);
+ // }
+ // }
+
+ return z[i][j]*std::exp(-log_eta[i][j]);
+}
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/sign.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/sign.hxx Fri Mar 30 09:31:59 2012 -0700
@@ -0,0 +1,8 @@
+#ifndef GAMR_SIGN_HXX
+#define GAMR_SIGN_HXX
+
+template <typename T> int sign(T val) {
+ return (val > T(0)) - (val < T(0));
+}
+
+#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/simplified_Rx.cxx
--- a/compute_coefficients/simplified_Rx.cxx Fri Mar 30 09:31:23 2012 -0700
+++ b/compute_coefficients/simplified_Rx.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -1,7 +1,7 @@
/* Compute only the parts of Rx that depend on zx, zy */
#include "../constants.hxx"
-#include "../compute_Cxyz.hxx"
+#include "compute_Cxyz.hxx"
double simplified_Rx(const int &i, const int &j,
const double zx[Nx+1][Ny],
diff -r da0a37975701 -r 7d5b5ce74b59 compute_coefficients/simplified_Ry.cxx
--- a/compute_coefficients/simplified_Ry.cxx Fri Mar 30 09:31:23 2012 -0700
+++ b/compute_coefficients/simplified_Ry.cxx Fri Mar 30 09:31:59 2012 -0700
@@ -1,7 +1,7 @@
/* Compute only the parts of Ry that depend on zx, zy */
#include "../constants.hxx"
-#include "../compute_Cxyz.hxx"
+#include "compute_Cxyz.hxx"
double simplified_Ry(const int &i, const int &j,
const double zx[Nx+1][Ny],
diff -r da0a37975701 -r 7d5b5ce74b59 compute_jumps.cxx
--- a/compute_jumps.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-#include "FTensor.hpp"
-#include "constants.hxx"
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump,
- FTensor::Tensor3_christof<double,2,2> &ddz_jump,
- double &p_jump,
- FTensor::Tensor1<double,2> &dp_jump)
-{
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
- const FTensor::Index<'c',2> c;
- const FTensor::Index<'d',2> d;
- const FTensor::Index<'e',2> e;
- const FTensor::Index<'f',2> f;
-
- z_jump(a)=eta_jump*v(b)*nt_to_xy(a,b);
-
- FTensor::Number<0> n;
- FTensor::Number<1> t;
- FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
-
- FTensor::Tensor2<double,2,2> dz_jump_nt;
- dz_jump_nt(a,n)=-eta_jump*dv(b)*ddel(a,b);
- dz_jump_nt(a,t)=eta_jump*dv(a);
- dz_jump(a,b)=dz_jump_nt(c,d)*nt_to_xy(a,c)*nt_to_xy(b,d);
-
- p_jump=-2*eta_jump*dv(t);
-
- /* TODO: dp_jump and ddz_jump both need to be corrected for
- curvature and density jumps. */
-
- FTensor::Tensor1<double,2> dp_jump_nt;
- dp_jump_nt(n)=2*eta_jump*ddv(n);
- dp_jump_nt(t)=-2*eta_jump*ddv(t);
- dp_jump(a)=nt_to_xy(a,b)*dp_jump_nt(b);
-
- FTensor::Tensor3_christof<double,2,2> ddz_jump_nt;
- ddz_jump_nt(a,t,t)=eta_jump*ddv(a);
- ddz_jump_nt(a,n,n)=-ddz_jump_nt(a,t,t) + dp_jump_nt(a);
- ddz_jump_nt(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
-
- /* It would be nice if there were a way to do this directly instead
- of iterating over all of the indices. */
- ddz_jump(a,n,n)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(n,f);
- ddz_jump(a,n,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(n,e)*nt_to_xy(t,f);
- ddz_jump(a,t,t)=ddz_jump_nt(d,e,f)*nt_to_xy(a,d)*nt_to_xy(t,e)*nt_to_xy(t,f);
-}
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump,
- FTensor::Tensor3_christof<double,2,2> &ddz_jump)
-{
- double p_jump;
- FTensor::Tensor1<double,2> dp_jump;
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
- p_jump,dp_jump);
-}
-
-void compute_jumps(const double &eta_jump,
- const FTensor::Tensor1<double,2> &v,
- const FTensor::Tensor1<double,2> &dv,
- const FTensor::Tensor1<double,2> &ddv,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &z_jump,
- FTensor::Tensor2<double,2,2> &dz_jump)
-{
- FTensor::Tensor3_christof<double,2,2> ddz_jump;
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
-}
-
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,19 +0,0 @@
-#ifndef GAMR_COMPUTE_V_ON_INTERFACE_HXX
-#define GAMR_COMPUTE_V_ON_INTERFACE_HXX
-
-#include "FTensor.hpp"
-
-void compute_v_on_interface(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1],
- const FTensor::Tensor1<double,2> &pos,
- const int &xyz,
- FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor2<double,2,2> &nt_to_xy);
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/Valid.hxx
--- a/compute_v_on_interface/Valid.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-#ifndef GAMR_VALID_HXX
-#define GAMR_VALID_HXX
-
-#include <limits>
-
-class Valid
-{
-public:
- int sign,i,j;
- bool valid;
- FTensor::Tensor1<double,2> diff;
- double r;
-
- Valid(): sign(0), valid(false), r(std::numeric_limits<double>::infinity()) {}
- Valid(const int &Sign, const bool &val, const int &I, const int &J,
- const FTensor::Tensor1<double,2> &Diff):
- sign(Sign), i(I), j(J), valid(val), diff(Diff),
- r(diff(0)*diff(0)+diff(1)*diff(1)) {}
- bool operator<(const Valid &v) const
- {
- return r<v.r;
- }
- bool operator>(const Valid &v) const
- {
- return r>v.r;
- }
-};
-
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/compute_1st_derivs.hxx
--- a/compute_v_on_interface/compute_1st_derivs.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-#ifndef GAMR_COMPUTE_1ST_DERIVS_HXX
-#define GAMR_COMPUTE_1ST_DERIVS_HXX
-
-#include <cassert>
-#include "Valid.hxx"
-#include <gsl/gsl_linalg.h>
-
-template<int ny>
-void compute_1st_derivs(const int &d, const std::vector<Valid> &derivs,
- const int &d0,
- const int &nx,
- const double z[][ny],
- const double log_eta[][ny],
- FTensor::Tensor2<double,2,2> &dv_pm)
-{
- /* Create and invert the Vandermonde matrix */
- assert(derivs.size()==3);
- double m[3*3], rhs[3];
- for(int i=0;i<3;++i)
- {
- switch(d0)
- {
- case 0:
- if(derivs[i].i==nx-1)
- {
- rhs[i]=0;
- }
- else
- {
- rhs[i]=vel(z,log_eta,derivs[i].i+1,derivs[i].j)
- - vel(z,log_eta,derivs[i].i,derivs[i].j);
- }
- break;
- case 1:
- if(derivs[i].j==ny-1)
- {
- rhs[i]=0;
- }
- else
- {
- rhs[i]=vel(z,log_eta,derivs[i].i,derivs[i].j+1)
- - vel(z,log_eta,derivs[i].i,derivs[i].j);
- }
- break;
- default:
- abort();
- }
- m[0+3*i]=1;
- m[1+3*i]=derivs[i].diff(0);
- m[2+3*i]=derivs[i].diff(1);
- }
-
- gsl_matrix_view mv=gsl_matrix_view_array(m,3,3);
- gsl_vector_view rhsv=gsl_vector_view_array(rhs,3);
- gsl_vector *x=gsl_vector_alloc(3);
- int s;
- gsl_permutation *perm=gsl_permutation_alloc(3);
- gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
- gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
-
- dv_pm(d,d0)=gsl_vector_get(x,0);
- gsl_vector_free(x);
- gsl_permutation_free(perm);
-}
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/compute_2nd_derivs.hxx
--- a/compute_v_on_interface/compute_2nd_derivs.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,74 +0,0 @@
-#ifndef GAMR_COMPUTE_2ND_DERIVS_HXX
-#define GAMR_COMPUTE_2ND_DERIVS_HXX
-
-#include "Valid.hxx"
-
-template<int ny>
-void compute_2nd_derivs(const int &d, const Valid &v,
- const int &d0, const int &d1,
- const int &nx,
- const double z[][ny],
- const double log_eta[][ny],
- FTensor::Tensor3_christof<double,2,2> &ddv_pm,
- double &eta_pm)
-{
- if(d0!=d1)
- {
- if(v.i==nx-1 || v.j==ny-1)
- {
- ddv_pm(d,d0,d1)=0;
- }
- else
- {
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j+1)
- - vel(z,log_eta,v.i,v.j+1)
- - vel(z,log_eta,v.i+1,v.j)
- + vel(z,log_eta,v.i,v.j);
- }
- }
- else
- {
- switch(d0)
- {
- case 0:
- if(v.i==0)
- {
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
- - vel(z,log_eta,v.i,v.j);
- }
- else if(v.i==nx-1)
- {
- ddv_pm(d,d0,d1)=
- - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
- }
- else
- {
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
- - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
- }
- break;
- case 1:
- if(v.j==0)
- {
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
- - vel(z,log_eta,v.i,v.j);
- }
- else if(v.j==ny-1)
- {
- ddv_pm(d,d0,d1)=
- - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
- }
- else
- {
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
- - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
- }
- break;
- default:
- abort();
- }
- }
- eta_pm=exp(log_eta[v.i][v.j]);
-}
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,413 +0,0 @@
-#include "../constants.hxx"
-#include <cmath>
-#include "FTensor.hpp"
-#include <vector>
-#include "../sign.hxx"
-#include "vel.hxx"
-#include "compute_2nd_derivs.hxx"
-#include "compute_1st_derivs.hxx"
-#include "compute_values.hxx"
-#include <tuple>
-#include <algorithm>
-#include <queue>
-
-void compute_dv_dtt(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1],
- const FTensor::Tensor1<double,2> &pos,
- const FTensor::Tensor1<double,2> &norm,
- const FTensor::Tensor1<double,2> &tangent,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &v)
-{
- double length=h/std::max(std::abs(norm(0)),std::abs(norm(1)));
-
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
- const FTensor::Index<'c',2> c;
- FTensor::Tensor3_christof<double,2,2> ddv_pm[2];
- ddv_pm[0](a,b,c)=0;
- ddv_pm[1](a,b,c)=0;
-
- FTensor::Tensor2<double,2,2> dv_pm[2];
- dv_pm[0](a,b)=0;
- dv_pm[1](a,b)=0;
-
- FTensor::Tensor1<double,2> v_pm[2];
- v_pm[0](a)=0;
- v_pm[1](a)=0;
-
- double eta_pm[2];
-
- /* 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));
-
- std::vector<Valid> valid((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);
-
- int starti((pos(0)-dx)/h), startj((pos(1)-dy)/h);
-
- /* d/dx^2 */
- 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), p_d(i*h+dx+h/2,j*h+dy),
- dd_diff, d_diff, diff;
- dd_diff(a)=p(a)-pos(a);
- d_diff(a)=p_d(a)-pos(a);
-
- if(d==0)
- {
- /* ddv */
- if(i>0 && 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)
- {
- d_diff(a)-=length*norm(a)*sign(distx[i][j]);
- 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);
- }
- /* v */
- diff(a)=dd_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
- {
- /* ddv */
- 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
- 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_diff(a)-=length*norm(a)*sign(disty[i][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);
-
- /* v */
- diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(disty[i][j]);
- valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),true,i,j,diff);
- }
- }
-
- /* d/dy^2 */
- 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), 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)
- {
- /* ddv */
- bool v;
- if(j==0)
- v=sign(distx[i][j])==sign(distx[i][j+1]);
- else if(j==max_y)
- v=sign(distx[i][j])==sign(distx[i][j-1]);
- else
- v=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,i,j,dd_diff);
-
- /* dv */
- bool v_d;
- if(j==max_y)
- v_d=true;
- else
- v_d=sign(distx[i][j])==sign(distx[i][j+1]);
-
- d_diff(a)-=length*norm(a)*sign(distx[i][j]);
- 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);
- }
- else
- {
- /* ddv */
- 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);
-
- /* dv */
- if(j<max_y)
- {
- d_diff(a)-=length*norm(a)*sign(disty[i][j]);
- 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);
- }
- }
- }
-
- /* 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_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<max_x)
- {
- bool v;
- if(j==max_y)
- v=sign(distx[i][j])==sign(distx[i+1][j]);
- else
- v=sign(distx[i][j])==sign(distx[i][j+1])
- && sign(distx[i][j])==sign(distx[i+1][j])
- && sign(distx[i][j])==sign(distx[i+1][j+1]);
-
- 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);
- }
- }
- else
- {
- if(j<max_y)
- {
- bool v;
- if(i==max_x)
- v=sign(distx[i][j])==sign(distx[i][j+1]);
- else
- v=sign(distx[i][j])==sign(distx[i][j+1])
- && sign(distx[i][j])==sign(distx[i+1][j])
- && sign(distx[i][j])==sign(distx[i+1][j+1]);
-
- 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);
- }
- }
- }
-
- /* Put the valid vectors into a priority queue so that we can
- take the sorted versions out as needed. A bit messier than
- just sorting the vector, but also significantly faster. */
-
- FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
- false,false);
-
- typedef std::priority_queue<Valid, std::vector<Valid>, std::greater<Valid> > mpq;
- mpq ddp_valid[]={{dd_valid(0,0).begin(),dd_valid(0,0).end()},
- {dd_valid(0,1).begin(),dd_valid(0,1).end()},
- {dd_valid(1,1).begin(),dd_valid(1,1).end()}};
-
- mpq dp_valid[]={{d_valid(0).begin(),d_valid(0).end()},
- {d_valid(1).begin(),d_valid(1).end()}};
- std::priority_queue<Valid, std::vector<Valid>, std::greater<Valid> >
- p_valid(valid.begin(),valid.end());
-
- /* Compute everything */
- /* v */
-
- std::vector<Valid> values[2];
- for(auto V=p_valid.top(); !p_valid.empty(); p_valid.pop(), V=p_valid.top())
- {
- if(values[0].size()==6 && values[1].size()==6)
- break;
-
- if(V.sign==0)
- continue;
- const int pm=(V.sign==-1 ? 0 : 1);
- switch(values[pm].size())
- {
- case 0:
- case 1:
- case 2:
- values[pm].push_back(V);
- break;
- case 6:
- break;
- default:
- int x(0),y(0);
- std::vector<int> x_coords={V.i};
- std::vector<int> y_coords={V.j};
- for(auto &value: values[pm])
- {
- if(V.i==value.i)
- ++x;
- if(V.j==value.j)
- ++y;
- if(values[pm].size()==5)
- {
- if(find(x_coords.begin(),x_coords.end(),value.i)==x_coords.end())
- x_coords.push_back(value.i);
- if(find(y_coords.begin(),y_coords.end(),value.j)==y_coords.end())
- y_coords.push_back(value.j);
- }
- }
- if(x<3 && y<3 && !(values[pm].size()==5 &&
- (x_coords.size()<3 || y_coords.size()<3)))
- values[pm].push_back(V);
- 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]);
- else
- compute_values(1,values[pm],Nx,zy,log_etay,
- norm,pm==0 ? -length : length,v_pm[pm]);
- }
- break;
- }
- }
-
- for(int d0=0;d0<2;++d0)
- {
- /* First derivative dv */
- std::vector<Valid> derivs[2];
- for(auto V=dp_valid[d0].top(); !dp_valid[d0].empty();
- dp_valid[d0].pop(), V=dp_valid[d0].top())
- {
- if(derivs[0].size()==3 && derivs[1].size()==3)
- break;
-
- if(V.sign==0)
- continue;
- if(V.valid)
- {
- const int pm=(V.sign==-1 ? 0 : 1);
- switch(derivs[pm].size())
- {
- default:
- break;
- case 0:
- case 1:
- derivs[pm].push_back(V);
- break;
- case 2:
- Valid &v0(derivs[pm][0]);
- Valid &v1(derivs[pm][1]);
- if(!((v0.i==v1.i && v1.i==V.i)
- || (v0.j==v1.j && v1.j==V.j)))
- {
- derivs[pm].push_back(V);
- if(d==0)
- compute_1st_derivs(0,derivs[pm],d0,Nx+1,
- zx,log_etax,dv_pm[pm]);
- else
- compute_1st_derivs(1,derivs[pm],d0,Nx,
- zy,log_etay,dv_pm[pm]);
- }
- break;
- }
- }
- }
- /* Second derivative ddv */
- for(int d1=d0;d1<2;++d1)
- {
- for(auto V=ddp_valid[2*d0+d1].top(); !ddp_valid[2*d0+d1].empty();
- ddp_valid[2*d0+d1].pop(), V=ddp_valid[2*d0+d1].top())
- {
- if(is_set(0,d0,d1) && is_set(1,d0,d1))
- break;
-
- for(int pm=0;pm<2;++pm)
- {
- if(V.sign==(pm==0 ? -1 : 1) && V.valid && !is_set(pm,d0,d1))
- {
- is_set(pm,d0,d1)=true;
- if(d==0)
- compute_2nd_derivs(0,V,d0,d1,Nx+1,zx,log_etax,
- ddv_pm[pm],eta_pm[pm]);
- else
- compute_2nd_derivs(1,V,d0,d1,Nx,zy,log_etay,
- ddv_pm[pm],eta_pm[pm]);
- }
- }
- }
- }
- }
- }
-
- /* ddv */
- double temp(0);
- FTensor::Tensor1<double,2> ddv_xy(0,0);
- for(int d=0;d<2;++d)
- {
- temp+=eta_pm[d]*eta_pm[d];
- ddv_xy(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
- }
- ddv_xy(a)/=temp*h*h;
-
- ddv(a)=nt_to_xy(b,a)*ddv_xy(b);
-
- /* dv */
- FTensor::Tensor1<double,2> dv_xy(0,0);
- for(int d=0;d<2;++d)
- dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
-
- dv(a)=nt_to_xy(b,a)*dv_xy(b);
-
- FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
- FTensor::Tensor1<double,2> ddz_nt_jump;
- double eta_jump(eta_pm[1]-eta_pm[0]);
- /* TODO: fix for curved surfaces */
- ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
-
- dv(a)-=h*length*ddz_nt_jump(a);
- dv(a)/=(eta_pm[0]+eta_pm[1])*h;
-
- /* v */
- FTensor::Tensor1<double,2> v_xy(0,0);
- for(int d=0;d<2;++d)
- v_xy(a)+=v_pm[d](a)*eta_pm[d];
-
- v(a)=nt_to_xy(b,a)*v_xy(b);
-
- /* TODO: fix for curved surfaces */
- FTensor::Tensor1<double,2> dz_n_jump;
- dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
-
- v(a)-=2*length*dz_n_jump(a)/3;
- v(a)/=(eta_pm[0] + eta_pm[1]);
-}
-
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-#include <algorithm>
-#include <cmath>
-#include "../constants.hxx"
-#include <limits>
-#include <iostream>
-#include "FTensor.hpp"
-#include "../compute_v_on_interface.hxx"
-#include "vel.hxx"
-
-void compute_dv_dtt(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1],
- const FTensor::Tensor1<double,2> &pos,
- const FTensor::Tensor1<double,2> &norm,
- const FTensor::Tensor1<double,2> &tangent,
- const FTensor::Tensor2<double,2,2> &nt_to_xy,
- FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &v);
-
-void compute_v_on_interface(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1],
- const FTensor::Tensor1<double,2> &pos,
- const int &xyz,
- FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &dv,
- FTensor::Tensor1<double,2> &ddv,
- 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));
-
- tangent(0)=-norm(1);
- tangent(1)=norm(0);
-
- FTensor::Number<0> n;
- FTensor::Number<1> t;
- nt_to_xy(a,n)=norm(a);
- nt_to_xy(a,t)=tangent(a);
-
- compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,nt_to_xy,
- ddv,dv,v);
-}
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/compute_values.hxx
--- a/compute_v_on_interface/compute_values.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,61 +0,0 @@
-#ifndef GAMR_COMPUTE_VALUES_HXX
-#define GAMR_COMPUTE_VALUES_HXX
-
-#include <cassert>
-#include "Valid.hxx"
-#include <gsl/gsl_linalg.h>
-
-template<int ny>
-void compute_values(const int &d, const std::vector<Valid> &values,
- const int &nx,
- const double z[][ny],
- const double log_eta[][ny],
- const FTensor::Tensor1<double,2> &norm,
- const double &signed_length,
- FTensor::Tensor1<double,2> &v_pm)
-{
- /* Create and invert the Vandermonde matrix */
- assert(values.size()==6);
- double m[6*6], rhs[6];
- for(int i=0;i<6;++i)
- {
- rhs[i]=vel(z,log_eta,values[i].i,values[i].j);
-
- m[0+6*i]=1;
- m[1+6*i]=values[i].diff(0);
- m[2+6*i]=values[i].diff(1);
- m[3+6*i]=values[i].diff(0)*values[i].diff(0)/2;
- m[4+6*i]=values[i].diff(1)*values[i].diff(1)/2;
- m[5+6*i]=values[i].diff(0)*values[i].diff(1);
- }
-
- gsl_matrix_view mv=gsl_matrix_view_array(m,6,6);
- gsl_vector_view rhsv=gsl_vector_view_array(rhs,6);
- gsl_vector *x=gsl_vector_alloc(6);
- int s;
- gsl_permutation *perm=gsl_permutation_alloc(6);
- gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
- gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
-
- double v=gsl_vector_get(x,0);
-
- FTensor::Tensor1<double,2> dv_xy;
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
- dv_xy(0)=gsl_vector_get(x,1);
- dv_xy(1)=gsl_vector_get(x,2);
- double dv=dv_xy(a)*norm(a)*signed_length/2;
-
- FTensor::Tensor2_symmetric<double,2> ddv_xy;
- ddv_xy(0,0)=gsl_vector_get(x,3);
- ddv_xy(1,1)=gsl_vector_get(x,4);
- ddv_xy(0,1)=gsl_vector_get(x,5);
- double ddv=ddv_xy(a,b)*norm(a)*norm(b)*signed_length*signed_length/4;
-
- v_pm(d)=(3*v - 5*dv + 1.5*ddv)/3;
-
- gsl_vector_free(x);
- gsl_permutation_free(perm);
-}
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 compute_v_on_interface/vel.hxx
--- a/compute_v_on_interface/vel.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-#ifndef GAMR_VEL_HXX
-#define GAMR_VEL_HXX
-
-inline double analytic(const double x, const double y, const bool return_ux)
-{
- const double eta_jump=max_eta-min_eta;
-
- double ux(1.1), ux_xp(1.2), ux_xxp(1.3), ux_xyp(1.4),
- ux_xm(1.5), ux_xxm(1.6), ux_xym(1.7),
- uy(2.1), uy_xp(2.2), uy_xxp(2.3), uy_xyp(2.4),
- uy_xm(2.5), uy_xxm(2.6), uy_xym(2.7);
- double ux_y(-(max_eta*uy_xp - min_eta*uy_xm)/eta_jump),
- ux_yy(-(max_eta*uy_xyp - min_eta*uy_xym)/eta_jump),
- uy_y(-(max_eta*ux_xp - min_eta*ux_xm)/eta_jump),
- uy_yy(-(max_eta*ux_xyp - min_eta*ux_xym)/eta_jump);
-
- // std::cout << "u y "
- // << ux << " "
- // << ux_y << " "
- // << ux_yy << " "
- // << uy << " "
- // << uy_y << " "
- // << uy_yy << " "
- // << "\n";
- // exit(0);
-
- if(return_ux)
- {
- if(x>0)
- return ux + x*ux_xp + x*x*ux_xxp/2 + x*y*ux_xyp
- + y*ux_y + y*y*ux_yy/2;
- else
- return ux + x*ux_xm + x*x*ux_xxm/2 + x*y*ux_xym
- + y*ux_y + y*y*ux_yy/2;
- }
- else
- {
- if(x>0)
- return uy + x*uy_xp + x*x*uy_xxp/2 + x*y*uy_xyp
- + y*uy_y + y*y*uy_yy/2;
- else
- return uy + x*uy_xm + x*x*uy_xxm/2 + x*y*uy_xym
- + y*uy_y + y*y*uy_yy/2;
- }
-
-}
-
-
-template<int ny>
-double vel(const double z[][ny], const double log_eta[][ny],
- const int i, const int j)
-{
- // if(j>ny/8-4 && j<ny/8+4)
- // {
- // double x,y;
- // if(ny%2==0)
- // {
- // x=i*h-middle;
- // y=(j+0.5)*h-1.0/8;
- // // y=j*h-1.0/8;
- // return analytic(x,y,true);
- // }
- // else
- // {
- // x=(i+0.5)*h-middle;
- // y=j*h-1.0/8;
- // // y=(j-0.5)*h-1.0/8;
- // return analytic(x,y,false);
- // }
- // }
-
- return z[i][j]*std::exp(-log_eta[i][j]);
-}
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 sign.hxx
--- a/sign.hxx Fri Mar 30 09:31:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,8 +0,0 @@
-#ifndef GAMR_SIGN_HXX
-#define GAMR_SIGN_HXX
-
-template <typename T> int sign(T val) {
- return (val > T(0)) - (val < T(0));
-}
-
-#endif
diff -r da0a37975701 -r 7d5b5ce74b59 wscript
--- a/wscript Fri Mar 30 09:31:23 2012 -0700
+++ b/wscript Fri Mar 30 09:31:59 2012 -0700
@@ -9,15 +9,15 @@ def build(bld):
source = ['main.cxx',
'initial.cxx',
'compute_coefficients/compute_coefficients.cxx',
+ 'compute_coefficients/compute_Cxyz.cxx',
+ 'compute_coefficients/compute_Cp.cxx',
'compute_coefficients/simplified_Rx.cxx',
'compute_coefficients/simplified_Ry.cxx',
'compute_coefficients/simplified_Rp.cxx',
- 'compute_v_on_interface/compute_v_on_interface.cxx',
- 'compute_v_on_interface/compute_dv_dtt.cxx',
- 'compute_Cxyz.cxx',
- 'compute_Cp.cxx',
- 'compute_jumps.cxx',
- 'Interface/Interface.cxx'],
+ 'compute_coefficients/compute_v_on_interface/compute_v_on_interface.cxx',
+ 'compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx',
+ 'compute_coefficients/compute_jumps.cxx',
+ 'compute_coefficients/Interface/Interface.cxx'],
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-Wall','-Drestrict=','-DFTENSOR_DEBUG'],
More information about the CIG-COMMITS
mailing list