[cig-commits] commit: Add support for anticorners. Change "intersect_xy" to
Mercurial
hg at geodynamics.org
Sun Mar 25 16:28:25 PDT 2012
changeset: 114:2690d006b243
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sun Mar 25 16:28:14 2012 -0700
files: Interface.hxx Interface/Interface.cxx compute_Cxyz.cxx compute_v_on_interface/compute_v_on_interface.cxx
description:
Add support for anticorners. Change "intersect_xy" to
"intersect_sides" and compute_xy to compute_corners. Simplify the
computation of the norm since it is only first order anyway.
diff -r 644a357572e3 -r 2690d006b243 Interface.hxx
--- a/Interface.hxx Fri Mar 23 09:50:30 2012 -0700
+++ b/Interface.hxx Sun Mar 25 16:28:14 2012 -0700
@@ -4,13 +4,16 @@
#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];
+ FTensor::Tensor1<double,2> grid_pos, pos[2], corner_pos[2][2],
+ anticorner_pos[2];
int corner_dir;
- bool intersect_dd[2], intersect_xy[2], intersect_corner[2][2];
+ 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,
@@ -19,54 +22,63 @@ public:
bool intersects() const
{
return intersect_dd[0] || intersect_dd[1]
- || intersect_xy[0] || intersect_xy[1] || intersect_dp
+ || 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_corners(const int &i, const int &j, const double dist[][ny])
+ template <int ny> void compute_xy(const int &i, const int &j,
+ const int ¢er_dist,
+ const double dist[][ny])
{
- intersect_xy[0]=(sign(dist[i][j])!=sign(dist[i+1][j])
- && sign(dist[i][j+1])!=sign(dist[i+1][j+1]));
- intersect_xy[1]=(sign(dist[i][j])!=sign(dist[i][j+1])
- && sign(dist[i+1][j])!=sign(dist[i+1][j+1]));
+ 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]=
- (sign(dist[i+n0][j+n1])!=sign(dist[i+(n0+1)%2][j+n1])
- && sign(dist[i+n0][j+n1])!=sign(dist[i][j+(n1+1)%2]));
+ 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;
- if(intersect_corner[n0][n1])
- {
- double dh[2]={sign(n0-0.5)*h/2, sign(n1-0.5)*h/2};
- double delta;
+ 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);
- /* TODO fix this for higher order curved surfaces */
-
- /* Correct in the direction that is closest to the
- corner. This should use the interface that is more
- perpendicular to the derivative direction. */
- if(std::fabs(dist[i+(n0+1)%2][j+n1])>std::fabs(dist[i+n0][j+(n1+1)%2]))
- {
- corner_dir=0;
- delta=dist[i+n0][j+(n1+1)%2]-dist[i+n0][j+n1];
- corner_pos[n0][n1](0)=grid_pos(0)+dh[0]-dist[i+n0][j+n1]*h/delta;
- corner_pos[n0][n1](1)=grid_pos(1)+dh[1];
- }
- else
- {
- corner_dir=1;
- delta=dist[i+n0][j+(n1+1)%2]-dist[i+n0][j+n1];
- corner_pos[n0][n1](0)=grid_pos(0)+dh[0];
- corner_pos[n0][n1](1)=grid_pos(1)+dh[1]-dist[i+n0][j+n1]*h/delta;
- }
-
-
- }
- }
- }
+ std::cerr << "Anticorners not tested yet"
+ << std::endl;
+ abort();
+ }
+ }
}
};
diff -r 644a357572e3 -r 2690d006b243 Interface/Interface.cxx
--- a/Interface/Interface.cxx Fri Mar 23 09:50:30 2012 -0700
+++ b/Interface/Interface.cxx Sun Mar 25 16:28:14 2012 -0700
@@ -59,7 +59,7 @@ Interface::Interface(const int &d, const
pos[1](1)=grid_pos(1)-distx[i][j]*h/delta;
}
- compute_corners(i-1,j,disty);
+ compute_xy(i-1,j,distx[i][j],disty);
}
else
{
@@ -102,6 +102,6 @@ Interface::Interface(const int &d, const
pos[1](1)=grid_pos(1)-disty[i][j]*h/delta;
}
- compute_corners(i,j-1,distx);
+ compute_xy(i,j-1,disty[i][j],distx);
}
}
diff -r 644a357572e3 -r 2690d006b243 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Fri Mar 23 09:50:30 2012 -0700
+++ b/compute_Cxyz.cxx Sun Mar 25 16:28:14 2012 -0700
@@ -24,6 +24,7 @@ void compute_Cxyz(const double zx[Nx+1][
const FTensor::Index<'c',2> c;
pos(a)=interface.grid_pos(a);
+ bool try_corners(true);
for(int dd=0;dd<2;++dd)
{
if(interface.intersect_dd[dd])
@@ -86,8 +87,9 @@ void compute_Cxyz(const double zx[Nx+1][
C+=(p_jump + dx*dp_jump(a)*dir(a))/h;
}
- if(interface.intersect_xy[dd])
+ 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;
@@ -99,19 +101,26 @@ void compute_Cxyz(const double zx[Nx+1][
+ dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
}
}
- // /* Handle corner cutting */
- // for(int ee=0;ee<2;++ee)
- // {
- // if(interface.intersect_corner[dd][ee])
- // {
- // FTensor::Tensor1<double,2> v, dv, ddv;
- // compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- // interface.corner_pos[dd][ee],d,
- // v,dv,ddv);
+ }
+ 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)
+ 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;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface.corner_pos[dd][ee],d,
+ v,dv,ddv);
- // double dh[2]={sign(dd-0.5)*h/2, sign(ee-0.5)*h/2};
-
- // }
- // }
- }
+ }
}
diff -r 644a357572e3 -r 2690d006b243 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Fri Mar 23 09:50:30 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Sun Mar 25 16:28:14 2012 -0700
@@ -51,24 +51,21 @@ void compute_v_on_interface(const double
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);
- FTensor::Tensor2<double,2,2> dx((pos(0)-ix*h)/h,(pos(0)-iy*h)/h-0.5,
- (pos(1)-jx*h)/h-0.5,(pos(1)-jy*h)/h);
+ // FTensor::Tensor2<double,2,2> dx((pos(0)-ix*h)/h,(pos(0)-iy*h)/h-0.5,
+ // (pos(1)-jx*h)/h-0.5,(pos(1)-jy*h)/h);
- /* TODO: fix for curved interfaces */
+ /* 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])*(1-dx(0,xyz))
- + (disty[iy+1][jy+1] - disty[iy][jy+1])*dx(0,xyz))/h;
- norm(1)=((disty[iy][jy+1] - disty[iy][jy])*(1-dx(1,xyz))
- + (disty[iy+1][jy+1] - disty[iy+1][jy])*dx(1,xyz))/h;
+ 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])*(1-dx(0,xyz))
- + (distx[ix+1][jx+1] - distx[ix][jx+1])*dx(0,xyz))/h;
- norm(1)=((distx[ix][jx+1] - distx[ix][jx])*(1-dx(1,xyz))
- + (distx[ix+1][jx+1] - distx[ix+1][jx])*dx(1,xyz))/h;
+ 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));
More information about the CIG-COMMITS
mailing list