add support for computing 0-origin grid point numbers across all patches

and converting these <--> (patch,irho,isigma) triples


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@614 f88db872-0e4f-0410-b76b-b9085cfa78c5

2 files changed, 141 insertions, 23 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 4a3e3d1..a11a9d3 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -14,6 +14,9 @@
 // patch_system::finish_interpatch_setup
 // patch_system::assert_all_ghost_zones_fully_setup
 //
+// patch_system::patch_irho_isigma_of_gpn
+// patch_system::ghosted_patch_irho_isigma_of_gpn
+//
 // patch_system::N_patches_of_type
 // patch_system::name_of_type
 // patch_system::type_of_name
@@ -91,6 +94,8 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 	  type_(type_in),
 	  N_patches_(N_patches_of_type(type_in)),
 	  all_patches_(N_patches_),
+	  starting_gpn_(N_patches_+1),
+	  ghosted_starting_gpn_(N_patches_+1),
 	  gridfn_storage_(NULL),	// set in setup_gridfn_storage()
 	  ghosted_gridfn_storage_(NULL)	// set in setup_gridfn_storage()
 {
@@ -188,6 +193,9 @@ delete gridfn_storage_;
 // of this patch system.  This function also correctly sets
 //	N_grid_points_
 //	N_ghosted_grid_points_
+//	all_patches_[]
+//	starting_gpn_[]
+//	ghosted_starting_gpn_[]
 // This function does *NOT* create any of the ghost zones, and does
 // *NOT* set up any gridfns.
 //
@@ -246,11 +254,19 @@ ghosted_N_grid_points_ = 0;
 			   pn, pi.ctype, pi.ctype);		/*NOTREACHED*/
 		}
 
-	all_patches_[pn] = p;
+	// these record number of grid points in *previous* patches,
+	// i.e. they do *not* include the number of grid points in this patch
+	starting_gpn_[pn] = N_grid_points_;
+	ghosted_starting_gpn_[pn] = ghosted_N_grid_points_;
 
 	N_grid_points_+= p->N_grid_points();
 	ghosted_N_grid_points_ += p->ghosted_N_grid_points();
+
+	all_patches_[pn] = p;
 	}
+
+starting_gpn_        [N_patches_] = N_grid_points_;
+ghosted_starting_gpn_[N_patches_] = ghosted_N_grid_points_;
 }
 
 //******************************************************************************
@@ -318,10 +334,10 @@ ghosted_gridfn_storage_ = new fp[ghosted_N_storage];
 
 // divide up the storage array among the patches
 // and set up the storage in the individual patches themselves
-int posn = 0;
-int ghosted_posn = 0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
+	const int         posn =         starting_gpn_[pn];
+	const int ghosted_posn = ghosted_starting_gpn_[pn];
 	const struct grid_arrays::gridfn_pars gridfn_pars
 		= {
 		  min_gfn_in, max_gfn_in,
@@ -337,21 +353,7 @@ int ghosted_posn = 0;
 
 	patch& p = ith_patch(pn);
 	p.setup_gridfn_storage(gridfn_pars, ghosted_gridfn_pars);
-
-	posn += p.N_grid_points();
-	ghosted_posn += p.ghosted_N_grid_points();
 	}
-if (! (    (posn == N_grid_points())
-	&& (ghosted_posn == ghosted_N_grid_points())    ) )
-   then error_exit(PANIC_EXIT,
-"***** patch_system::setup_gridfn_storage():\n"
-"        inconsistent number of points in patch systems!\n"
-"        (this should never happen!)\n"
-"        posn=%d N_grid_points()=%d (these should be equal)\n"
-"        ghosted_posn=%d ghosted_N_grid_points()=%d (these should be equal)\n"
-,
-		   posn, N_grid_points(),
-		   ghosted_posn, ghosted_N_grid_points());	/*NOTREACHED*/
 
 CCTK_VInfo(CCTK_THORNSTRING,
 	   "   checking that storage is partitioned properly");
@@ -859,6 +861,99 @@ void patch_system::assert_all_ghost_zones_fully_setup() const
 //******************************************************************************
 
 //
+// This function decodes a 0-origin grid point number into a
+// (patch,irho,isigma) triple.
+//
+// Arguments:
+// gpn = The grid point number to decode.
+// (irho,isigma) = (out) The decoded patch coordinates.
+//
+// Results:
+// This function returns a reference to the decoded patch.  (An alternative
+// design would be to return this via a  patch*&  argument, but design here
+// avoids fiddling with pointers.)
+//
+patch& patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
+	const
+{
+assert( gpn >= 0 );
+assert( gpn < N_grid_points() );
+
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	// n.b. [pn+1] is ok since starting_gpn_[] has size N_patches()+1
+	if ((gpn >= starting_gpn_[pn])  && (gpn < starting_gpn_[pn+1]))
+	   then {
+		patch& p = ith_patch(pn);
+		const int gpn_in_patch = gpn - starting_gpn_[pn];
+		assert( gpn_in_patch >= 0 );
+		assert( gpn_in_patch < p.N_grid_points() );
+		p.irho_isigma_of_gpn(gpn_in_patch, irho,isigma);
+		return p;
+		}
+	}
+
+error_exit(PANIC_EXIT,
+"***** patch_system::patch_irho_isigma_of_gpn(gpn=%d):\n"
+"        couldn't find any patch! (this should never happen!)\n"
+"        N_grid_points()=%d\n"
+,
+	   gpn,
+	   N_grid_points());					/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
+// This function decodes a 0-origin grid point number into a *ghosted*
+// (patch,irho,isigma) triple.
+//
+// Arguments:
+// gpn = The grid point number to decode.
+// (irho,isigma) = (out) The decoded patch coordinates.
+//
+// Results:
+// This function returns a reference to the decoded patch.  (An alternative
+// design would be to return this via a  patch*&  argument, but design here
+// avoids fiddling with pointers.)
+//
+patch& patch_system::ghosted_patch_irho_isigma_of_gpn
+	(int gpn, int& irho, int& isigma)
+	const
+{
+assert( gpn >= 0 );
+assert( gpn < ghosted_N_grid_points() );
+
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	// n.b. [pn+1] is ok since ghosted_starting_gpn_[]
+	//      has size N_patches()+1
+	if (    (gpn >= ghosted_starting_gpn_[pn])
+	     && (gpn <  ghosted_starting_gpn_[pn+1]))
+	   then {
+		patch& p = ith_patch(pn);
+		const int gpn_in_patch = gpn - ghosted_starting_gpn_[pn];
+		assert( gpn_in_patch >= 0 );
+		assert( gpn_in_patch < p.ghosted_N_grid_points() );
+		p.ghosted_irho_isigma_of_gpn(gpn_in_patch, irho,isigma);
+		return p;
+		}
+	}
+
+error_exit(PANIC_EXIT,
+"***** patch_system::ghosted_patch_irho_isigma_of_gpn(gpn=%d):\n"
+"        couldn't find any patch! (this should never happen!)\n"
+"        ghosted_N_grid_points()=%d\n"
+,
+	   gpn,
+	   ghosted_N_grid_points());				/*NOTREACHED*/
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function decodes a patch system's type into N_patches.
 //
 //static
@@ -1011,10 +1106,7 @@ void patch_system::synchronize_ghost_zones(int ghosted_min_gfn_to_sync,
 		ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho);
 		if (gz.is_interpatch())
 		   then gz.synchronize(ghosted_min_gfn_to_sync,
-				       ghosted_max_gfn_to_sync,
-				       true,	// want min-par corner?
-				       true,	// want non-corner?
-				       true);	// want max-par corner?
+				       ghosted_max_gfn_to_sync);
 		}
 		}
 	}
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 8fd9541..38302bf 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -106,6 +106,26 @@ public:
 	int N_grid_points() const { return N_grid_points_; }
 	int ghosted_N_grid_points() const { return ghosted_N_grid_points_; }
 
+	// convert (patch,irho,isigma) <--> 1-D 0-origin grid point number (gpn)
+	int gpn_of_patch_irho_isigma(const patch& p, int irho, int isigma)
+		const
+		{
+		return starting_gpn_[p.patch_number()]
+		       + p.gpn_of_irho_isigma(irho,isigma);
+		}
+	int ghosted_gpn_of_patch_irho_isigma(const patch& p,
+					     int irho, int isigma)
+		const
+		{
+		return ghosted_starting_gpn_[p.patch_number()]
+		       + p.ghosted_gpn_of_irho_isigma(irho,isigma);
+		}
+	// ... n.b. we return patch as a reference via the function result;
+	//     an alternative would be to have a patch*& argument
+	patch& patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
+		const;
+	patch& ghosted_patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
+		const;
 
 	//
 	// ***** meta-info about gridfns *****
@@ -121,7 +141,7 @@ public:
 
 
 	//
-	// ***** access to gridfns as 1-D arrays *****
+	// ***** access to gridfns as 1-D arrays[gpn] *****
 	// ... n.b. this interface implicitly assumes that gridfn data
 	//     arrays are contiguous across patches; this is ensured by
 	//     setup_gridfn_storage() (called by our constructor)
@@ -321,9 +341,15 @@ private:
 	int N_patches_;
 	int N_grid_points_, ghosted_N_grid_points_;
 
-	// pointers to individual patches
+	// [pn] = --> individual patches
 	std::vector<patch *> all_patches_;
 
+	// [pn] = starting grid point number of individual patches
+	// ... arrays are actually of size N_patches_+1, the [N_patches_]
+	//     entries are == N_grid_points_ and ghosted_N_grid_points_
+	std::vector<int> starting_gpn_;
+	std::vector<int> ghosted_starting_gpn_;
+
 	// pointers to storage blocks for all gridfns
 	// ... patches point into these, but we own the storage blocks
 	fp* gridfn_storage_;