// patch.hh -- describes a coordinate/grid patch // $Id$ // // ***** how patch boundaries are handled ***** // patch - abstract base class to describe a coordinate/grid patch // // z_patch - derived class for a +/- z patch // x_patch - derived class for a +/- x patch // y_patch - derived class for a +/- y patch // // // prerequisites: // // // // "jt/util.hh" // "jt/array.hh" // "jt/linear_map.hh" // "fp.hh" // "coords.hh" // "grid.hh" // "fd_grid.hh" // //***************************************************************************** //***************************************************************************** //***************************************************************************** // // ***** how patch boundaries are handled ***** // // // Basically, we handle patch boundaries using the usual "ghost zone" // technique, interpolating values from neighboring patches as necessary. // // In more detail, we use the following interrelated types of objects // to handle patch boundaries: // // A patch_edge object represents the basic geometry of a min/max // rho/sigma side of a patch, i.e. it provides which-side-am-I predicates, // coordinate conversions between (perp,par) and (rho,sigma), etc. // Every patch has (points to) 4 patch_edge objects, one for each of // the patch's sides. // // A ghost_zone object describes a patch's ghost zone, and knows how // to fill in gridfns there based on either the patch system's symmetry // or interpolation from a neighboring patch. ghost_zone is an abstract // base class, from which we derive two classes: // * A symmetry_ghost_zone object describes a ghost zone which is a // (discrete) symmetry of spacetime, either mirror-image or periodic. // Such an object knows how to fill in ghost-zone gridfn data from // the "other side" of the symmetry. // * An interpatch_ghost_zone object describes a ghost zone which // overlaps another patch. Such an object knows how to get ghost // zone gridfn data from the other patch. More accurately, it gets // the data by asking (calling) the appropriate one of the other // patch's patch_frontier objects. // Every patch has (points to) 4 ghost_zone objects, one for each of // the patch's sides. // // A patch_frontier object represents the part of a patch *inside* // its nominal grid from which data is interpolated to compute gridfns // in another patch's ghost zone. A patch_frontier object stores // information about the interpolator (which interpolator to use, any // options for it) and any cached data for the interpolator. Every // patch has (points to) between 0 and 4 patch_frontier objects, // one for each of the patch's sides where there's an adjacent patch // and thus from which interpolation is done. // // For example, suppose we have two patches p and q with a common // angular boundary. Then the desired network of pointers looks like // this (omitting the patch_edge objects for simplicity): // // +-----+ +-----+ // | | <--> p.interpatch_ghost_zone --> q.patch_frontier <--> | | // | p | | q | // | | <--> p.patch_frontier <-- q.interpatch_ghost_zone <--> | | // +-----+ +-----+ // // Because of the mutual pointers, we can't easily construct (say) // p's interpatch_ghost_zone until after q itself has been constructed, // and vice versa. Moreover, the patch_frontier:: constructor // needs the adjacent-side ghost_zone objects to already exist, // and it needs to know the iperp range of the frontier, which can // only be computed from the adjacent-patch interpatch_ghost_zone // object. // // The solution adopted here is to use a 4-phase algorithm, ultimately // driver by the patch_system constructor: // * The patch constructors themselves construct the patch_edge objects // and links them to/from the patches. // * The patch_system constructor calls the appropriate functions // patch::setup_mirror_symmetry_ghost_zone() // patch::setup_periodic_symmetry_ghost_zone() // patch::setup_interpatch_ghost_zone() // to construct the ghost_zone objects and link them to/from the // patches. // * The patch_system constructor calls the function // patch::setup_patch_frontier() // to construct the patch_frontier objects and link them to/from the // patches. // * The patch_system constructor calls the function // interpatch_ghost_zone::set_other_frontier() // to link the patch_frontier objects from the interpatch_ghost_zone // objects. // //***************************************************************************** //***************************************************************************** //***************************************************************************** // // patch - abstract base class to describe a generic coordinate/grid patch // // // There are 3 types of patches, z, x, and y. Each type uses two of // (mu,nu,phi) as its angular coordinates (rho,sigma); the remaining // "unused" one of (mu,nu,phi) is tau. // // z patch ==> (rho,sigma) = (mu,nu) tau = phi // x patch ==> (rho,sigma) = (nu,phi) tau = mu // y patch ==> (rho,sigma) = (mu,phi) tau = nu // // forward declarations class patch_edge; class ghost_zone; class symmetry_ghost_zone; class interpatch_ghost_zone; class patch_frontier; class patch_system; // // const qualifiers refer to the gridfn values // class patch : public fd_grid { public: // // ***** patch system, type, and coordinate metadata ***** // // to which patch system do we belong? patch_system& my_patch_system() const { return my_patch_system_; } // each patch has a 0-origin small-integer patch number, // usually denoted pn int patch_number() const { return patch_number_; } // each patch has a human-readable patch name for debugging etc // ... this should be unique among all patches const char *name() const { return name_; } // typically "+z" etc // are we a +[xyz] or -[xyz] patch? bool is_plus() const { return is_plus_; } // ... values for the is_plus_in constructor argument static const bool patch_is_plus = true; static const bool patch_is_minus = false; // are we a (+/-) x or y or z patch? // ... n.b. type is `char' because this is handy for both // switch() and human-readable printing char ctype() const { return ctype_; } // 'z' or 'x' or 'y' // (rho,sigma,tau) coordinates as singleton coordinate sets local_coords::coords_set coords_set_rho() const { return coords_set_rho_; } local_coords::coords_set coords_set_sigma() const { return coords_set_sigma_; } local_coords::coords_set coords_set_tau() const { return coords_set_tau_; } // {rho,sigma} coordinate set local_coords::coords_set coords_set_rho_sigma() const { return coords_set_rho() | coords_set_sigma(); } // (rho,sigma) coordinates as human-readable character strings // (for labelling output files etc) virtual const char *name_of_rho() const = 0; virtual const char *name_of_sigma() const = 0; // // ***** (rho,sigma,tau) coordinate conversions ***** // // convert (rho,sigma) --> tau virtual fp tau_of_rho_sigma(fp rho, fp sigma) const = 0; // convert (rho,sigma) --> (mu,nu,phi) virtual fp mu_of_rho_sigma(fp rho, fp sigma) const = 0; virtual fp nu_of_rho_sigma(fp rho, fp sigma) const = 0; virtual fp phi_of_rho_sigma(fp rho, fp sigma) const = 0; // convert (rho,sigma) <--> usual polar spherical (theta,phi) virtual void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi) const = 0; virtual void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma) const = 0; // convert (r,rho,sigma) <--> local (x,y,z) virtual void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z) const = 0; virtual fp rho_of_xyz(fp x, fp y, fp z) const = 0; virtual fp sigma_of_xyz(fp x, fp y, fp z) const = 0; // plotting coordinates (dpx,dpy) // ... character string describing how (dpx,dpy) are // defined in terms of (mu,nu,phi), eg "90 - dphi" // (used for labelling output files) virtual const char *name_of_dpx() const = 0; virtual const char *name_of_dpy() const = 0; // ... (irho,isimga) --> (px,py) virtual fp dpx_of_irho_isigma(int irho, int isigma) const = 0; virtual fp dpy_of_irho_isigma(int irho, int isigma) const = 0; // print a gridfn (via C stdio) in ASCII format // output is assumed to already be open // output format is suitable for gnuplot 'splot' void print_gridfn(int gfn, bool want_ghost_zones, FILE *output_fp = stdout) const; // // ***** patch edges **** // const patch_edge& min_rho_patch_edge() const { return min_rho_patch_edge_; } const patch_edge& max_rho_patch_edge() const { return max_rho_patch_edge_; } const patch_edge& min_sigma_patch_edge() const { return min_sigma_patch_edge_; } const patch_edge& max_sigma_patch_edge() const { return max_sigma_patch_edge_; } const patch_edge& minmax_ang_patch_edge(bool want_min, bool want_rho) const { return want_min ? (want_rho ? min_rho_patch_edge() : min_sigma_patch_edge()) : (want_rho ? max_rho_patch_edge() : max_sigma_patch_edge()); } // find which patch edge is adjacent to neighboring patch q, // or error_exit() if it's not actually a neighboring patch // ... computation done using only (rho,sigma) coordinate sets // and min/max dang bounds ==> ok to use in setting up ghost zones // ... N_overlap_points = number of grid points (grid spacings // in the perpendicular direction) these patches overlap // ... if this is nonzero, then these patches must have // the *same* grid spacing in the perpendicular direction // ... e.g. delta_dang = 5, this patch max_dang = 50, // other patch min_dang = 40 ==> N_overlap_points = 2 const patch_edge& edge_adjacent_to_patch(const patch& q, int N_overlap_points = 0) const; // // ***** ghost zones ***** // ghost_zone& min_rho_ghost_zone() const { return *min_rho_ghost_zone_; } ghost_zone& max_rho_ghost_zone() const { return *max_rho_ghost_zone_; } ghost_zone& min_sigma_ghost_zone() const { return *min_sigma_ghost_zone_; } ghost_zone& max_sigma_ghost_zone() const { return *max_sigma_ghost_zone_; } ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho) const { return want_min ? (want_rho ? min_rho_ghost_zone() : min_sigma_ghost_zone()) : (want_rho ? max_rho_ghost_zone() : max_sigma_ghost_zone()); } ghost_zone& ghost_zone_on_edge(const patch_edge &e) const; private: // helper function for setup_*_ghost_zone(): // assert() that ghost zone pointer on specified edge is NULL // (i.e. that we haven't already setup this ghost zone), // then assign new value to it void set_ghost_zone(ghost_zone* gzp, const patch_edge& edge); public: // get *interpatch* ghost zone on specified edge // ... verifies that ghost zone is indeed interpatch interpatch_ghost_zone& interpatch_ghost_zone_on_edge (const patch_edge &e) const; // // ***** set up ghost zone and frontier subobjects // // assert() that this ghost zone hasn't been set up yet, // then set it up as mirror-symmetry // ... return reference to newly-set-up ghost zone object symmetry_ghost_zone& setup_mirror_symmetry_ghost_zone (const patch_edge& edge); // assert() that this ghost zone hasn't been set up yet, // then set it up as periodic-symmetry // ... return reference to newly-set-up ghost zone object symmetry_ghost_zone& setup_periodic_symmetry_ghost_zone (const patch_edge& my_edge, const patch_edge& symmetry_edge, bool ipar_map_sign); // one of cpm_map::{plus,minus}_map // assert() that this ghost zone hasn't been set up yet, // then set it up as interpatch // ... but don't link it to the other patch's frontier // (which doesn't exist yet) // ... return reference to newly-set-up ghost zone object interpatch_ghost_zone& setup_interpatch_ghost_zone (const patch_edge& my_edge, const patch_edge& other_edge, int N_overlap_points); // assert() that this frontier hasn't been set up yet, then set it up // ... return reference to new frontier patch_frontier& setup_patch_frontier(const patch_edge& my_edge_in, int min_iperp_in, int max_iperp_in, int interpolator_order); // assert() that all ghost zones (and frontiers, where applicable) // are fully setup void assert_all_ghost_zones_fully_setup() const; protected: // ... used only from derived classes // ... doesn't set up ghost zone or frontier info, since this // depends on knowing our neighbouring patches, which may // not exist yet patch(patch_system &my_patch_system_in, int patch_number_in, const char *name_in, bool is_plus_in, char ctype_in, local_coords::coords_set coords_set_rho_in, local_coords::coords_set coords_set_sigma_in, local_coords::coords_set coords_set_tau_in, const grid_arrays::grid_array_pars& grid_array_pars_in, const grid::grid_pars& grid_pars_in, int min_gfn_in, int max_gfn_in); public: // destructor must be virtual to allow destruction // of derived classes via ptr/ref to this class virtual ~patch(); private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them patch(const patch& rhs); patch& operator=(const patch& rhs); private: // // ***** data members ***** // // type/coordinate metadata patch_system &my_patch_system_; const int patch_number_; const char *name_; const bool is_plus_; const char ctype_; const local_coords::coords_set coords_set_rho_, coords_set_sigma_, coords_set_tau_; // edges // ... references are initialized in constructor const patch_edge& min_rho_patch_edge_, & max_rho_patch_edge_, & min_sigma_patch_edge_, & max_sigma_patch_edge_; // ghost zones // ... pointers are initialized by // setup_mirror_symmetry_ghost_zone() // setup_periodic_symmetry_ghost_zone() // setup_interpatch_ghost_zone() ghost_zone *min_rho_ghost_zone_, *max_rho_ghost_zone_, *min_sigma_ghost_zone_, *max_sigma_ghost_zone_; // frontiers (NULL pointers if no frontier for this ghost zone) // ... pointers are initialized by setup_interpatch_ghost_zone() patch_frontier *min_rho_patch_frontier_, *max_rho_patch_frontier_, *min_sigma_patch_frontier_, *max_sigma_patch_frontier_; }; //***************************************************************************** //***************************************************************************** //***************************************************************************** // // This class describes a +/- z patch. It doesn't define any new // functions not already present in class patch ; it "just" defines // non-virtual versions of all the pure virtual functions defined there; // // z patch ==> (rho,sigma) = (mu,nu) tau = phi // class z_patch : public patch { public: // human-readable names of (rho,sigma) const char *name_of_rho() const { return "mu"; } const char *name_of_sigma() const { return "nu"; } // convert (rho,sigma) --> tau fp tau_of_rho_sigma(fp rho, fp sigma) const { return local_coords::phi_of_mu_nu(rho, sigma); } // convert (rho,sigma) --> (mu,nu,phi) fp mu_of_rho_sigma(fp rho, fp sigma) const { return rho; } fp nu_of_rho_sigma(fp rho, fp sigma) const { return sigma; } fp phi_of_rho_sigma(fp rho, fp sigma) const { return local_coords::phi_of_mu_nu(rho, sigma); } // convert (rho,sigma) <--> usual polar spherical (theta,phi) void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi) const { local_coords::theta_phi_of_mu_nu(rho, sigma, ps_theta, ps_phi); } void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma) const { local_coords::mu_nu_of_theta_phi(ps_theta, ps_phi, rho, sigma); } // convert (r,rho,sigma) <--> (x,y,z) void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z) 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); } fp sigma_of_xyz(fp x, fp y, fp z) const { return local_coords::nu_of_xz(x, z); } // plotting coordinates (px,py) // ... character string describing how (dpx,dpy) are // defined in terms of (mu,nu,phi), eg "90 - dphi" // (used for labelling output files) const char *name_of_dpx() const { return "dnu"; } const char *name_of_dpy() const { return is_plus() ? "dmu" : "dmu - 180"; } // ... (irho,isimga) --> (px,py) fp dpx_of_irho_isigma(int irho, int isigma) const { return jtutil::degrees_of_radians(sigma_of_isigma(isigma)); } fp dpy_of_irho_isigma(int irho, int sigma) const { const fp dmu = jtutil::degrees_of_radians(rho_of_irho(irho)); return is_plus() ? dmu : dmu - 180.0; } // constructor, destructor z_patch(patch_system &my_patch_system_in, int patch_number_in, const char *name_in, bool is_plus_in, const grid_arrays::grid_array_pars& grid_array_pars_in, const grid::grid_pars& grid_pars_in, int min_gfn_in, int max_gfn_in); ~z_patch() { } private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them z_patch(const z_patch& rhs); z_patch& operator=(const z_patch& rhs); }; //***************************************************************************** // // This class describes a +/- x patch. It doesn't define any new // functions not already present in class patch ; it "just" defines // non-virtual versions of all the pure virtual functions defined there; // // x patch ==> (rho,sigma) = (nu,phi) tau = mu // class x_patch : public patch { public: // human-readable names of (rho,sigma) const char *name_of_rho() const { return "nu"; } const char *name_of_sigma() const { return "phi"; } // convert (rho,sigma) --> tau fp tau_of_rho_sigma(fp rho, fp sigma) const { return local_coords::mu_of_nu_phi(rho, sigma); } // convert (rho,sigma) --> (mu,nu,phi) fp nu_of_rho_sigma(fp rho, fp sigma) const { return rho; } fp phi_of_rho_sigma(fp rho, fp sigma) const { return sigma; } fp mu_of_rho_sigma(fp rho, fp sigma) const { return local_coords::mu_of_nu_phi(rho, sigma); } // convert (rho,sigma) <--> usual polar spherical (theta,phi) void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi) const { local_coords::theta_phi_of_nu_phi(rho, sigma, ps_theta, ps_phi); } void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma) const { local_coords::nu_phi_of_theta_phi(ps_theta, ps_phi, rho, sigma); } // convert (r,rho,sigma) <--> (x,y,z) void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z) 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); } fp sigma_of_xyz(fp x, fp y, fp z) const { return local_coords::phi_of_xy(x, y); } // plotting coordinates (px,py) // ... character string describing how (dpx,dpy) are // defined in terms of (mu,nu,phi), eg "90 - dphi" // (used for labelling output files) const char *name_of_dpx() const { return "dnu"; } const char *name_of_dpy() const { return is_plus() ? "dphi" : "dphi - 180"; } // ... (irho,isimga) --> (px,py) fp dpx_of_irho_isigma(int irho, int isigma) const { return jtutil::degrees_of_radians(rho_of_irho(irho)); } fp dpy_of_irho_isigma(int irho, int isigma) const { const fp dphi = jtutil::degrees_of_radians(sigma_of_isigma(isigma)); return is_plus() ? dphi : dphi - 180.0; } // constructor, destructor x_patch(patch_system &my_patch_system_in, int patch_number_in, const char *name_in, bool is_plus_in, const grid_arrays::grid_array_pars& grid_array_pars_in, const grid::grid_pars& grid_pars_in, int min_gfn_in, int max_gfn_in); ~x_patch() { } private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them x_patch(const x_patch& rhs); x_patch& operator=(const x_patch& rhs); }; //***************************************************************************** // // This class describes a +/- y patch. It doesn't define any new // functions not already present in class patch ; it "just" defines // non-virtual versions of all the pure virtual functions defined there; // // y patch ==> (rho,sigma) = (mu,phi) tau = nu // class y_patch : public patch { public: // human-readable names of (rho,sigma) const char *name_of_rho() const { return "mu"; } const char *name_of_sigma() const { return "phi"; } // convert (rho,sigma) --> tau fp tau_of_rho_sigma(fp rho, fp sigma) const { return local_coords::nu_of_mu_phi(rho, sigma); } // convert (rho,sigma) --> (mu,nu,phi) fp mu_of_rho_sigma(fp rho, fp sigma) const { return rho; } fp phi_of_rho_sigma(fp rho, fp sigma) const { return sigma; } fp nu_of_rho_sigma(fp rho, fp sigma) const { return local_coords::nu_of_mu_phi(rho, sigma); } // convert (rho,sigma) <--> usual polar spherical (theta,phi) void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi) const { local_coords::theta_phi_of_mu_phi(rho, sigma, ps_theta, ps_phi); } void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma) const { local_coords::mu_phi_of_theta_phi(ps_theta, ps_phi, rho, sigma); } // convert (r,rho,sigma) <--> (x,y,z) void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z) 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); } fp sigma_of_xyz(fp x, fp y, fp z) const { return local_coords::phi_of_xy(x, y); } // plotting coordinates (px,py) // ... character string describing how (dpx,dpy) are // defined in terms of (mu,nu,phi), eg "90 - dphi" // (used for labelling output files) const char *name_of_dpx() const { return is_plus() ? "90 - dphi" : "90 + dphi"; } const char *name_of_dpy() const { return "dmu"; } // ... (rho,simga) --> (px,py) fp dpx_of_irho_isigma(int irho, int isigma) const { const fp dphi = jtutil::degrees_of_radians(sigma_of_isigma(isigma)); return is_plus() ? 90.0 - dphi : 90.0 + dphi; } fp dpy_of_irho_isigma(int irho, int isigma) const { return jtutil::degrees_of_radians(rho_of_irho(irho)); } // constructor, destructor y_patch(patch_system &my_patch_system_in, int patch_number_in, const char *name_in, bool is_plus_in, const grid_arrays::grid_array_pars& grid_array_pars_in, const grid::grid_pars& grid_pars_in, int min_gfn_in, int max_gfn_in); ~y_patch() { } private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them y_patch(const y_patch& rhs); y_patch& operator=(const y_patch& rhs); };