diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-03-17 12:06:49 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-03-17 12:06:49 +0000 |
commit | 8d274d170dde9d858e0a13a43762d3976ea4e3b1 (patch) | |
tree | d394b0725cf7af6883c57b8bde59f0ffef479d6c /src/patch | |
parent | 498e26dd4ee819db7a0c4930e73624af936b7d47 (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.hh | 17 | ||||
-rw-r--r-- | src/patch/ghost_zone.cc | 2 | ||||
-rw-r--r-- | src/patch/ghost_zone.hh | 2 | ||||
-rw-r--r-- | src/patch/patch.hh | 12 | ||||
-rw-r--r-- | src/patch/patch_interp.hh | 11 | ||||
-rw-r--r-- | src/patch/patch_system.cc | 369 | ||||
-rw-r--r-- | src/patch/patch_system.hh | 69 |
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_; }; |