merge in changes from laptop:

- fix nasty bug in interpatch_ghost_zone::min_ipar()
  and interpatch_ghost_zone::max_ipar()
  where we used the corners when we shouldn't
- fix similar bug in interpatch_ghost_zone::finish_setup()
  where we used *our* adjacent ghost zone types to decide what to tell
  the other patch's patch_interp object about whether or not to use
  its min/max par corners -- this should be the *other* patch's adjacent
  ghost zones
- rename that to min/max par ghost zones
- rename some other functions
- finish implementing patch_system::synchronize_Jacobian()
  to properly take into account the 3-phase algorithm of
  patch_system::synchronize()


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

9 files changed, 585 insertions, 196 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 1b4ba11..3ace6e3 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -2,6 +2,9 @@
 // $Id$
 
 //
+// ghost_zone::cast_to_symmetry_ghost_zone
+// ghost_zone::cast_to_interpatch_ghost_zone
+//
 // symmetry_ghost_zone::symmetry_ghost_zone (mirror symmetry)
 // symmetry_ghost_zone::symmetry_ghost_zone (periodic BC)
 // symmetry_ghost_zone::~symmetry_ghost_zone
@@ -47,12 +50,49 @@ using jtutil::error_exit;
 //******************************************************************************
 
 //
+// These functions verify (assert()) that a ghost zone is indeed of
+// the specified type, then static_cast to the appropriate derived class.
+//
+
+const symmetry_ghost_zone& ghost_zone::cast_to_symmetry_ghost_zone()
+	const
+{
+assert( is_symmetry() );
+return static_cast<const symmetry_ghost_zone &>(*this);
+}
+
+symmetry_ghost_zone& ghost_zone::cast_to_symmetry_ghost_zone()
+{
+assert( is_symmetry() );
+return static_cast<symmetry_ghost_zone &>(*this);
+}
+
+//**************************************
+
+const interpatch_ghost_zone& ghost_zone::cast_to_interpatch_ghost_zone()
+	const
+{
+assert( is_interpatch() );
+return static_cast<const interpatch_ghost_zone &>(*this);
+}
+
+interpatch_ghost_zone& ghost_zone::cast_to_interpatch_ghost_zone()
+{
+assert( is_interpatch() );
+return static_cast<interpatch_ghost_zone &>(*this);
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function constructs a mirror-symmetry ghost zone object
 //
 symmetry_ghost_zone::symmetry_ghost_zone(const patch_edge& my_edge_in)
-	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
-	  symmetry_patch_(my_edge_in.my_patch()),
-	  symmetry_edge_(my_edge_in)
+	: ghost_zone(my_edge_in,
+		     my_edge_in,	// other edge == my edge
+		     ghost_zone_is_symmetry)
 {
 // iperp_map: i --> (i of ghost zone) - i
 iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
@@ -68,36 +108,35 @@ ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar());
 // This function constructs a periodic-symmetry ghost zone object.
 //
 symmetry_ghost_zone::symmetry_ghost_zone
-	(const patch_edge& my_edge_in, const patch_edge& symmetry_edge_in,
-	 int my_edge_sample_ipar,      int symmetry_edge_sample_ipar,
+	(const patch_edge& my_edge_in, const patch_edge& other_edge_in,
+	 int my_edge_sample_ipar,      int other_edge_sample_ipar,
 	 bool ipar_map_is_plus)
-	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
-	  symmetry_patch_(symmetry_edge_in.my_patch()),
-	  symmetry_edge_(symmetry_edge_in)
+	: ghost_zone(my_edge_in,
+		     other_edge_in,
+		     ghost_zone_is_symmetry)
 {
 //
 // perpendicular map
 //
-const fp fp_my_period_plane_iperp = my_edge().fp_grid_outer_iperp();
-const fp fp_symmetry_period_plane_iperp = symmetry_edge().fp_grid_outer_iperp();
+const fp fp_my_period_plane_iperp    = my_edge()   .fp_grid_outer_iperp();
+const fp fp_other_period_plane_iperp = other_edge().fp_grid_outer_iperp();
 
 // iperp mapping must be outside --> inside
 // i.e. if both edges have iperp as the same min/max "direction",
 //	then the mapping is  iperp increasing --> iperp decreasing
 //      (i.e. the map's sign is -1)
 const bool is_iperp_map_plus
-	= ! (my_edge().is_min() == symmetry_edge().is_min());
+	= ! (my_edge().is_min() == other_edge().is_min());
 iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 				     fp_my_period_plane_iperp,
-				     fp_symmetry_period_plane_iperp,
+				     fp_other_period_plane_iperp,
 				     is_iperp_map_plus);
 
 //
 // parallel map
 //
 ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar(),
-				    my_edge_sample_ipar,
-				    symmetry_edge_sample_ipar,
+				    my_edge_sample_ipar, other_edge_sample_ipar,
 				    ipar_map_is_plus);
 }
 
@@ -122,7 +161,7 @@ delete iperp_map_;
 //
 void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 				      bool want_corners = true,
-				      bool want_non_corner = true)
+				      bool want_noncorner = true)
 {
 	for (int gfn = ghosted_min_gfn ; gfn <= ghosted_max_gfn ; ++gfn)
 	{
@@ -131,17 +170,17 @@ void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
 	{
 	// do we want to do this point?
-	if (! my_edge().ipar_is_in_selected_part(want_corners, want_non_corner,
+	if (! my_edge().ipar_is_in_selected_part(want_corners, want_noncorner,
 						 ipar) )
 	   then continue;				// *** LOOP CONTROL ***
 
 	const int sym_iperp = iperp_map_of_iperp(iperp);
 	const int sym_ipar  = ipar_map_of_ipar  (ipar );
-	const int sym_irho = symmetry_edge()
+	const int sym_irho = other_edge()
 			     .irho_of_iperp_ipar  (sym_iperp,sym_ipar);
-	const int sym_isigma = symmetry_edge()
+	const int sym_isigma = other_edge()
 			       .isigma_of_iperp_ipar(sym_iperp,sym_ipar);
-	const fp sym_gridfn = symmetry_patch()
+	const fp sym_gridfn = other_patch()
 			      .ghosted_gridfn(gfn, sym_irho,sym_isigma);
 
 	const int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
@@ -162,9 +201,9 @@ void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
 					     const patch_edge& other_edge_in,
 					     int N_overlap_points)
-	: ghost_zone(my_edge_in, ghost_zone_is_interpatch),
-	  other_patch_(other_edge_in.my_patch()),
-	  other_edge_(other_edge_in),
+	: ghost_zone(my_edge_in,
+		     other_edge_in,
+		     ghost_zone_is_interpatch),
 	  // remaining pointers are all set up properly by finish_setup()
 	  other_patch_interp_(NULL),
 	  other_iperp_(NULL),
@@ -342,7 +381,7 @@ delete other_patch_interp_;
 // a given iperp, taking into account how we treat the corners
 // (cf. the example in the header comments in "ghost_zone.hh"):
 //
-// If an adjacent ghost zone is mirror-symmetry,
+// If an adjacent ghost zone is symmetry,
 //    we do not include that corner;
 // If an adjacent ghost zone is interpatch,
 //    we include up to the diagonal line, and if we are a rho ghost zone,
@@ -358,7 +397,7 @@ delete other_patch_interp_;
 int interpatch_ghost_zone::min_ipar(int iperp) const
 {
 return min_par_adjacent_ghost_zone().is_symmetry()
-       ? my_edge().min_ipar_with_corners()
+       ? my_edge().min_ipar_without_corners()
        : my_edge().min_ipar_without_corners()
 	 - iabs(iperp - my_edge().nominal_grid_outer_iperp())
 	 + (is_rho() ? 0 : 1);
@@ -367,7 +406,7 @@ return min_par_adjacent_ghost_zone().is_symmetry()
 int interpatch_ghost_zone::max_ipar(int iperp) const
 {
 return max_par_adjacent_ghost_zone().is_symmetry()
-       ? my_edge().max_ipar_with_corners()
+       ? my_edge().max_ipar_without_corners()
        : my_edge().max_ipar_without_corners()
 	 + iabs(iperp - my_edge().nominal_grid_outer_iperp())
 	 - (is_rho() ? 0 : 1);
@@ -393,9 +432,8 @@ max_other_iperp_ = std::max(other_iperp(min_iperp()), other_iperp(max_iperp()));
 
 //
 // set up arrays giving actual [min,max] ipar that we'll use
-// at each other_iperp
-// ... we will pass these arrays by reference to the other patch's
-//     patch_interp object, with ipar being parindex there
+// at each other_iperp (later on we will pass these arrays to the
+// other patch's  patch_interp  object, with ipar being parindex there
 //
 min_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
 max_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
@@ -447,24 +485,29 @@ interp_result_buffer_
 
 //
 // construct the patch_interp object to interpolate from the *other* patch
-// ... tell it to use corners if and only if the adjacent ghost zones
-//     are symmetry; cf  min_ipar()  and  max_ipar()  functions above,
-//     and header comments in this file
-//
-const bool ok_to_use_min_par_corner
-	= min_par_adjacent_ghost_zone().is_symmetry()
+// ... the patch_interp should use gridfn data from it's (the other patch's)
+//     min/max par ghost zones if those (adjacent) adjacent ghost zones
+//     are symmetry, but not if they're interpatch,
+//     cf the header comments in "ghost_zone.hh"
+//
+const ghost_zone& other_ghost_zone = other_patch()
+				     .ghost_zone_on_edge(other_edge());
+const bool ok_to_use_min_par_ghost_zone
+	= other_ghost_zone.min_par_adjacent_ghost_zone()
+			  .is_symmetry()
 	  ? true
 	  : false;
-const bool ok_to_use_max_par_corner
-	= max_par_adjacent_ghost_zone().is_symmetry()
+const bool ok_to_use_max_par_ghost_zone
+	= other_ghost_zone.max_par_adjacent_ghost_zone()
+			  .is_symmetry()
 	  ? true
 	  : false;
 other_patch_interp_ = new patch_interp(other_edge(),
 				       min_other_iperp_, max_other_iperp_,
 				       *min_ipar_used_, *max_ipar_used_,
 				       *other_par_,
-				       ok_to_use_min_par_corner,
-				       ok_to_use_max_par_corner,
+				       ok_to_use_min_par_ghost_zone,
+				       ok_to_use_max_par_ghost_zone,
 				       interp_handle, interp_par_table_handle);
 }
 
@@ -487,7 +530,7 @@ assert( other_patch()
 	.ghost_zone_on_edge(other_edge())
 	.is_interpatch() );
 assert( my_patch() == other_patch()
-		      .interpatch_ghost_zone_on_edge(other_edge())
+		      .ghost_zone_on_edge(other_edge())
 		      .other_patch() );
 }
 
@@ -505,17 +548,17 @@ assert( my_patch() == other_patch()
 void interpatch_ghost_zone::synchronize
 	(int ghosted_min_gfn, int ghosted_max_gfn,
 	 bool want_corners = true,
-	 bool want_non_corner = true)
+	 bool want_noncorner = true)
 {
 // make sure the caller wants the entire ghost zone
-if (! (want_corners && want_non_corner))
+if (! (want_corners && want_noncorner))
    then error_exit(ERROR_EXIT,
 "***** interpatch_ghost_zone::synchronize():\n"
 "        we only support operating on the *entire* ghost zone,\n"
 "        but we were passed flags specifying a proper subset!\n"
-"        want_corners=(int)%d want_non_corner=(int)%d\n"
+"        want_corners=(int)%d want_noncorner=(int)%d\n"
 ,
-		   want_corners, want_non_corner);		/*NOTREACHED*/
+		   want_corners, want_noncorner);		/*NOTREACHED*/
 
 // do the interpolation into our result buffer
 other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
@@ -555,18 +598,18 @@ other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
 void interpatch_ghost_zone::compute_Jacobian
 	(int ghosted_min_gfn, int ghosted_max_gfn,
 	 bool want_corners = true,
-	 bool want_non_corner = true)
+	 bool want_noncorner = true)
 	const
 {
 // make sure the caller wants the entire ghost zone
-if (! (want_corners && want_non_corner))
+if (! (want_corners && want_noncorner))
    then error_exit(ERROR_EXIT,
 "***** interpatch_ghost_zone::compute_Jacobian():\n"
 "        we only support operating on the *entire* ghost zone,\n"
 "        but we were passed flags specifying a proper subset!\n"
-"        want_corners=(int)%d want_non_corner=(int)%d\n"
+"        want_corners=(int)%d want_noncorner=(int)%d\n"
 ,
-		   want_corners, want_non_corner);		/*NOTREACHED*/
+		   want_corners, want_noncorner);		/*NOTREACHED*/
 
 assert(other_patch_interp_ != NULL);
 other_patch_interp_->verify_Jacobian_sparsity_pattern_ok();
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 6033d32..cb9e328 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -271,10 +271,7 @@ public:
 	// The API in the remaining functions implicitly assumes that
 	// the Jacobian is independent of  ghosted_gfn , and also that
 	// the structure of the Jacobian is such that the set of y points
-	// (with nonzero Jacobian values) in a single row of the Jacobian
-	// matrix (i.e. the set of points on which a single ghost-zone point
-	// depends),
-	// - lies entirely within a single y patch
+	// on which a single ghost-zone point depends,
 	// - has a single yiperp value (depending on our iperp, of course)
 	// - have a contiguous interval of yipar (depending on our iperp
 	//   and ipar, of course), whose size is
@@ -286,10 +283,6 @@ public:
 	//   of our iperp and ipar
 	//
 
-	// to which patch/edge do the y points in this Jacobian row belong?
-	virtual patch&            Jacobian_y_patch() const = 0;
-	virtual const patch_edge& Jacobian_y_edge()  const = 0;
-
 	// what is the [min,max] range of m for this ghost zone?
 	virtual int Jacobian_min_y_ipar_m() const = 0;
 	virtual int Jacobian_max_y_ipar_m() const = 0;
@@ -316,8 +309,12 @@ public:
 	//
 
 	// to which patch/edge do we belong?
-	patch& my_patch() const { return my_patch_; }
-	const patch_edge& my_edge() const { return my_edge_; }
+	patch&            my_patch() const { return my_patch_; }
+	const patch_edge& my_edge()  const { return my_edge_; }
+
+	// from which patch/edge do we get data?
+	patch&            other_patch() const { return other_patch_; }
+	const patch_edge& other_edge()  const { return other_edge_; }
 
 	// what type of ghost zone are we?
 	bool is_interpatch() const { return  is_interpatch_; }
@@ -340,10 +337,8 @@ public:
 		}
 
 	// inner/outer iperp of the ghost zone wrt our patch
-	int inner_iperp() const
-		{ return is_min() ? max_iperp() : min_iperp(); }
-	int outer_iperp() const
-		{ return is_min() ? min_iperp() : max_iperp(); }
+	int inner_iperp() const { return is_min() ? max_iperp() : min_iperp(); }
+	int outer_iperp() const { return is_min() ? min_iperp() : max_iperp(); }
 
 	// extreme min/max ipar that might possibly be part of this ghost zone
 	// (derived classes may actually use a subset of this)
@@ -380,6 +375,17 @@ public:
 		}
 
 	//
+	// ***** safely cast to derived classes *****
+	//
+
+	// assert that gz is of specified type,
+	// then static_cast to derive type
+	const symmetry_ghost_zone& cast_to_symmetry_ghost_zone() const;
+	      symmetry_ghost_zone& cast_to_symmetry_ghost_zone();
+	const interpatch_ghost_zone& cast_to_interpatch_ghost_zone() const;
+	      interpatch_ghost_zone& cast_to_interpatch_ghost_zone();
+
+	//
 	// ***** constructor, finish setup, destructor *****
 	//
 protected:
@@ -391,9 +397,12 @@ protected:
 	// ... only used in implementing our derived classes;
 	//     the rest of the world constructs our derived classes instead
 	ghost_zone(const patch_edge& my_edge_in,
+		   const patch_edge& other_edge_in,
 		   bool is_interpatch_in)
 		: my_patch_(my_edge_in.my_patch()),
 		  my_edge_(my_edge_in),
+		  other_patch_(other_edge_in.my_patch()),
+		  other_edge_(other_edge_in),
 		  is_interpatch_(is_interpatch_in)
 		{ }
 public:
@@ -419,6 +428,8 @@ private:
 private:
 	patch& my_patch_;
 	const patch_edge& my_edge_;
+	patch& other_patch_;
+	const patch_edge& other_edge_;
 	const bool is_interpatch_;
 	};
 
@@ -503,10 +514,6 @@ public:
 		const
 		{ }
 
-	// to which patch/edge do the y points in this Jacobian row belong?
-	patch&            Jacobian_y_patch() const { return symmetry_patch_; }
-	const patch_edge& Jacobian_y_edge()  const { return symmetry_edge_; }
-
 	// what is the [min,max] range of m for this ghost zone?
 	int Jacobian_min_y_ipar_m() const { return 0; }
 	int Jacobian_max_y_ipar_m() const { return 0; }
@@ -535,8 +542,6 @@ public:
 	//
 
 	// symmetry-map coordinates
-	patch& symmetry_patch() const { return symmetry_patch_; }
-	const patch_edge& symmetry_edge() const { return symmetry_edge_; }
 	int iperp_map_of_iperp(int iperp) const
 		{ return iperp_map_->map(iperp); }
 	int ipar_map_of_ipar(int ipar) const
@@ -556,15 +561,14 @@ public:
 	// ***** constructors, destructor *****
 	//
 public:
-
 	// constructor for mirror-symmetry ghost zone
 	symmetry_ghost_zone(const patch_edge& my_edge_in);
 
 	// constructor for periodic-symmetry ghost zone
 	// ... ipar mapping specified by giving sample point and mapping sign
 	symmetry_ghost_zone
-	    (const patch_edge& my_edge_in, const patch_edge& symmetry_edge_in,
-	     int my_edge_sample_ipar,      int symmetry_edge_sample_ipar,
+	    (const patch_edge& my_edge_in, const patch_edge& other_edge_in,
+	     int my_edge_sample_ipar,      int other_edge_sample_ipar,
 	     bool ipar_map_is_plus);
 
 	~symmetry_ghost_zone();
@@ -577,10 +581,6 @@ private:
 	symmetry_ghost_zone& operator=(const symmetry_ghost_zone& rhs);
 
 private:
-	// symmetry mapping --> interior of which patch?  which edge?
-	patch& symmetry_patch_;
-	const patch_edge& symmetry_edge_;
-
 	// symmetry mappings for (iperp,ipar)
 	// ... we own these objects
 	const jtutil::cpm_map<fp>* iperp_map_;
@@ -626,7 +626,7 @@ public:
 	// (flags specify which part(s) of the ghost zone we want)
 	//
 	// ... the present implementation only supports the case where
-	//     all of the flags are set
+	//     both flags are set
 	//
 	void synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 			 bool want_corners = true,
@@ -644,17 +644,13 @@ public:
 	// allocate internal buffers, compute Jacobian
 	//
 	// ... the present implementation only supports the case where
-	//     all of the flags are set
+	//     both flags are set
 	//
 	void compute_Jacobian(int ghosted_min_gfn, int ghosted_max_gfn,
 			      bool want_corners = true,
 			      bool want_noncorner = true)
 		const;
 
-	// to which patch/edge do the y points in this Jacobian row belong?
-	patch&            Jacobian_y_patch() const { return other_patch_; }
-	const patch_edge& Jacobian_y_edge()  const { return other_edge_; }
-
 	// what is the [min,max] range of m for this ghost zone?
 	int Jacobian_min_y_ipar_m() const { return Jacobian_min_y_ipar_m_; }
 	int Jacobian_max_y_ipar_m() const { return Jacobian_max_y_ipar_m_; }
@@ -694,10 +690,6 @@ public:
 	//
 
 public:
-	// basic connectivity info
-	patch& other_patch() const { return other_patch_; }
-	const patch_edge& other_edge() const { return other_edge_; }
-
 	// check consistency of this and the other patch's ghost zones
 	// and patch_interp objects
 	void assert_fully_setup() const;
@@ -741,9 +733,6 @@ private:
 	interpatch_ghost_zone& operator=(const interpatch_ghost_zone& rhs);
 
 private:
-	patch& other_patch_;
-	const patch_edge& other_edge_;
-
 	//
 	// all the remaining pointers are initialized to NULL pointers
 	// in our constructor, then finally allocated and set up by
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index fe520b0..92262ef 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -10,8 +10,7 @@
 //
 // patch::ghost_zone_on_edge
 // patch::corner_ghost_zone_containing_point
-// patch::symmetry_ghost_zone_on_edge
-// patch::interpatch_ghost_zone_on_edge
+// patch::ghost_zone_containing_point
 // patch::create_mirror_symmetry_ghost_zone
 // patch::create_periodic_symmetry_ghost_zone
 // patch::create_interpatch_ghost_zone
@@ -214,33 +213,43 @@ return is_in_rho_ghost_zone ? rho_gz : sigma_gz;
 //******************************************************************************
 
 //
-// These functions verify that the ghost zone on a specified edge is
-// indeed of the specified type, and return a reference to it as the
-// appropriate derived class.
+// This function determines which ghost zone contains a specified point.
+// For a corner point between two symmetry ghost zones, it's unspecified
+// which ghost zone will be chosen.
 //
-symmetry_ghost_zone& patch::symmetry_ghost_zone_on_edge
-	(const patch_edge &e)
-	const
-{
-ghost_zone &g = ghost_zone_on_edge(e);
-assert(g.is_symmetry());
-return static_cast<symmetry_ghost_zone &>(g);
-}
-
-//**************************************
-
+// If the point isn't in any ghost zone of this patch, an error_exit()
+// is done.
 //
-// This function verifies that the ghost zone on a specified edge
-// is indeed interpatch, and returns a reference to it as an
-//  interpatch_ghost_zone .
+// Arguments:
+// irho,isigma = Specify the point.
 //
-interpatch_ghost_zone& patch::interpatch_ghost_zone_on_edge
-	(const patch_edge &e)
+// Results:
+// This function returns (a reference to) the desired ghost zone.
+ghost_zone& patch::ghost_zone_containing_point(int irho, int isigma)
 	const
 {
-ghost_zone &g = ghost_zone_on_edge(e);
-assert(g.is_interpatch());
-return static_cast<interpatch_ghost_zone &>(g);
+	// n.b. these loops must use _int_ variables for the loop
+	//      to terminate!
+	for (int want_min = false ; want_min <= true ; ++want_min)
+	{
+	for (int want_rho = false ; want_rho <= true ; ++want_rho)
+	{
+	const patch_edge& e = minmax_ang_patch_edge(want_min, want_rho);
+	const int iperp = e.iperp_of_irho_isigma(irho, isigma);
+	const int ipar  = e.ipar_of_irho_isigma (irho, isigma);
+
+	ghost_zone& gz = minmax_ang_ghost_zone(want_min, want_rho);
+	if (gz.is_in_ghost_zone(iperp, ipar))
+	   then return gz;
+	}
+	}
+
+error_exit(ERROR_EXIT,
+"***** patch::ghost_zone_containing_point():\n"
+"        no ghost zone contains point (this should never happen)!\n"
+"        patch=%s irho=%d isigma=%d\n"
+,
+	   name(), irho, isigma);				/*NOTREACHED*/
 }
 
 //******************************************************************************
@@ -267,20 +276,20 @@ set_ghost_zone(my_edge, temp);
 // ghost zone and properly links this to/from the patch.
 //
 void patch::create_periodic_symmetry_ghost_zone
-	(const patch_edge& my_edge, const patch_edge& symmetry_edge,
+	(const patch_edge& my_edge, const patch_edge& other_edge,
 	 bool is_ipar_map_plus)
 {
 // make sure we belong to the right patch
 assert(& my_edge.my_patch() == this);
 
 int my_sample_ipar = my_edge.min_ipar_without_corners();
-int symmetry_sample_ipar
-	= is_ipar_map_plus ? symmetry_edge.min_ipar_without_corners()
-			   : symmetry_edge.max_ipar_without_corners();
+int other_sample_ipar = is_ipar_map_plus
+			? other_edge.min_ipar_without_corners()
+			: other_edge.max_ipar_without_corners();
 
 symmetry_ghost_zone *temp
-	= new symmetry_ghost_zone(my_edge,        symmetry_edge,
-				  my_sample_ipar, symmetry_sample_ipar,
+	= new symmetry_ghost_zone(my_edge,        other_edge,
+				  my_sample_ipar, other_sample_ipar,
 				  is_ipar_map_plus);
 set_ghost_zone(my_edge, temp);
 }
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 3724188..8a90d56 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -163,6 +163,8 @@ public:
 	// n.b. this does *not* compare any of the gridfn data!
 	bool operator==(const patch& other_patch) const
 		{ return patch_number() == other_patch.patch_number(); }
+	bool operator!=(const patch& other_patch) const
+		{ return ! operator==(other_patch); }
 
 	// (rho,sigma,tau) coordinates as singleton coordinate sets
 	local_coords::coords_set coords_set_rho() const
@@ -312,9 +314,6 @@ public:
 		assert(max_sigma_ghost_zone_ != NULL);
 		return *max_sigma_ghost_zone_;
 		}
-
-	// ... these fns are defined explicitly (unlike for grid:: stuff)
-	//     because they're actually used by Jacobian code in ../gr/
 	ghost_zone& minmax_rho_ghost_zone(bool want_min)
 		const
 		{
@@ -337,21 +336,15 @@ public:
 
 	ghost_zone& ghost_zone_on_edge(const patch_edge &e) const;
 
-	// which of the two ghost zones at a specifieid corner,
+	// which of the two ghost zones at a specified corner,
 	// contains a specified point?
 	ghost_zone& corner_ghost_zone_containing_point
 		(bool rho_is_min, bool sigma_is_min,	// specifies corner
 		 int irho, int isigma)			// specifies point
 		const;
 
-	// get ghost zone on specified edge,
-	// assert() that it is indeed of the specified type,
-	// static_cast<> to {symmetry,interpatch}_ghost_zone&
-	symmetry_ghost_zone& symmetry_ghost_zone_on_edge
-		(const patch_edge &e)
-		const;
-	interpatch_ghost_zone& interpatch_ghost_zone_on_edge
-		(const patch_edge &e)
+	// which ghost zone contains a specified point?
+	ghost_zone& ghost_zone_containing_point(int irho, int isigma)
 		const;
 
 
@@ -367,7 +360,7 @@ public:
 	// assert() that this ghost zone hasn't been set up yet,
 	// then set it up as periodic-symmetry
 	void create_periodic_symmetry_ghost_zone
-		(const patch_edge& my_edge, const patch_edge& symmetry_edge,
+		(const patch_edge& my_edge, const patch_edge& other_edge,
 		 bool ipar_map_is_plus);
 
 	// assert() that this ghost zone hasn't been set up yet,
diff --git a/src/patch/patch_edge.hh b/src/patch/patch_edge.hh
index 4538fda..b1586b5 100644
--- a/src/patch/patch_edge.hh
+++ b/src/patch/patch_edge.hh
@@ -68,6 +68,8 @@ public:
 		       && (   is_rho() == other_edge.  is_rho() )
 		       && (   is_min() == other_edge.  is_min() );
 		}
+	bool operator!=(const patch_edge& other_edge) const
+		{ return ! operator==(other_edge); }
 
 
 	//
@@ -85,6 +87,11 @@ public:
 		return my_patch()
 		       .minmax_ang_patch_edge(grid::side_is_max, par_is_rho());
 		}
+	const patch_edge& minmax_par_adjacent_edge(bool want_min) const
+		{
+		return want_min ? min_par_adjacent_edge()
+				: max_par_adjacent_edge();
+		}
 
 
 	//
@@ -222,7 +229,7 @@ public:
 		return    ipar_is_in_min_ipar_corner(ipar)
 		       || ipar_is_in_max_ipar_corner(ipar);
 		}
-	bool ipar_is_in_non_corner(int ipar) const
+	bool ipar_is_in_noncorner(int ipar) const
 		{
 		return    (ipar >= min_ipar_without_corners())
 		       && (ipar <= max_ipar_without_corners());
@@ -231,12 +238,12 @@ public:
 	// convenience function selecting amongst the above
 	// membership predicates
 	bool ipar_is_in_selected_part(bool want_corners,
-				      bool want_non_corner,
+				      bool want_noncorner,
 				      int ipar)
 		const
 		{
-		return    (want_corners    && ipar_is_in_corner    (ipar))
-		       || (want_non_corner && ipar_is_in_non_corner(ipar));
+		return    (want_corners   && ipar_is_in_corner   (ipar))
+		       || (want_noncorner && ipar_is_in_noncorner(ipar));
 		}
 
 	// outer (farthest from patch center) iperp of nominal grid
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index c19dd42..2005d89 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -55,18 +55,20 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 			   const jtutil::array1d<int>& min_parindex_array_in,
 			   const jtutil::array1d<int>& max_parindex_array_in,
 			   const jtutil::array2d<fp>& interp_par_in,
-			   bool ok_to_use_min_par_corner,
-			   bool ok_to_use_max_par_corner,
+			   bool ok_to_use_min_par_ghost_zone,
+			   bool ok_to_use_max_par_ghost_zone,
 			   int interp_handle_in, int interp_par_table_handle_in)
 	: my_patch_(my_edge_in.my_patch()),
 	  my_edge_(my_edge_in),
 	  min_gfn_(my_patch().ghosted_min_gfn()),
 	  max_gfn_(my_patch().ghosted_max_gfn()),
+	  ok_to_use_min_par_ghost_zone_(ok_to_use_min_par_ghost_zone),
+	  ok_to_use_max_par_ghost_zone_(ok_to_use_max_par_ghost_zone),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
-	  min_ipar_(ok_to_use_min_par_corner
+	  min_ipar_(ok_to_use_min_par_ghost_zone
 		    ? my_edge_in.min_ipar_with_corners()
 		    : my_edge_in.min_ipar_without_corners()),
-	  max_ipar_(ok_to_use_max_par_corner
+	  max_ipar_(ok_to_use_max_par_ghost_zone
 		    ? my_edge_in.max_ipar_with_corners()
 		    : my_edge_in.max_ipar_without_corners()),
 	  min_parindex_array_(min_parindex_array_in),
@@ -220,13 +222,75 @@ const CCTK_INT *gridfn_type_codes_ptr
 					 N_gridfns,
 					    gridfn_type_codes_ptr,
 					    interp_buffer);
+
+	if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
+	   then {
+		// look in interpolator output table entries
+		// to see *which* point is out-of-range
+		CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end;
+		if (    (Util_TableGetInt(interp_par_table_handle_,
+					  &out_of_range_pt,
+					  "out_of_range_pt") < 0)
+		     || (Util_TableGetInt(interp_par_table_handle_,
+					  &out_of_range_end,
+					  "out_of_range_end") < 0)    )
+		   then error_exit(ERROR_EXIT,
+"***** patch_interp::interpolate():\n"
+"        point out of range in interpatch interpolation!\n"
+"        ==> something's wrong with the patch system!\n"
+"        (unable to get info about which point is out of range:\n"
+"         maybe an interpolator problem?)\n"
+"        iperp=%d of [%d,%d]\n"
+"        my_patch()=\"%s\" my_edge()=\"%s\"\n"
+,
+				   iperp, min_iperp(), max_iperp(),
+				   my_patch().name(), my_edge().name());
+								/*NOTREACHED*/
+
+		assert(out_of_range_pt >= 0);
+		assert(out_of_range_pt < N_interp_points);
+		const fp out_of_range_interp_par
+		   = interp_coords_ptr[out_of_range_pt];
+
+		const fp gridfn_min_par = gridfn_coord_origin_;
+		const fp gridfn_max_par
+		   = gridfn_coord_origin_
+		     + (N_gridfn_data_points-1)*gridfn_coord_delta_;
+
+		error_exit(ERROR_EXIT,
+"***** patch_interp::interpolate():\n"
+"        point out of range in interpatch interpolation!\n"
+"        ==> something's wrong with the patch system!\n"
+"        iperp=%d of [%d,%d]\n"
+"        my_patch()=\"%s\" my_edge()=\"%s\"\n"
+"        bad point is pt=%d of N_interp_points=%d, end=%d\n"
+"        ==> interp_par=%g\n"
+"        gridfn_coord_origin_=%g gridfn_coord_delta_=%g\n"
+"        ok_to_use_[min,max]_par_ghost_zone_=(int)[%d,%d]\n"
+"        [min,max]_iperp_=[%d,%d] [min,max]_ipar_=[%d,%d]\n"
+"        ==> gridfn_[min,max]_par=[%g,%g]\n"
+,
+			   iperp, min_iperp(), max_iperp(),
+			   my_patch().name(), my_edge().name(),
+			   out_of_range_pt, N_interp_points, out_of_range_end,
+			   double(out_of_range_interp_par),
+			   double(gridfn_coord_origin_),
+			   double(gridfn_coord_delta_),
+			   int(ok_to_use_min_par_ghost_zone_),
+			   int(ok_to_use_max_par_ghost_zone_),
+			   min_iperp_, max_iperp_, min_ipar_, max_ipar_,
+			   double(gridfn_min_par), double(gridfn_max_par));
+								/*NOTREACHED*/
+		}
+
 	if (status < 0)
 	   then error_exit(ERROR_EXIT,
 "***** patch_interp::interpolate():\n"
 "        error return %d from interpolator at iperp=%d of [%d,%d]!\n"
+"        my_patch()=\"%s\" my_edge()=\"%s\"\n"
 ,
-			   status, iperp, min_iperp(), max_iperp());
-								/*NOTREACHED*/
+			   status, iperp, min_iperp(), max_iperp(),
+			   my_patch().name(), my_edge().name()); /*NOTREACHED*/
 	}
 }
 
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 5693a4a..8ad3357 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -31,8 +31,8 @@
 // from its owning patch for use by another patch's ghost_zone object
 // (in setting up the gridfn in the other ghost zone).  A patch_interp
 // object deals only in its own patch's coordinates; other code elsewhere
-// (in practice in ghost_zone::) is responsible for translating other
-// patch's coordinates into our coordinates.
+// (in practice in interpatch_ghost_zone::) is responsible for translating
+// other patch's coordinates into our coordinates.
 //
 
 //
@@ -60,10 +60,10 @@
 //              1    2    3    4    5    6    7    8    9
 // iperp=10           (2a)   (3b)   (4c)
 // iperp=11          (2d)   (3e)  (4f)   (5g)
-// where the (2a) etc are the interpolation points, with 2 etc being
-// the parindex values and a etc being unique identifiers used in our
-// description below.  In terms of our member data, this interpolation
-// region would be described by
+// where the (2a)-(5g) are the interpolation points, with 2-5 being the
+// parindex values and a-g being unique identifiers used in our description
+// below.  In terms of our member data, this interpolation region would
+// be described by
 //	[min,max]_iperp_=[10,11]
 //	[min,max]_ipar_=[1,9]
 //	[min,max]_parindex_array_(10)=[2,5]
@@ -141,9 +141,11 @@ public:
 	void verify_Jacobian_sparsity_pattern_ok() const;
 
 	//
-	// the API for the remaining Jacobian functions implicitly
+	// The API for the remaining Jacobian functions implicitly
 	// assumes that the Jacobian sparsity pattern is "ok" as
-	// verified by  verify_Jacobian_sparsity_pattern_ok()
+	// verified by  verify_Jacobian_sparsity_pattern_ok() ,
+	// and in particular that  [min,max]_ipar_m  are independent
+	// of iperp and parindex.
 	//
 
 	// get [min,max] ipar m coordinates of interpolation molecules
@@ -212,9 +214,9 @@ public:
 	//	  are unused.  We keep a reference to this array, so it
 	//	  should have a lifetime at last as long as that of this
 	//	  object.
-	// ok_to_use_[min,max]_par_corner
+	// ok_to_use_[min,max]_par_ghost_zone
 	//	= Boolean flags saying whether or not we should use gridfn
-	//	  data from the [min,max]_par corners in the interpolation.
+	//	  data from the [min,max]_par ghost zones in the interpolation.
 	// interp_handle_in = Cactus handle to the interpolation operator.
 	// interp_par_table_handle_in
 	//	= Cactus handle to a Cactus key/value table giving
@@ -223,17 +225,7 @@ public:
 	//	  the clone, so the original table is not modified by
 	//	  any actions of this class.
 	//
-	// Note that this constructor looks at what type of ghost zones
-	// are on the adjacent edges to decide how to handle the corners:
-	// * if an adjacent ghost zone is a symmetry ghost zone
-	//   we assume it's already been mirrored by the time
-	//   we get set up ==> we can (and do) use the data from it
-	// * if an adjacent ghost zone is an interpatch ghost zone,
-	//   we can't assume it's already been interpatch-interpolated
-	//   by the time we're called ==> we don't use data from it
-	// Thus the adjacent-side  ghost_zone  objects must already exist!
-	//
-	// This constructor also requires that this patch's gridfns already
+	// This constructor requires that this patch's gridfns already
 	// exist, since we size various arrays based on the patch's min/max
 	// ghosted gfn.
 	//
@@ -242,8 +234,8 @@ public:
 		     const jtutil::array1d<int>& min_parindex_array_in,
 		     const jtutil::array1d<int>& max_parindex_array_in,
 		     const jtutil::array2d<fp>& interp_par_in,
-		     bool ok_to_use_min_par_corner,
-		     bool ok_to_use_max_par_corner,
+		     bool ok_to_use_min_par_ghost_zone,
+		     bool ok_to_use_max_par_ghost_zone,
 		     int interp_handle_in, int interp_par_table_handle_in);
 	~patch_interp();
 
@@ -266,6 +258,10 @@ private:
 	// (any given interpolate() call may specify a subrange)
 	const int min_gfn_, max_gfn_;
 
+	// these are strictly speaking redundant
+	// but we keep them for use in debugging
+	bool ok_to_use_min_par_ghost_zone_, ok_to_use_max_par_ghost_zone_;
+
 	// patch interpolation region,
 	// i.e. range of (iperp,ipar) in this patch from which
 	// we will use gridfn data in interpolation
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index f4881e9..ac42d6e 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -5,11 +5,12 @@
 // patch_system::~patch_system
 // patch_system::create_patches
 // patch_system::setup_gridfn_storage
-// patch_system::setup_full_sphere_patch_system
-// patch_system::setup_plus_z_hemisphere_patch_system
-// patch_system::setup_plus_xy_quadrant_patch_system
-// patch_system::setup_plus_xz_quadrant_patch_system
-// patch_system::setup_octant_patch_system
+// patch_system::interlink_full_sphere_patch_system
+// patch_system::interlink_plus_z_hemisphere_patch_system
+// patch_system::interlink_plus_xy_quadrant_patch_system
+// patch_system::interlink_plus_xz_quadrant_patch_system
+// patch_system::interlink_plus_xyz_octant_patch_system
+// patch_system::create_periodic_symmetry_ghost_zones
 // patch_system::create_interpatch_ghost_zones
 // patch_system::finish_interpatch_setup
 // patch_system::assert_all_ghost_zones_fully_setup
@@ -32,6 +33,8 @@
 // patch_system::compute_synchronize_Jacobian
 // patch_system::synchronize_Jacobian_global_minmax_ym
 // patch_system::synchronize_Jacobian
+/// patch_system::fold_Jacobian
+/// patch_system::ghost_zone_Jacobian
 //
 
 #include <cstdio>
@@ -41,7 +44,7 @@
 #include <assert.h>
 #include <limits.h>
 using std::printf;
-using std::strcmp;
+using std::fprintf;
 
 #include "cctk.h"
 
@@ -846,9 +849,11 @@ py.create_interpatch_ghost_zone(ey, ex, N_overlap_points);
 {
 const patch_edge& ex = px.edge_adjacent_to_patch(py, N_overlap_points);
 const patch_edge& ey = py.edge_adjacent_to_patch(px, N_overlap_points);
-px.interpatch_ghost_zone_on_edge(ex)
+px.ghost_zone_on_edge(ex)
+  .cast_to_interpatch_ghost_zone()
   .finish_setup(interp_handle, interp_par_table_handle);
-py.interpatch_ghost_zone_on_edge(ey)
+py.ghost_zone_on_edge(ey)
+  .cast_to_interpatch_ghost_zone()
   .finish_setup(interp_handle, interp_par_table_handle);
 }
 
@@ -1455,7 +1460,8 @@ max_ym = global_max_ym_;
 
 //
 // Given that  compute_synchronize_Jacobian()  has been called, this
-// function computes a single row of the Jacobian:
+// function computes a single row of the Jacobian, taking into account
+// synchronize()'s 3-phase algorithm:
 // - It returns the edge to which the y point belongs (the caller can get
 //   the patch from this edge).
 // - It stores y_iperp and y_posn and min/max ym in the named arguments.
@@ -1463,10 +1469,32 @@ max_ym = global_max_ym_;
 //	partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
 //	-------------------------------------------------------------
 //	     partial gridfn(ghosted_gfn, py, y_iperp, y_posn+ym)
-//   (taking into account synchronize()'s full 3-phase algorithm)
 //   in the caller-supplied buffer
 //	Jacobian_buffer(ym)
-//   for each  ym  in the min/max ym range
+//   for each  ym  in the min/max ym range.
+//
+// In practice, the main task of this function is taking into account
+// synchronize()'s 3-phase algorithm.  There are several cases:
+// - ghost zone is symmetry && x point is in non-corner
+//   ==> x value was computed by a phase 1 symmetry operation,
+//       using (only) nominal-grid data
+//   ==> overall Jacobian(x,y) = ghost zone Jacobian(x,y)
+// - ghost zone is symmetry && x point is in corner
+//   --> x value was computed by a phase 3 symmetry operation,
+//       from some point (call it z),
+//       ==> overall Jacobian(x,y) = overall Jacobian(z,y)
+//       ==> call this function recursively to get z's Jacobian
+//           (z must be in the noncorner part of some ghost zone,
+//	      so this won't lead to infinite recursion)
+// - ghost zone is interpatch
+//   ==> x value was computed by a phase 2 interpatch interpolation
+//       - using (only) nominal-grid data
+//         ==> overall Jacobian(x,y) = ghost zone Jacobian(x,y)
+//       - using a mixture of nominal-grid data
+//         and data computed by a phase 1 symmetry operation
+//         ==> overall Jacobian(x,y) = "fold" ghost zone Jacobian(x,y)
+//                                     to take the phase 1 symmetry
+//                                     operation into account
 //
 const patch_edge&
   patch_system::synchronize_Jacobian(const ghost_zone& xgz,
@@ -1476,31 +1504,265 @@ const patch_edge&
 				     jtutil::array1d<fp>& Jacobian_buffer)
 	const
 {
+const patch_edge& xe = xgz.my_edge();
+
+if (xgz.is_symmetry() && xe.ipar_is_in_noncorner(x_ipar))
+   then {
+	// ghost zone is symmetry && x point is in non-corner
+	// ==> x value was computed by a phase 1 symmetry operation,
+	//     using (only) nominal-grid data
+	// ==> overall Jacobian(x,y) = ghost zone Jacobian(x,y)
+	return ghost_zone_Jacobian(xgz,
+				   x_iperp, x_ipar,
+				   y_iperp,
+				   y_posn, min_ym, max_ym,
+				   Jacobian_buffer);
+	}
+
+else if (xgz.is_symmetry() && xe.ipar_is_in_corner(x_ipar))
+   then {
+	// ghost zone is symmetry && x point is in corner
+	// --> x value was computed by a phase 3 symmetry operation,
+	//     from some point (call it z),
+	//     ==> overall Jacobian(x,y) = overall Jacobian(z,y)
+	//     ==> call this function recursively to get z's Jacobian
+	//         (z must be in the noncorner part of some ghost zone,
+	//          so this won't lead to infinite recursion)
+
+	const patch&      zp = xgz.other_patch();
+	const patch_edge& ze = xgz.other_edge();
+	const symmetry_ghost_zone& xsgz = xgz.cast_to_symmetry_ghost_zone();
+	const int z_iperp = xsgz.iperp_map_of_iperp(x_iperp);
+	const int z_ipar  = xsgz.ipar_map_of_ipar  (x_ipar);
+
+	//
+	// Computing z's edge/ghost zone is tricky.  For example:
+	//	                        |
+	//	           p1         e3|e4          p2
+	//	                        |
+	//	                        |   z
+	//	-----------e1-----------+------------e2---------
+	//	                        |   x
+	//	                        |
+	// Here the point x in the corner of p1's e1 ghost zone,
+	// is computed by the phase 3 symmetry operation (a reflection
+	// about e1) from z.  Thus zp == p1 and ze == e1.
+	//
+	// But we need to "turn the corner" to compute z's "true" edge
+	// e3 (so we can recursively call this function to compute z's
+	// Jacobian).  Thus we explicitly check which ghost zone of p1
+	// (here the e3 ghost zone) contains the point z.
+	//
+	const int z_irho   = ze.  irho_of_iperp_ipar(z_iperp, z_ipar);
+	const int z_isigma = ze.isigma_of_iperp_ipar(z_iperp, z_ipar);
+	const ghost_zone& true_zgz
+		= zp.ghost_zone_containing_point(z_irho, z_isigma);
+	const patch_edge& true_ze = true_zgz.my_edge();
+	const int true_z_iperp = true_ze.iperp_of_irho_isigma(z_irho, z_isigma);
+	const int true_z_ipar  = true_ze. ipar_of_irho_isigma(z_irho, z_isigma);
+
+	// make sure we have the right ghost zone!
+	assert( true_zgz.is_in_ghost_zone(true_z_iperp, true_z_ipar) );
+
+	return synchronize_Jacobian(true_zgz,
+				    true_z_iperp, true_z_ipar,
+				    y_iperp,
+				    y_posn, min_ym, max_ym,
+				    Jacobian_buffer);
+	}
+
+else if (xgz.is_interpatch())
+   then {
+	// ghost zone is interpatch
+	// ==> x value was computed by a phase 2 interpatch interpolation
+	//     - using (only) nominal-grid data
+	//       ==> overall Jacobian(x,y) = ghost zone Jacobian(x,y)
+	//     - using a mixture of nominal-grid data
+	//       and data computed by a phase 1 symmetry operation
+	//       ==> overall Jacobian(x,y) = "fold" ghost zone Jacobian(x,y)
+	//                                   to take the phase 1 symmetry
+	//                                   operation into account
+	//
+	// For example,
+	//	                        |
+	//	           xp         xe|ye a        yp
+	//	                        |   b
+	//	                        |  xc
+	//	----------xse-----------+---d-------yse----------
+	//	                        |   e
+	//	                        |
+	// here point x is computed by interpatch-interpolating in the
+	// par direction from the 5 y points abcde.  e is outside the
+	// nominal grid, so its Jacobian must be "folded" over to c.
+	// Notice that this "folding" must be done about the edge yse,
+	// *not* about ye itself.
+
+	// Jacobian of the phase 2 interpatch interpolation
+	const patch_edge& ye = ghost_zone_Jacobian(xgz,
+						   x_iperp, x_ipar,
+						   y_iperp,
+						   y_posn, min_ym, max_ym,
+						   Jacobian_buffer);
+	const patch& yp = ye.my_patch();
+	const int min_y_ipar = y_posn + min_ym;
+	const int max_y_ipar = y_posn + max_ym;
+
+	// fold any points in the Jacobian outside the nominal grid
+	if (ye.ipar_is_in_min_ipar_corner(min_y_ipar))
+	   then {
+		fold_Jacobian(ye, ye.min_par_adjacent_edge(),
+			      y_iperp,
+			      y_posn, min_ym, max_ym,
+			      min_ym, ye.min_ipar_corner__max_ipar() - y_posn,
+			      Jacobian_buffer);
+		min_ym = ye.min_ipar_without_corners() - y_posn;
+		}
+	if (ye.ipar_is_in_max_ipar_corner(max_y_ipar))
+	   then {
+		fold_Jacobian(ye, ye.max_par_adjacent_edge(),
+			      y_iperp,
+			      y_posn, min_ym, max_ym,
+			      ye.max_ipar_corner__min_ipar() - y_posn, max_ym,
+			      Jacobian_buffer);
+		max_ym = ye.max_ipar_without_corners() - y_posn;
+		}
+
+	return ye;
+	}
+
+else	error_exit(PANIC_EXIT,
+"***** patch_system::synchronize_Jacobian():\n"
+"        don't know what to do with ghost zone (this should never happen)!\n"
+"        xgz.my_patch()=\"%s\"    xe=xgz.my_edge()=\"%s\"\n"
+"        xgz.other_patch()=\"%s\"    xgz.other_edge()=\"%s\"\n"
+"        xgz.is_symmetry()=(int)%d xgz.is_interpatch()=(int)%d\n"
+"        x_iperp=%d x_ipar=%d\n"
+"        xe.ipar_is_in_{min,max}_ipar_corner(x_ipar)=(int){%d,%d}\n"
+"        xe.ipar_is_in_{corner,noncorner}(x_ipar)=(int){%d,%d}\n"
+,
+		   xgz.my_patch().name(),    xe.name(),
+		   xgz.other_patch().name(), xgz.other_edge().name(),
+		   int(xgz.is_symmetry()), int(xgz.is_interpatch()),
+		   x_iperp, x_ipar,
+		   xe.ipar_is_in_min_ipar_corner(x_ipar),
+		   xe.ipar_is_in_max_ipar_corner(x_ipar),
+		   xe.ipar_is_in_corner(x_ipar),
+		   xe.ipar_is_in_noncorner(x_ipar));		/*NOTREACHED*/
+}
+
+//******************************************************************************
+
 //
-// Implementation notes:
+// This function "folds" part of a(n interpatch) Jacobian row to take
+// a symmetry operation into account.  For example:
+//	         |
+//	         |e_Jac
+//	         |             p
+//	         |   a
+//	         |   b
+//	         |   c=y
+//	---------+---d-------e_fold-------
+//	         |   e=x    sgz_fold
+//	         |
+// Here the Jacobian abcde is to be "folded", because e is outside the
+// nominal grid (its Jacobian must be "folded" over to c).
+//
+// Notice that the folding (about the edge e_fold) is in the par direction
+// with respect to e_Jac, but the perp direction with respect to e_fold.
+// Since e_fold and e_Jac are adjacent edges, 
+//	e_Jac (iperp,ipar) == e_fold (ipar,iperp)
 //
-// There are several cases for the 3-phase algorithm.
-// - ghost zone is symmetry && x point is in non-corner
-//   --> x value was computed by a phase 1 symmetry operation,
-//       from nominal-grid data
-//   ==> use the ghost zone's Jacobian unchanged
-// - ghost zone is symmetry && x point is in corner
-//   --> x value was computed by a phase 3 symmetry operation,
-//       - from a symmetry ghost zone,
-//         which must have been computed in phase 1 from nominal-grid data
-//         ==> use the two symmetry ghost zones' combined Jacobian
-//       - from an interpatch ghost zone
-//         ==> use that interpatch ghost zone's Jacobian
-// - ghost zone is interpatch
-//   --> x value was computed by a phase 2 interpatch interpolation
-//       - from (only) nominal-grid data
-//         ==> use that interpatch ghost zone's Jacobian unchanged
-//       - from a mixture of nominal-grid data
-//         and data computed by a phase 1 symmetry operation
-//         ==> "fold" the interpatch ghost zone's Jacobian
-//             to reflect the phase 1 symmetry operation
+// Arguments:
+// e_Jac = edge which the Jacobian lies along
+// e_fold = edge about which to fold; the corresponding ghost zone must be
+//	    symmetry ghost zone, and at present we only support the case
+//	    where this is a "local" (mirror-image) symmetry ghost zone
+// iperp = iperp-wrt-e_Jac coordinate of Jacobian
+// posn = ipar-wrt-e_Jac coordinate of Jacobian molecule reference point
+// [min,max]_m = range of ipar-wrt-e_Jac molecule m in Jacobian
+// [min,max]_fold_m = range of ipar-wrt-e_Jac molecule m which is to folded;
+//		      this must be a subrange of [min,max]_m
+//
+void patch_system::fold_Jacobian(const patch_edge& e_Jac,
+				 const patch_edge& e_fold,
+				 int iperp,
+				 int posn, int min_m, int max_m,
+				 int min_fold_m, int max_fold_m,
+				 jtutil::array1d<fp>& Jacobian_buffer)
+	const
+{
+// check that [min,max]_fold_m is a subrange of [min,max]_m
+assert( min_fold_m >= min_m );
+assert( min_fold_m <= max_m );
+assert( max_fold_m >= min_m );
+assert( max_fold_m <= max_m );
+
+const patch& p = e_fold.my_patch();
+assert( e_Jac.my_patch() == p );
+
+const symmetry_ghost_zone& sgz_fold = p.ghost_zone_on_edge(e_fold)
+				       .cast_to_symmetry_ghost_zone();
+
 //
+// At present we only handle the case show in the comments above,
+// where sgz_fold is a local (mirror-image) symmetry, i.e. where
+// y is guaranteed to be within the molecule abcde.
+//
+if (sgz_fold.other_edge() != e_fold)
+   then error_exit(ERROR_EXIT,
+"***** patch_system::fold_Jacobian()\n"
+"        implementation restriction: at present we only handle folding\n"
+"        via \"local\" (mirror-image) symmetries!\n"
+"        p=\"%s\" e_Jac=\"%s\" e_fold=\"%s\"\n"
+"        but sgz_fold.other_edge()=\"%s\" != e_fold\n"
+,
+		   p.name(), e_Jac.name(), e_fold.name(),
+		   sgz_fold.other_edge().name());		/*NOTREACHED*/
+
+	for (int xm = min_fold_m ; xm <= max_fold_m ; ++xm)
+	{
+	const int x_Jac_ipar = posn + xm;	// x ipar wrt e_Jac
+	const int x_fold_iperp = x_Jac_ipar;	// ... == iperp wrt e_fold
+
+	const int y_fold_iperp = sgz_fold.iperp_map_of_iperp(x_fold_iperp);
+						// y iperp wrt e_fold
+	const int y_Jac_ipar = y_fold_iperp;	// ... == ipar wrt e_Jac
+	const int ym = y_Jac_ipar - posn;
+
+	// check that y is indeed within the molecule
+	assert( ym >= min_m );
+	assert( ym <= max_m );
+
+	// actually "fold" the molecule
+	Jacobian_buffer(ym) += Jacobian_buffer(xm);
+	}
+}
+
+//******************************************************************************
 
+//
+// Given that  compute_synchronize_Jacobian()  has been called, this
+// function computes a single row of the Jacobian of a given ghost zone,
+// *not* taking into account synchronize()'s 3-phase algorithm:
+// - It returns the edge to which the y point belongs (the caller can get
+//   the patch from this edge).
+// - It stores y_iperp and y_posn and min/max ym in the named arguments.
+// - It stores the Jacobian elements
+//	partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
+//	-------------------------------------------------------------
+//	     partial gridfn(ghosted_gfn, py, y_iperp, y_posn+ym)
+//   in the caller-supplied buffer
+//	Jacobian_buffer(ym)
+//   for each  ym  in the min/max ym range
+//
+const patch_edge&
+  patch_system::ghost_zone_Jacobian(const ghost_zone& xgz,
+				    int x_iperp, int x_ipar,
+				    int& y_iperp,
+				    int& y_posn, int& min_ym, int& max_ym,
+				    jtutil::array1d<fp>& Jacobian_buffer)
+	const
+{
 y_iperp = xgz.Jacobian_y_iperp(x_iperp);
 
 y_posn = xgz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
@@ -1512,5 +1774,5 @@ max_ym = xgz.Jacobian_max_y_ipar_m();
 	Jacobian_buffer(ym) = xgz.Jacobian(x_iperp, x_ipar, ym);
 	}
 
-return xgz.Jacobian_y_edge();
+return xgz.other_edge();
 }
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 6c42fa8..8700f14 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -206,6 +206,32 @@ public:
 			       jtutil::array1d<fp>& Jacobian_buffer)
 		const;
 
+	// helper functions for synchronize_Jacobian():
+private:
+	// "fold" (part of) a Jacobian row
+	// to take a symmetry operation into acount
+	// e_Jac = edge which the Jacobian lies along
+	// e_fold = edge about which to fold
+	// [min,max]_m = range of m in the Jacobian
+	// [min,max]_fold_m = range of m to fold
+	//		      (must be a subrange of {min,max}_m)
+	void fold_Jacobian(const patch_edge& e_Jac, const patch_edge& e_fold,
+			   int iperp,
+			   int posn, int min_m, int max_m,
+			   int min_fold_m, int max_fold_m,
+			   jtutil::array1d<fp>& Jacobian_buffer)
+		const;
+
+	// compute the Jacobian of ghost zone's synchronize()
+	// *without* taking into account 3-phase algorithm
+	const patch_edge&
+	  ghost_zone_Jacobian(const ghost_zone& xgz,
+			      int x_iperp, int x_ipar,
+			      int& y_iperp,
+			      int& y_posn, int& min_ym, int& max_ym,
+			      jtutil::array1d<fp>& Jacobian_buffer)
+		const;
+
 
 	//
 	// ***** gridfn operations *****