aboutsummaryrefslogtreecommitdiff
path: root/src/patch
diff options
context:
space:
mode:
authorjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2003-03-17 12:06:49 +0000
committerjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2003-03-17 12:06:49 +0000
commit8d274d170dde9d858e0a13a43762d3976ea4e3b1 (patch)
treed394b0725cf7af6883c57b8bde59f0ffef479d6c /src/patch
parent498e26dd4ee819db7a0c4930e73624af936b7d47 (diff)
add 1st draft of support for setting an excision mask
-- alas this doesn't work properly (yet) for multiple processors git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@975 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src/patch')
-rw-r--r--src/patch/coords.hh17
-rw-r--r--src/patch/ghost_zone.cc2
-rw-r--r--src/patch/ghost_zone.hh2
-rw-r--r--src/patch/patch.hh12
-rw-r--r--src/patch/patch_interp.hh11
-rw-r--r--src/patch/patch_system.cc369
-rw-r--r--src/patch/patch_system.hh69
7 files changed, 324 insertions, 158 deletions
diff --git a/src/patch/coords.hh b/src/patch/coords.hh
index 8c16f5b..560d1d6 100644
--- a/src/patch/coords.hh
+++ b/src/patch/coords.hh
@@ -249,6 +249,7 @@ inline
class global_coords
{
public:
+ #ifdef NOT_USED
// global (x,y,z) --> local (x,y,z)
fp local_x_of_global_x(fp global_x) const
{ return global_x - origin_x_; }
@@ -256,7 +257,9 @@ public:
{ return global_y - origin_y_; }
fp local_z_of_global_z(fp global_z) const
{ return global_z - origin_z_; }
+ #endif /* NOT_USED */
+ #ifdef NOT_USED
// local (x,y,z) --> global (x,y,z)
fp global_x_of_local_x(fp local_x) const
{ return origin_x_ + local_x; }
@@ -264,12 +267,26 @@ public:
{ return origin_y_ + local_y; }
fp global_z_of_local_z(fp local_z) const
{ return origin_z_ + local_z; }
+ #endif /* NOT_USED */
// get global (x,y,z) coordinates of local origin point
fp origin_x() const { return origin_x_; }
fp origin_y() const { return origin_y_; }
fp origin_z() const { return origin_z_; }
+ // radius of given (x,y,z) with respect to local origin point
+ #ifdef NOT_USED
+ fp radius_of_local_xyz(fp local_x, fp local_y, fp local_z) const
+ { return jtutil::hypot3(local_x, local_y, local_z); }
+ fp radius_of_global_xyz(fp global_x, fp global_y, fp global_z)
+ const
+ {
+ return radius_of_local_xyz(local_x_of_global_x(global_x),
+ local_y_of_global_y(global_y),
+ local_z_of_global_z(global_z));
+ }
+ #endif /* NOT_USED */
+
// constructor: specify global (x,y,z) coordinates of local origin point
global_coords(fp origin_x_in, fp origin_y_in, fp origin_z_in)
: origin_x_(origin_x_in),
diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 922f0ff..fc61220 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -417,7 +417,7 @@ return max_par_adjacent_ghost_zone().is_symmetry()
// This function finishes the construction/setup of an interpatch_ghost_zone
// object. It
// - sets up the par coordinate mapping information
-// - sets up the interpolator data pointer and result arrays
+// - sets up the interpatch interpolator data pointer and result arrays
// - constructs the patch_interp object to interpolate from the *other* patch
//
// We use our ipar as the patch_interp's parindex.
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 6402a4b..ae98286 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -725,7 +725,7 @@ public:
// finish setup (requires adjacent-side ghost_zone objects
// to exist, though not to have finish_setup() called):
// - setup par coordinate mapping information
- // - setup interpolator data pointers & result buffer
+ // - setup interpatch interpolator data pointers & result buffer
// - create patch_interp object to interpolate from *other* patch
void finish_setup(int interp_handle, int interp_par_table_handle);
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 8060f5e..5002aa7 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -661,9 +661,9 @@ public:
const
{ local_coords::xyz_of_r_mu_nu(r,rho,sigma, x,y,z); }
fp rho_of_xyz(fp x, fp y, fp z) const
- { return local_coords::mu_of_yz(y,z); }
+ { return modulo_reduce_rho(local_coords::mu_of_yz(y,z)); }
fp sigma_of_xyz(fp x, fp y, fp z) const
- { return local_coords::nu_of_xz(x,z); }
+ { return modulo_reduce_sigma(local_coords::nu_of_xz(x,z)); }
// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
// with respect to the local coordinate system
@@ -806,9 +806,9 @@ public:
const
{ local_coords::xyz_of_r_nu_phi(r, rho, sigma, x, y, z); }
fp rho_of_xyz(fp x, fp y, fp z) const
- { return local_coords::nu_of_xz(x, z); }
+ { return modulo_reduce_rho(local_coords::nu_of_xz(x, z)); }
fp sigma_of_xyz(fp x, fp y, fp z) const
- { return local_coords::phi_of_xy(x, y); }
+ { return modulo_reduce_sigma(local_coords::phi_of_xy(x, y)); }
// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
// with respect to the local coordinate system
@@ -953,9 +953,9 @@ public:
const
{ local_coords::xyz_of_r_mu_phi(r, rho, sigma, x, y, z); }
fp rho_of_xyz(fp x, fp y, fp z) const
- { return local_coords::mu_of_yz(y, z); }
+ { return modulo_reduce_rho(local_coords::mu_of_yz(y, z)); }
fp sigma_of_xyz(fp x, fp y, fp z) const
- { return local_coords::phi_of_xy(x, y); }
+ { return modulo_reduce_sigma(local_coords::phi_of_xy(x, y)); }
// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
// with respect to the local coordinate system
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index b4bff04..e252693 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -217,13 +217,14 @@ public:
// ok_to_use_[min,max]_par_ghost_zone
// = Boolean flags saying whether or not we should use gridfn
// data from the [min,max]_par ghost zones in the interpolation.
- // interp_handle_in = Cactus handle to the interpolation operator.
+ // interp_handle_in = Cactus handle to the interpatch interpolation
+ // operator.
// interp_par_table_handle_in
// = Cactus handle to a Cactus key/value table giving
- // parameters (eg order) for the interpolation operator.
- // This class internally clones this table and modifies
- // the clone, so the original table is not modified by
- // any actions of this class.
+ // parameters (eg order) for the interpatch interpolation
+ // operator. This class internally clones this table and
+ // modifies the clone, so the original table is not modified
+ // by any actions of this class.
//
// This constructor requires that this patch's gridfns already
// exist, since we size various arrays based on the patch's min/max
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index d67ca02..10cbd5b 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -24,8 +24,6 @@
// patch_system::plus_or_minus_xyz_patch
// patch_system::patch_number_of_name
//
-// patch_system::patch_containing_local_xyz
-//
// patch_system::patch_irho_isigma_of_gpn
// patch_system::ghosted_patch_irho_isigma_of_gpn
//
@@ -37,6 +35,9 @@
//
// patch_system::circumference
// patch_system::integrate_gridfn
+//
+// patch_system::patch_containing_local_xyz
+// patch_system::radius_in_local_xyz_direction
//
// patch_system::print_unknown_gridfn
// patch_system::read_unknown_gridfn
@@ -104,6 +105,18 @@ using jtutil::error_exit;
// It's a fatal error (error_exit()) if this
// doesn't evenly divide the grid sizes in both
// directions.
+// ip_interp_handle = Cactus handle to the interpatch interpolation operator;
+// this must be a 1-D interpolator
+// ip_interp_par_table_handle = Cactus handle to the parameter table for the
+// interpatch interpolation operator
+// surface_interp_handle = Cactus handle to the surface interpolation
+// operator; this is optional, and is only used by
+// radius_in_{local,global}_xyz_direction()
+// If this is used, it must be a 2-D interpolator
+// surface_interp_par_table_handle = Cactus handle to the parameter table
+// for the surface interpolation operator;
+// this is optional, and is only used by
+// radius_in_local_xyz_direction()
// print_summary_msg_flag = true to print 2 lines of CCTK_VInfo messages
// giving the patch system type and origin
// false to skip this
@@ -118,7 +131,9 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
int N_zones_per_right_angle,
int min_gfn_in, int max_gfn_in,
int ghosted_min_gfn_in, int ghosted_max_gfn_in,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
+ int surface_interp_handle_in,
+ int surface_interp_par_table_handle_in,
bool print_summary_msg_flag,
bool print_detailed_msg_flag)
@@ -128,8 +143,12 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
all_patches_(new patch*[N_patches_]),
starting_gpn_(new int[N_patches_+1]),
ghosted_starting_gpn_(new int[N_patches_+1]),
- gridfn_storage_(NULL), // set in setup_gridfn_storage()
- ghosted_gridfn_storage_(NULL) // set in setup_gridfn_storage()
+ gridfn_storage_(NULL), // set in setup_gridfn_storage()
+ ghosted_gridfn_storage_(NULL), // set in setup_gridfn_storage()
+ global_min_ym_(0), global_max_ym_(0),
+ // set in compute_synchronize_Jacobian()
+ surface_interp_handle_ (surface_interp_handle_in),
+ surface_interp_par_table_handle_(surface_interp_par_table_handle_in)
{
if (! jtutil::is_odd(patch_overlap_width))
then error_exit(ERROR_EXIT,
@@ -174,8 +193,8 @@ case patch_system__full_sphere:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__full_sphere(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -188,8 +207,8 @@ case patch_system__plus_z_hemisphere:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_z_hemisphere(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -202,8 +221,8 @@ case patch_system__plus_xy_quadrant_mirrored:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_xy_quadrant_mirrored(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -216,8 +235,8 @@ case patch_system__plus_xy_quadrant_rotating:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_xy_quadrant_rotating(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -230,8 +249,8 @@ case patch_system__plus_xz_quadrant_rotating:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_xz_quadrant_rotating(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -244,8 +263,8 @@ case patch_system__plus_xyz_octant_mirrored:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_xyz_octant_mirrored(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -258,8 +277,8 @@ case patch_system__plus_xyz_octant_rotating:
ghosted_min_gfn_in, ghosted_max_gfn_in,
print_detailed_msg_flag);
setup_ghost_zones__plus_xyz_octant_rotating(patch_overlap_width,
- interp_handle,
- interp_par_table_handle,
+ ip_interp_handle,
+ ip_interp_par_table_handle,
print_detailed_msg_flag);
break;
@@ -620,7 +639,7 @@ const patch& plast = ith_patch(N_patches()-1);
//
void patch_system::setup_ghost_zones__full_sphere
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -657,40 +676,40 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, mx,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(py, mx,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mx, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(my, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, mx,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -703,7 +722,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_z_hemisphere
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -739,28 +758,28 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, mx,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(py, mx,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mx, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(my, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -773,7 +792,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_xy_quadrant_mirrored
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -807,19 +826,19 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -832,7 +851,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_xy_quadrant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -869,19 +888,19 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(mz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -894,7 +913,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_xz_quadrant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -931,19 +950,19 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, my,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -956,7 +975,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_xyz_octant_mirrored
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -987,13 +1006,13 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -1006,7 +1025,7 @@ assert_all_ghost_zones_fully_setup();
//
void patch_system::setup_ghost_zones__plus_xyz_octant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag)
{
if (print_msg_flag)
@@ -1039,13 +1058,13 @@ if (print_msg_flag)
" finishing interpatch setup");
finish_interpatch_setup(pz, px,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(pz, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
finish_interpatch_setup(px, py,
patch_overlap_width,
- interp_handle, interp_par_table_handle);
+ ip_interp_handle, ip_interp_par_table_handle);
assert_all_ghost_zones_fully_setup();
}
@@ -1103,16 +1122,16 @@ py.create_interpatch_ghost_zone(ey, ex, patch_overlap_width);
void patch_system::finish_interpatch_setup
(patch &px, patch &py,
int patch_overlap_width,
- int interp_handle, int interp_par_table_handle)
+ int ip_interp_handle, int ip_interp_par_table_handle)
{
const patch_edge& ex = px.edge_adjacent_to_patch(py, patch_overlap_width);
const patch_edge& ey = py.edge_adjacent_to_patch(px, patch_overlap_width);
px.ghost_zone_on_edge(ex)
.cast_to_interpatch_ghost_zone()
- .finish_setup(interp_handle, interp_par_table_handle);
+ .finish_setup(ip_interp_handle, ip_interp_par_table_handle);
py.ghost_zone_on_edge(ey)
.cast_to_interpatch_ghost_zone()
- .finish_setup(interp_handle, interp_par_table_handle);
+ .finish_setup(ip_interp_handle, ip_interp_par_table_handle);
}
//******************************************************************************
@@ -1270,57 +1289,6 @@ error_exit(ERROR_EXIT,
//******************************************************************************
//
-// This function finds what patch contains (the ray from the origin to)
-// a given local (x,y,z) position.
-//
-// If there are multiple patches containing the position, we return the
-// one which would still contain it if patches didn't overlap; if multiple
-// patches satisfy this criterion then it's arbitrary which one we return.
-//
-// If no patch contains the position (this can only happen for a
-// non--full-sphere patch system), then we do an error_exit() (and
-// don't return to the caller).
-//
-// Arguments:
-// (x,y,z) = The local coordinates to be converted.
-//
-// Results:
-// This function returns a reference to the containing patch.
-//
-const patch& patch_system::patch_containing_local_xyz(fp x, fp y, fp z)
- const
-{
-if ((x == 0.0) && (y == 0.0) && (z == 0.0))
- then error_exit(ERROR_EXIT,
-"***** patch_system::patch_containing_local_xyz(): position is at the origin!");
- /*NOTREACHED*/
-
-// to which axis is (x,y,z) closest?
-// ... or equivalently, which of |x|, |y|, and |z| is largest?
-const fp abs_x = jtutil::abs(x);
-const fp abs_y = jtutil::abs(y);
-const fp abs_z = jtutil::abs(z);
-
-if ((abs_z >= abs_x) && (abs_z >= abs_y))
- then return plus_or_minus_xyz_patch(z > 0.0, 'z'); // +/- z patch
-else if ((abs_x >= abs_y) && (abs_x >= abs_z))
- then return plus_or_minus_xyz_patch(x > 0.0, 'x'); // +/- x patch
-else if ((abs_y >= abs_x) && (abs_y >= abs_z))
- then return plus_or_minus_xyz_patch(y > 0.0, 'y'); // +/- y patch
-else error_exit(ERROR_EXIT,
-"***** patch_system::patch_containing_local_xyz():\n"
-" unknown (wierd!) ordering of |x|, |y|, and |z|!\n"
-" (this should never happen!)\n"
-" (x,y,z)=(%g,%g,%g)\n"
-,
- double(x), double(y), double(z));
-}
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
// This function decodes a 0-origin grid point number into a
// (patch,irho,isigma) triple.
//
@@ -1709,6 +1677,159 @@ return integral;
//******************************************************************************
//
+// This function finds what patch contains (the ray from the origin to)
+// a given local (x,y,z) position.
+//
+// If there are multiple patches containing the position, we return the
+// one which would still contain it if patches didn't overlap; if multiple
+// patches satisfy this criterion then it's arbitrary which one we return.
+//
+// If no patch contains the position (this can only if the point as at
+// the local coordinate origin, or for a non--full-sphere patch system),
+// then we return a NULL pointer.
+//
+// Arguments:
+// (x,y,z) = The local coordinates to be converted.
+//
+// Results:
+// This function returns a reference to the containing patch.
+//
+const patch* patch_system::patch_containing_local_xyz(fp x, fp y, fp z)
+ const
+{
+if ((x == 0.0) && (y == 0.0) && (z == 0.0))
+ then return NULL;
+
+// to which axis is (x,y,z) closest?
+// ... or equivalently, which of |x|, |y|, and |z| is largest?
+const fp abs_x = jtutil::abs(x);
+const fp abs_y = jtutil::abs(y);
+const fp abs_z = jtutil::abs(z);
+
+if ((abs_z >= abs_x) && (abs_z >= abs_y))
+ then return &plus_or_minus_xyz_patch(z > 0.0, 'z'); // +/- z patch
+else if ((abs_x >= abs_y) && (abs_x >= abs_z))
+ then return &plus_or_minus_xyz_patch(x > 0.0, 'x'); // +/- x patch
+else if ((abs_y >= abs_x) && (abs_y >= abs_z))
+ then return &plus_or_minus_xyz_patch(y > 0.0, 'y'); // +/- y patch
+else error_exit(ERROR_EXIT,
+"***** patch_system::patch_containing_local_xyz():\n"
+" unknown (wierd!) ordering of |x|, |y|, and |z|!\n"
+" (this should never happen!)\n"
+" [local] (x,y,z)=(%g,%g,%g)\n"
+,
+ double(x), double(y), double(z));
+}
+
+//******************************************************************************
+
+//
+// This function computes the radius of a patch-system 2-surface in the
+// direction of a specified (x,y,z) point. (If the point coincides with
+// the local origin, we return the dummy value 1.0.)
+//
+// Bugs:
+// Due to the surface-interpolator overhead, repeatedly calling this
+// function is rather inefficient.
+//
+fp patch_system::radius_in_local_xyz_direction(int ghosted_radius_gfn,
+ fp x, fp y, fp z)
+ const
+{
+if ((x == 0.0) && (y == 0.0) && (z == 0.0))
+ then return 1.0;
+
+const patch* p_ptr = patch_containing_local_xyz(x, y, z);
+if (p_ptr == NULL)
+ then error_exit(ERROR_EXIT,
+"***** patch_system::radius_in_local_xyz_direction():\n"
+" can't find containing patch!\n"
+" [local] (x,y,z)=(%g,%g,%g)\n"
+ ,
+ double(x), double(y), double(z));
+
+const patch& p = *p_ptr;
+const fp rho = p. rho_of_xyz(x,y,z);
+const fp sigma = p.sigma_of_xyz(x,y,z);
+
+
+//
+// Set up the surface interpolator to interpolate the surface radius
+// gridfn to the (rho,sigma) coordinates:
+//
+// Notes on the interpolator setup:
+// * The interpolator assumes Fortran subscripting, so we take the
+// coordinates in the order (sigma,rho) to match our C subscripting
+// in the patch system.
+// * To avoid having to set up min/max array subscripts in the parameter
+// table, we treat the patch as using 0-origin (integer) array subscripts
+// (irho - ghosted_min_irho(), isigma - ghosted_min_isigma()). However,
+// we use the usual floating-point coordinates.
+//
+
+const int N_dims = 2;
+const CCTK_REAL coord_origin[N_dims]
+ = { p.ghosted_min_sigma(), p.ghosted_min_rho() };
+const CCTK_REAL coord_delta[N_dims]
+ = { p.delta_sigma(), p.delta_rho() };
+
+const int N_interp_points = 1;
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const void* const interp_coords[N_dims]
+ = { static_cast<const void*>(&sigma), static_cast<const void*>(&rho) };
+
+const int N_input_arrays = 1;
+const CCTK_INT input_array_dims[N_dims]
+ = { p.ghosted_N_isigma(), p.ghosted_N_irho() };
+const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
+const void* const input_arrays[N_input_arrays]
+ = {
+ static_cast<const void*>(
+ p.ghosted_gridfn_data_array(ghosted_radius_gfn)
+ )
+ };
+
+const int N_output_arrays = 1;
+const CCTK_INT output_array_type_codes[N_output_arrays]
+ = { CCTK_VARIABLE_REAL };
+fp xyz_radius;
+void* const output_arrays[N_output_arrays]
+ = { static_cast<void*>(&xyz_radius) };
+
+const int status
+ = CCTK_InterpLocalUniform(N_dims,
+ surface_interp_handle_,
+ surface_interp_par_table_handle_,
+ coord_origin, coord_delta,
+ N_interp_points, interp_coords_type_code,
+ interp_coords,
+ N_input_arrays, input_array_dims,
+ input_array_type_codes,
+ input_arrays,
+ N_output_arrays, output_array_type_codes,
+ output_arrays);
+if (status < 0)
+ then error_exit(ERROR_EXIT,
+"***** patch_system::radius_in_local_xyz_direction():\n"
+" error return (status=%d) from surface interpolator!\n"
+" [local] (x,y,z)=(%g,%g,%g)\n"
+" %s patch (rho,sigma)=(%g,%g)\n"
+" (drho,dsigma)=(%g,%g)\n"
+ ,
+ status,
+ double(x), double(y), double(z),
+ p.name(), rho, sigma,
+ double(p.drho_of_rho(rho)),
+ double(p.dsigma_of_sigma(sigma))); /*NOTREACHED*/
+
+return xyz_radius;
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
// This function prints an unknown-grid gridfn in ASCII format to a
// named output file. The output format is suitable for a gnuplot
// 'splot' command. (Individual patches may be selected with the
@@ -1799,9 +1920,9 @@ fprintf(output_fp, "\n");
fp local_x, local_y, local_z;
p.xyz_of_r_rho_sigma(r,rho,sigma,
local_x,local_y,local_z);
- const fp global_x = global_x_of_local_x(local_x);
- const fp global_y = global_y_of_local_y(local_y);
- const fp global_z = global_z_of_local_z(local_z);
+ const fp global_x = origin_x() + local_x;
+ const fp global_y = origin_y() + local_y;
+ const fp global_z = origin_z() + local_z;
fprintf(output_fp,
"\t%#.10g\t%#.10g\t%#.10g",
global_x, global_y, global_z);
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 413ecb0..5f517a3 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -74,6 +74,7 @@ public:
// ***** coordinates *****
//
public:
+ #ifdef NOT_USED
// global (x,y,z) --> local (x,y,z)
fp local_x_of_global_x(fp global_x) const
{ return global_coords_.local_x_of_global_x(global_x); }
@@ -81,7 +82,9 @@ public:
{ return global_coords_.local_y_of_global_y(global_y); }
fp local_z_of_global_z(fp global_z) const
{ return global_coords_.local_z_of_global_z(global_z); }
+ #endif /* NOT_USED */
+ #ifdef NOT_USED
// local (x,y,z) --> global (x,y,z)
fp global_x_of_local_x(fp local_x) const
{ return global_coords_.global_x_of_local_x(local_x); }
@@ -89,23 +92,13 @@ public:
{ return global_coords_.global_y_of_local_y(local_y); }
fp global_z_of_local_z(fp local_z) const
{ return global_coords_.global_z_of_local_z(local_z); }
+ #endif /* NOT_USED */
// get global (x,y,z) coordinates of local origin point
fp origin_x() const { return global_coords_.origin_x(); }
fp origin_y() const { return global_coords_.origin_y(); }
fp origin_z() const { return global_coords_.origin_z(); }
- // find patch containing (ray from origin to) given local (x,y,z)
- // ... if there are multiple patches containing the position,
- // we return the one which would still contain it if patches
- // didn't overlap; if multiple patches satisfy this criterion
- // then it's arbitrary which one we return
- // ... if no patch contains the position (for a non--full-sphere
- // patch system), or the position is at the origin, then
- // we do an error_exit() and don't return to the caller
- const patch& patch_containing_local_xyz(fp x, fp y, fp z)
- const;
-
//
// ***** meta-info about the entire patch system *****
@@ -280,7 +273,34 @@ public:
//
- // ***** line/surface integrals *****
+ // ***** testing (x,y,z) point position versus a surface *****
+ //
+
+ // find patch containing (ray from origin to) given local (x,y,z)
+ // ... if there are multiple patches containing the position,
+ // we return the one which would still contain it if patches
+ // didn't overlap; if multiple patches satisfy this criterion
+ // then it's arbitrary which one we return
+ // ... if no patch contains the position (for a non--full-sphere
+ // patch system), or the position is at the origin, then
+ // we return a NULL pointer
+ const patch* patch_containing_local_xyz(fp x, fp y, fp z)
+ const;
+
+ // radius of surface in direction of an (x,y,z) point
+ // or dummy value 1.0 if point is identical to local origin
+ //
+ // FIXME:
+ // We should provide another API to compute this for a whole
+ // batch of points at once, since this would be more efficient
+ // (the interpolator overhead would be amortized over the whole batch)
+ fp radius_in_local_xyz_direction(int ghosted_radius_gfn,
+ fp x, fp y, fp z)
+ const;
+
+
+ //
+ // ***** line/surface operations *****
//
// compute the circumference of a surface in the {xy, xz, yz} plane
@@ -323,6 +343,7 @@ public:
const;
+
//
// ***** I/O *****
//
@@ -460,7 +481,9 @@ public:
int N_zones_per_right_angle,
int min_gfn_in, int max_gfn_in,
int ghosted_min_gfn_in, int ghosted_max_gfn_in,
- int interp_handle_in, int interp_par_table_handle_in,
+ int ip_interp_handle_in, int ip_interp_par_table_handle_in,
+ int surface_interp_handle_in,
+ int surface_interp_par_table_handle_in,
bool print_summary_msg_flag, bool print_detailed_msg_flag);
~patch_system();
@@ -487,31 +510,31 @@ private:
// setup (create/interlink) all ghost zones
void setup_ghost_zones__full_sphere
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_z_hemisphere
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_xy_quadrant_mirrored
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_xy_quadrant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_xz_quadrant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_xyz_octant_mirrored
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
void setup_ghost_zones__plus_xyz_octant_rotating
(int patch_overlap_width,
- int interp_handle, int interp_par_table_handle,
+ int ip_interp_handle, int ip_interp_par_table_handle,
bool print_msg_flag);
// create/interlink a pair of periodic-symmetry ghost zones
@@ -532,7 +555,7 @@ private:
void finish_interpatch_setup
(patch &px, patch &py,
int patch_overlap_width,
- int interp_handle, int interp_par_table_handle);
+ int ip_interp_handle, int ip_interp_par_table_handle);
// assert() that all ghost zones of all patches are fully setup
void assert_all_ghost_zones_fully_setup() const;
@@ -574,4 +597,8 @@ private:
// min/max m over all ghost zone points
mutable int global_min_ym_, global_max_ym_;
+
+ // info about the surface interpolator
+ // ... used only by radius_in_local_xyz_direction()
+ int surface_interp_handle_, surface_interp_par_table_handle_;
};