* merge {xy,xz,yz}_circumference() into a single circumference() function

* add flags to this and to  integrate_gridfn()  so it knows about
  the gridfn symmetry and automagically calculates what the integral
  would have been over the full 2-sphere


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

2 files changed, 114 insertions, 169 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 3adaaed..b72fdaa 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -35,10 +35,7 @@
 // patch_system::gridfn_norms
 // patch_system::ghosted_gridfn_norms
 //
-// patch_system::xy_circumference
-// patch_system::xz_circumference
-// patch_system::yz_circumference
-// patch_system::arc_length_sum
+// patch_system::circumference
 // patch_system::integrate_gridfn
 //
 // patch_system::print_unknown_gridfn
@@ -1546,189 +1543,104 @@ norms.reset();
 
 //
 // This function computes an approximation to the circumference of a
-// surface in the xy plane.  That is, we integrate in the phi direction
-// around the equator (the +/- x and +/- y patches).  Note that we
-// copute the full circumference all around the sphere, even if the
-// patch system only covers a proper subset of this.
+// surface in the xy, xz, or yz plane.  Note that we compute the full
+// circumference all around the sphere, even if the patch system only
+// covers a proper subset of this.
+//
+// We assume that adjacent patches are butt-joined, i.e. that their
+// nominal boundaries just touch.
 //
 // Arguments:
+// plane[] = (in) "xy", "xz", or "yz" to specify the integration plane.
 // ghosted_radius_gfn = (in) The surface radius.
 // g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
 // method = (in) Selects the integration scheme.
 //
-fp patch_system::xy_circumference(int ghosted_radius_gfn,
-				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-						int g_yy_gfn, int g_yz_gfn,
-							      int g_zz_gfn,
-				  enum patch::integration_method method)
+fp patch_system::circumference(const char plane[],
+			       int ghosted_radius_gfn,
+			       int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					     int g_yy_gfn, int g_yz_gfn,
+							   int g_zz_gfn,
+			       enum patch::integration_method method)
 	const
 {
-fp circumference = 0.0;
-
-    for (int pn = 0 ; pn < N_patches() ; ++pn)
-    {
-    const patch& p = ith_patch(pn);
-    if ((p.ctype() == 'x') || (p.ctype() == 'y'))
-       then circumference += p.plane_arc_length('z',
-						ghosted_radius_gfn,
-						g_xx_gfn, g_xy_gfn, g_xz_gfn,
-							  g_yy_gfn, g_yz_gfn,
-								    g_zz_gfn,
-						method);
-    }
-
-switch	(type())
-	{
-case patch_system__full_sphere:						break;
-case patch_system__plus_z_hemisphere:					break;
-case patch_system__plus_xy_quadrant_mirrored:	circumference *= 4.0;	break;
-case patch_system__plus_xy_quadrant_rotating:	circumference *= 4.0;	break;
-case patch_system__plus_xz_quadrant_rotating:	circumference *= 2.0;	break;
-case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
-case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
-default:
-	error_exit(PANIC_EXIT,
-"***** patch_system::xy_circumference(): unknown type=(int)%d!\n"
-"                                        (this should never happen!)\n",
-		   int(type()));				/*NOTREACHED*/
-	}
-
-return circumference;
-}
-
-//******************************************************************************
-
 //
-// This function coputes an approximation to the circumference of a
-// surface in the xz plane.  That is, we integrate in the nu direction
-// around the y=0 plane (the +/- x and +/- z patches).  Note that we
-// copute the full circumference all around the sphere, even if the
-// patch system only covers a proper subset of this.
+// compute arc length around the patch system
 //
-// Arguments:
-// ghosted_radius_gfn = (in) The surface radius.
-// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
-// method = (in) Selects the integration scheme.
-//
-fp patch_system::xz_circumference(int ghosted_radius_gfn,
-				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-						int g_yy_gfn, int g_yz_gfn,
-							      int g_zz_gfn,
-				  enum patch::integration_method method)
-	const
-{
-fp circumference = 0.0;
-
+fp arc_length = 0.0;
     for (int pn = 0 ; pn < N_patches() ; ++pn)
     {
     const patch& p = ith_patch(pn);
-    if ((p.ctype() == 'x') || (p.ctype() == 'z'))
-       then circumference += p.plane_arc_length('y',
-						ghosted_radius_gfn,
-						g_xx_gfn, g_xy_gfn, g_xz_gfn,
-							  g_yy_gfn, g_yz_gfn,
-								    g_zz_gfn,
-						method);
+    if ((p.ctype() == plane[0]) || (p.ctype() == plane[1]))
+       then arc_length += p.plane_arc_length(plane,
+					     ghosted_radius_gfn,
+					     g_xx_gfn, g_xy_gfn, g_xz_gfn,
+						       g_yy_gfn, g_yz_gfn,
+								 g_zz_gfn,
+					     method);
     }
 
-switch	(type())
-	{
-case patch_system__full_sphere:						break;
-case patch_system__plus_z_hemisphere:		circumference *= 2.0;	break;
-case patch_system__plus_xy_quadrant_mirrored:	circumference *= 2.0;	break;
-case patch_system__plus_xy_quadrant_rotating:	circumference *= 2.0;	break;
-case patch_system__plus_xz_quadrant_rotating:	circumference *= 4.0;	break;
-case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
-case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
-default:
-	error_exit(PANIC_EXIT,
-"***** patch_system::xz_circumference(): unknown type=(int)%d!\n"
-"                                        (this should never happen!)\n",
-		   int(type()));				/*NOTREACHED*/
-	}
-
-return circumference;
-}
-
-//******************************************************************************
-
-//
-// This function coputes an approximation to the circumference of a
-// surface in the yz plane.  That is, we integrate in the mu direction
-// around the x=0 plane (the +/- y and +/- z patches).  Note that we
-// copute the full circumference all around the sphere, even if the
-// patch system only covers a proper subset of this.
 //
-// Arguments:
-// ghosted_radius_gfn = (in) The surface radius.
-// g_{xx,xy,yz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
-// method = (in) Selects the integration scheme.
+// correct the arc length
+// for the fact that the patch system may not cover the full 2-sphere
 //
-fp patch_system::yz_circumference(int ghosted_radius_gfn,
-				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-						int g_yy_gfn, int g_yz_gfn,
-							      int g_zz_gfn,
-				  enum patch::integration_method method)
-	const
-{
-fp circumference = 0.0;
-
-    for (int pn = 0 ; pn < N_patches() ; ++pn)
-    {
-    const patch& p = ith_patch(pn);
-    if ((p.ctype() == 'y') || (p.ctype() == 'z'))
-       then circumference += p.plane_arc_length('x',
-						ghosted_radius_gfn,
-						g_xx_gfn, g_xy_gfn, g_xz_gfn,
-							  g_yy_gfn, g_yz_gfn,
-								    g_zz_gfn,
-						method);
-    }
-
 switch	(type())
 	{
-case patch_system__full_sphere:						break;
-case patch_system__plus_z_hemisphere:		circumference *= 2.0;	break;
-case patch_system__plus_xy_quadrant_mirrored:	circumference *= 2.0;	break;
-case patch_system__plus_xy_quadrant_rotating:	circumference *= 2.0;	break;
-case patch_system__plus_xz_quadrant_rotating:	circumference *= 2.0;	break;
-case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
-case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
+case patch_system__full_sphere:
+	break;
+case patch_system__plus_z_hemisphere:
+	arc_length *= ((plane[0] == 'x') && (plane[1] == 'y')) ? 1.0 : 2.0;
+	break;
+case patch_system__plus_xy_quadrant_mirrored:
+case patch_system__plus_xy_quadrant_rotating:
+	arc_length *= ((plane[0] == 'x') && (plane[1] == 'y')) ? 4.0 : 2.0;
+	break;
+case patch_system__plus_xz_quadrant_rotating:
+	arc_length *= ((plane[0] == 'x') && (plane[1] == 'z')) ? 4.0 : 2.0;
+	break;
+case patch_system__plus_xyz_octant_mirrored:
+case patch_system__plus_xyz_octant_rotating:
+	arc_length *= 4.0;
+	break;
 default:
 	error_exit(PANIC_EXIT,
-"***** patch_system::xz_circumference(): unknown type=(int)%d!\n"
-"                                        (this should never happen!)\n",
+"***** patch_system::circumference(): unknown patch system type()=(int)%d!\n"
+"                                     (this should never happen!)\n",
 		   int(type()));				/*NOTREACHED*/
 	}
 
-return circumference;
+return arc_length;
 }
 
 //******************************************************************************
 
 //
 // This function computes an approximation to the (surface) integral of
-// a gridfn over the patch system,
+// a gridfn over the 2-sphere
 //	$\int f(\rho,\sigma) \, dA$
 //		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
 // where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma).
 //
-// FIXME:
-// It would be better to compute the full-sphere integral, but that would
-// require knowing the gridfn's symmetries.
-//
-// FIXME:
-// We silently assume that the patches butt-join, i.e. that they don't
-// overlap, i.e. that patch_overlap_width=1
+// We assume that adjacent patches are butt-joined, i.e. that their
+// nominal boundaries just touch.
 //
 // Arguments:
 // unknown_src_gfn = (in) The gridfn to be integrated.  This may be
 //			  either nominal-grid or ghosted-grid.
+// src_gfn_is_even_across_{xy,xz,yz}_plane
+//	= (in) Boolean flags specifying whether the gridfn to be integrated
+//	       is even (true) or odd (false) across the corresponding planes.
+//	       Only the flags corresponding to boundaries of the patch system
+//	       are used.  For example, for a  plus_z_hemisphere  patch system,
+//	       only the  src_gfn_is_even_across_xy_plane  flag is used.
 // ghosted_radius_gfn = (in) The surface radius.
 // g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
 // method = (in) Selects the integration scheme.
 //
 fp patch_system::integrate_gridfn(int unknown_src_gfn,
+				  bool src_gfn_is_even_across_xy_plane,
+				  bool src_gfn_is_even_across_xz_plane,
+				  bool src_gfn_is_even_across_yz_plane,
 				  int ghosted_radius_gfn,
 				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 						int g_yy_gfn, int g_yz_gfn,
@@ -1736,6 +1648,9 @@ fp patch_system::integrate_gridfn(int unknown_src_gfn,
 				  enum patch::integration_method method)
 	const
 {
+//
+// compute integral over patch system
+//
 fp integral = 0.0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
@@ -1747,6 +1662,40 @@ fp integral = 0.0;
 							   g_zz_gfn,
 				       method);
 	}
+
+//
+// correct the integral
+// for the fact that the patch system may not cover the full 2-sphere
+//
+switch	(type())
+	{
+case patch_system__full_sphere:
+	break;
+case patch_system__plus_z_hemisphere:
+	integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0;
+	break;
+case patch_system__plus_xy_quadrant_mirrored:
+case patch_system__plus_xy_quadrant_rotating:
+	integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0;
+	integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0;
+	break;
+case patch_system__plus_xz_quadrant_rotating:
+	integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0;
+	integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0;
+	break;
+case patch_system__plus_xyz_octant_mirrored:
+case patch_system__plus_xyz_octant_rotating:
+	integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0;
+	integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0;
+	integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0;
+	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** patch_system::integrate_gridfn(): bad patch system type()=(int)%d!\n"
+"                                        (this should never happen!)\n",
+		   int(type()));				/*NOTREACHED*/
+	}
+
 return integral;
 }
 
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index ce7b100..413ecb0 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -286,39 +286,35 @@ public:
 	// compute the circumference of a surface in the {xy, xz, yz} plane
 	// ... note this is the full circumference all around the sphere,
 	//     even if the patch system only covers a proper subset of this
-	fp xy_circumference(int ghosted_radius_gfn,
-			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-					  int g_yy_gfn, int g_yz_gfn,
-							int g_zz_gfn,
-			    enum patch::integration_method method)
-		const;
-	fp xz_circumference(int ghosted_radius_gfn,
-			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-					  int g_yy_gfn, int g_yz_gfn,
-							int g_zz_gfn,
-			    enum patch::integration_method method)
-		const;
-	fp yz_circumference(int ghosted_radius_gfn,
-			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
-					  int g_yy_gfn, int g_yz_gfn,
-							int g_zz_gfn,
-			    enum patch::integration_method method)
+	// ... the implementation assumes adjacent patches are butt-joined
+	// ... plane must be one of "xy", "xz", or "yz"
+	fp circumference(const char plane[],
+			 int ghosted_radius_gfn,
+			 int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+				       int g_yy_gfn, int g_yz_gfn,
+						     int g_zz_gfn,
+			 enum patch::integration_method method)
 		const;
 
-	// compute the surface integral of a gridfn over the patch system,
+	// compute the surface integral of a gridfn over the 2-sphere
 	//	$\int f(\rho,\sigma) \, dA$
 	//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
 	// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma)
 	// ... integration method selected by  method  argument
 	// ... src gridfn may be either nominal-grid or ghosted-grid
-	// ... note integral is only over patch system itself,
-	//     even if this covers a proper subset of the full sphere
-	//     FIXME: it would be better to compute the full-sphere integral
-	//            (but this would require knowing the gridfn's symmetries)
-	//     FIXME: current implementation assumes patches' nominal grids
-	//            form a partition of the region covered, i.e. they
-	//            just touch each other
+	// ... Boolean flags  src_gfn_is_even_across_{xy,xz,yz}_planes
+	//     specify whether the gridfn to be integrated is even (true)
+	//     or odd (false) across the corresponding planes.  Only the
+	//     flags corresponding to boundaries of the patch system are
+	//     used.  For example, for a  plus_z_hemisphere  patch system,
+	//     only the  src_gfn_is_even_across_xy_plane  flag is used.
+	// ... note integral is over the full 2-sphere,
+	//     even if the patch system only covers a proper subset of this
+	// ... the implementation assumes adjacent patches are butt-joined
 	fp integrate_gridfn(int unknown_src_gfn,
+			    bool src_gfn_is_even_across_xy_plane,
+			    bool src_gfn_is_even_across_xz_plane,
+			    bool src_gfn_is_even_across_yz_plane,
 			    int ghosted_radius_gfn,
 			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
 					  int g_yy_gfn, int g_yz_gfn,