add code to compute surface integrals of gridfns over patches

and over the whole patch system
-- note the volume element isn't included yet


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

4 files changed, 458 insertions, 83 deletions

diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index 92262ef..dde6dcd 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -8,6 +8,9 @@
 // x_patch::x_patch
 // y_patch::y_patch
 //
+// patch::decode_integration_method
+// patch::integrate_gridfn
+// patch::integration_coeff
 // patch::ghost_zone_on_edge
 // patch::corner_ghost_zone_containing_point
 // patch::ghost_zone_containing_point
@@ -21,9 +24,11 @@
 
 #include <cstdio>
 #include <cmath>
+#include <cstring>
 #include <assert.h>
 using std::fprintf;
 using std::printf;
+using std::strcmp;
 
 #include "cctk.h"
 
@@ -156,6 +161,137 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
 //******************************************************************************
 
 //
+// This function decodes the character-string name of an integration method
+// into an  enum integration_method .  See the comments in "patch.hh" on the
+// declaration of  enum integration_method  for details on the methods and
+// their character-string names.
+//
+//static
+  enum patch::integration_method
+    patch::decode_integration_method(const char method_string[])
+{
+if	(    STRING_EQUAL(method_string, "trapezoid")
+	  || STRING_EQUAL(method_string, "trapezoid rule")    )
+   then return integration_method__trapezoid;
+else if (    STRING_EQUAL(method_string, "Simpson")
+	  || STRING_EQUAL(method_string, "Simpson's rule")    )
+   then return integration_method__Simpson;
+else if (    STRING_EQUAL(method_string, "Simpson (variant)")
+	  || STRING_EQUAL(method_string, "Simpson's rule (variant)")    )
+   then return integration_method__Simpson_variant;
+else	error_exit(ERROR_EXIT,
+"***** patch::decode_integration_method():\n"
+"        unknown method_string=\"%s\"!\n"
+,
+		   method_string);				/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
+// This function computes an approximation to the (surface) integral of
+// a gridfn over the patch's nominal area,
+//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
+//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
+// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in
+// $(\rho,sigma)$ coordinates.
+//
+// The integration scheme is selected based on the  method  argument.
+//
+// Bugs:
+// The way the integration coefficients are computed is somewhat inefficient.
+//
+fp patch::integrate_gridfn(int src_gfn,
+			   bool area_weighting_flag,
+			   enum integration_method method)
+	const
+{
+fp sum = 0.0;
+	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
+	{
+	for (int isigma = min_isigma() ; isigma <= max_isigma() ; ++isigma)
+	{
+	sum += gridfn(src_gfn, irho,isigma)
+	       * integration_coeff(method,
+				   max_irho()-min_irho(),
+				   irho      -min_irho())
+	       * integration_coeff(method,
+				   max_isigma()-min_isigma(),
+				   isigma      -min_isigma());
+	}
+	}
+return delta_rho() * delta_sigma() * sum;
+}
+
+//******************************************************************************
+
+//
+// This function computes the integration coefficients for
+//  integrate_gridfn() .  That is, if we write
+//	$\int_{x_0}^{x_N} f(x) \, dx
+//		\approx \Delta x \, \sum_{i=0}^N c_i f(x_i)$
+// then this function computes $c_i$.
+//
+// Arguments:
+// method = Specifies the integration method.
+// N = The number of integration intervals.
+// i = Specifies the point at which the coefficient is desired.
+//
+//static
+  fp patch::integration_coeff(enum integration_method method, int N, int i)
+{
+assert(i >= 0);
+assert(i <= N);
+
+switch	(method)
+	{
+case integration_method__trapezoid:
+	if ((i == 0) || (i == N))
+	   then return 0.5;
+	   else return 1.0;
+
+case integration_method__Simpson:
+	if ((N % 2) != 0)
+	   then error_exit(ERROR_EXIT,
+"***** patch::integration_coeff():\n"
+"        Simpson's rule requires N to be even, but N=%d!\n",
+			   N);					/*NOTREACHED*/
+	if	((i == 0) || (i == N))
+	   then return 1.0/3.0;
+	else if ((i % 2) == 0)
+	   then return 2.0/3.0;
+	else	return 4.0/3.0;
+
+case integration_method__Simpson_variant:
+	if (N < 7)
+	   then error_exit(ERROR_EXIT,
+"***** patch::integration_coeff():\n"
+"        Simpson's rule (variant) requires N >= 7, but N=%d!\n",
+			   N);					/*NOTREACHED*/
+	if	((i == 0) || (i == N))
+	   then return 17.0/48.0;
+	else if ((i == 1) || (i == N-1))
+	   then return 59.0/48.0;
+	else if ((i == 2) || (i == N-2))
+	   then return 43.0/48.0;
+	else if ((i == 3) || (i == N-3))
+	   then return 49.0/48.0;
+	else	return 1.0;
+
+default:
+	error_exit(ERROR_EXIT,
+"***** patch::integration_coeff(): unknown method=(int)%d!\n"
+"                                  (this should never happen!)\n"
+,
+		   int(method));				/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function returns a reference to the ghost zone on a specified
 // edge, after first assert()ing that the edge belongs to this patch.
 //
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 8a90d56..950a7f6 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -252,6 +252,107 @@ public:
 
 
 	//
+	// ***** gridfn operations *****
+	//
+public:
+
+	//
+	// The following enum describes the integration methods supported
+	// by  integrate_gridfn() .
+	//
+	// For convenience of exposition we describe the methods as if for
+	// 1-D integration, but  integrate_gridfn()  actually does 2-D
+	// (surface) integration over the patch.
+	//
+	// Suppose we're computing $\int_{x_0}^{x^N} f(x) \, dx$, using the
+	// equally spaced integration points $f_0$, $f_1$, \dots, $f_N$,
+	// spaced $\Delta x$ apart.  Then the integration methods are as
+	// follows, with the convention that $\langle X \rangle$ denotes
+	// indefinite repetition of the "X" terms, depending on N:
+	//
+	enum	integration_method
+		{
+		// Trapezoid rule
+		// ... character-string name "trapezoid" or "trapezoid rule"
+		// ... 2nd order accurate for smooth functions
+		// ... requires N >= 1
+		// $$
+		// \Delta x \left[
+		//            \half f_0
+		//          + \langle
+		//            f_k
+		//            \rangle
+		//          + \half f_N
+		//          \right]
+		// $$
+		integration_method__trapezoid,
+
+		// Simpson's rule
+		// ... character-string name "Simpson" or "Simpson's rule"
+		// ... 4th order accurate for smooth functions
+		// ... requires N >= 2 and N even
+		// $$
+		// \Delta x \left[
+		//            \frac{1}{3} f_0
+		//          + \frac{4}{3} f_1
+		//          + \langle
+		//            \frac{2}{3} f_{2k} + \frac{4}{3} f_{2k+1}
+		//            \rangle
+		//          + \frac{1}{3} f_N
+		//          \right]
+		// $$
+		integration_method__Simpson,
+
+		// Simpson's rule, variant form
+		// ... characgter-string name "Simpson (variant)"
+		//     or "Simpson's rule (variant)"
+		// ... described in Numerical Recipes 1st edition (4.1.14)
+		// ... 4th order accurate for smooth functions
+		// ... requires N >= 7
+		// $$
+		// \Delta x \left[
+		//            \frac{17}{48} f_0
+		//          + \frac{59}{48} f_1
+		//          + \frac{43}{48} f_2
+		//          + \frac{49}{48} f_3
+		//          + \langle
+		//            f_k
+		//            \rangle
+		//          + \frac{49}{48} f_{N-3}
+		//          + \frac{43}{48} f_{N-2}
+		//          + \frac{59}{48} f_{N-1}
+		//          + \frac{17}{48} f_N
+		//          \right]
+		// $$
+		integration_method__Simpson_variant // no comma here!
+		};
+
+	// decode character string name into internal enum
+	static
+	  enum integration_method
+	    decode_integration_method(const char method_string[]);
+
+
+	// integrate a gridfn: computes an approximation to the surface
+	// integral
+	//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
+	//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
+	// where $dA$ is the area element in $(\rho,sigma)$ coordinates
+	// ... integration method selected by  method  argument
+	//     FIXME: right now this is ignored :( :(
+	fp integrate_gridfn(int src_gfn,
+			    bool area_weighting_flag,
+			    enum integration_method method)
+		const;
+private:
+	// compute integration coefficient $c_i$ where
+	// $\int_{x_0}^{x_N} f(x) \, dx
+	//	\approx \Delta x \, \sum_{i=0}^N c_i f(x_i)$
+	static
+	  fp integration_coeff(enum integration_method method, int N, int i);
+
+
+	//
 	// ***** patch edges ****
 	//
 public:
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 575556f..c8ed7d0 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -23,10 +23,12 @@
 // patch_system::patch_irho_isigma_of_gpn
 // patch_system::ghosted_patch_irho_isigma_of_gpn
 //
+// patch_system::set_gridfn_to_constant
 // patch_system::gridfn_copy
 // patch_system::add_to_ghosted_gridfn
 // patch_system::gridfn_norms
 // patch_system::ghosted_gridfn_norms
+// patch_system::integrate_gridfn
 //
 // patch_system::print_unknown_gridfn
 // patch_system::read_unknown_gridfn
@@ -93,6 +95,9 @@ using jtutil::error_exit;
 //		      be 1; in general
 //			N_overlap_points = 2*N_extend_points + 1
 // delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
+// print_msg_flag = true to print CCTK_VInfo messages describing the
+//			 construction process,
+//		    false to skip this (and be silent except for errors)
 //
 patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 			   enum patch_system_type type_in,
@@ -100,7 +105,8 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 			   fp delta_drho_dsigma,
 			   int min_gfn_in, int max_gfn_in,
 			   int ghosted_min_gfn_in, int ghosted_max_gfn_in,
-			   int interp_handle, int interp_par_table_handle)
+			   int interp_handle, int interp_par_table_handle,
+			   bool print_msg_flag)
 
 	: global_coords_(origin_x_in, origin_y_in, origin_z_in),
 	  type_(type_in),
@@ -124,53 +130,63 @@ switch	(type_in)
 	{
 case full_sphere_patch_system:
 	create_patches(patch_system_info::full_sphere::patch_info_array,
-		       N_ghost_points, N_extend_points,
-		       delta_drho_dsigma);
+		       N_ghost_points, N_extend_points, delta_drho_dsigma,
+		       print_msg_flag);
 	setup_gridfn_storage(min_gfn_in, max_gfn_in,
-			     ghosted_min_gfn_in, ghosted_max_gfn_in);
+			     ghosted_min_gfn_in, ghosted_max_gfn_in,
+			     print_msg_flag);
 	interlink_full_sphere_patch_system(N_overlap_points,
 					   interp_handle,
-					   interp_par_table_handle);
+					   interp_par_table_handle,
+					   print_msg_flag);
 	break;
 case plus_z_hemisphere_patch_system:
 	create_patches(patch_system_info::plus_z_hemisphere::patch_info_array,
-		       N_ghost_points, N_extend_points,
-		       delta_drho_dsigma);
+		       N_ghost_points, N_extend_points, delta_drho_dsigma,
+		       print_msg_flag);
 	setup_gridfn_storage(min_gfn_in, max_gfn_in,
-			     ghosted_min_gfn_in, ghosted_max_gfn_in);
+			     ghosted_min_gfn_in, ghosted_max_gfn_in,
+			     print_msg_flag);
 	interlink_plus_z_hemisphere_patch_system(N_overlap_points,
 						 interp_handle,
-						 interp_par_table_handle);
+						 interp_par_table_handle,
+						 print_msg_flag);
 	break;
 case plus_xy_quadrant_patch_system:
 	create_patches(patch_system_info::plus_xy_quadrant::patch_info_array,
-		       N_ghost_points, N_extend_points,
-		       delta_drho_dsigma);
+		       N_ghost_points, N_extend_points, delta_drho_dsigma,
+		       print_msg_flag);
 	setup_gridfn_storage(min_gfn_in, max_gfn_in,
-			     ghosted_min_gfn_in, ghosted_max_gfn_in);
+			     ghosted_min_gfn_in, ghosted_max_gfn_in,
+			     print_msg_flag);
 	interlink_plus_xy_quadrant_patch_system(N_overlap_points,
 						interp_handle,
-						interp_par_table_handle);
+						interp_par_table_handle,
+						print_msg_flag);
 	break;
 case plus_xz_quadrant_patch_system:
 	create_patches(patch_system_info::plus_xz_quadrant::patch_info_array,
-		       N_ghost_points, N_extend_points,
-		       delta_drho_dsigma);
+		       N_ghost_points, N_extend_points, delta_drho_dsigma,
+		       print_msg_flag);
 	setup_gridfn_storage(min_gfn_in, max_gfn_in,
-			     ghosted_min_gfn_in, ghosted_max_gfn_in);
+			     ghosted_min_gfn_in, ghosted_max_gfn_in,
+			     print_msg_flag);
 	interlink_plus_xz_quadrant_patch_system(N_overlap_points,
 						interp_handle,
-						interp_par_table_handle);
+						interp_par_table_handle,
+						print_msg_flag);
 	break;
 case plus_xyz_octant_patch_system:
 	create_patches(patch_system_info::plus_xyz_octant::patch_info_array,
-		       N_ghost_points, N_extend_points,
-		       delta_drho_dsigma);
+		       N_ghost_points, N_extend_points, delta_drho_dsigma,
+		       print_msg_flag);
 	setup_gridfn_storage(min_gfn_in, max_gfn_in,
-			     ghosted_min_gfn_in, ghosted_max_gfn_in);
+			     ghosted_min_gfn_in, ghosted_max_gfn_in,
+			     print_msg_flag);
 	interlink_plus_xyz_octant_patch_system(N_overlap_points,
 					       interp_handle,
-					       interp_par_table_handle);
+					       interp_par_table_handle,
+					       print_msg_flag);
 	break;
 default:
 	error_exit(ERROR_EXIT,
@@ -213,23 +229,28 @@ delete gridfn_storage_;
 //
 void patch_system::create_patches(const struct patch_info patch_info_in[],
 				  int N_ghost_points, int N_extend_points,
-				  fp delta_drho_dsigma)
+				  fp delta_drho_dsigma,
+				  bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "   constructing %s patch system",
-	   name_of_type(type()));
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "                at origin (%g,%g,%g)",
-	   double(origin_x()), double(origin_y()), double(origin_z()));
+if (print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   constructing %s patch system",
+		   name_of_type(type()));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "                at origin (%g,%g,%g)",
+		   double(origin_x()), double(origin_y()), double(origin_z()));
+	}
 N_grid_points_ = 0;
 ghosted_N_grid_points_ = 0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
 	const struct patch_info& pi = patch_info_in[pn];
 
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "      constructing %s patch",
-		   pi.name);
+	if (print_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "      constructing %s patch",
+			   pi.name);
 
 	struct patch *p;
 	switch	(pi.ctype)
@@ -322,7 +343,8 @@ ghosted_starting_gpn_[N_patches_] = ghosted_N_grid_points_;
 //
 void patch_system::setup_gridfn_storage
 	(int min_gfn_in, int max_gfn_in,
-	 int ghosted_min_gfn_in, int ghosted_max_gfn_in)
+	 int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+	 bool print_msg_flag)
 {
 const int N_gridfns_in = jtutil::how_many_in_range(min_gfn_in, max_gfn_in);
 const int ghosted_N_gridfns_in
@@ -334,14 +356,17 @@ const int ghosted_gfn_stride = ghosted_N_grid_points();
 const int N_storage = gfn_stride * N_gridfns_in;
 const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in;
 
-CCTK_VInfo(CCTK_THORNSTRING, "   setting up gridfn storage");
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      gfn=[%d,%d] ghosted_gfn=[%d,%d]",
-	   min_gfn_in, max_gfn_in,
-	   ghosted_min_gfn_in, ghosted_max_gfn_in);
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      N_grid_points()=%d ghosted_N_grid_points()=%d",
-	   N_grid_points(), ghosted_N_grid_points());
+if (print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING, "   setting up gridfn storage");
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "      gfn=[%d,%d] ghosted_gfn=[%d,%d]",
+		   min_gfn_in, max_gfn_in,
+		   ghosted_min_gfn_in, ghosted_max_gfn_in);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "      N_grid_points()=%d ghosted_N_grid_points()=%d",
+		   N_grid_points(), ghosted_N_grid_points());
+	}
 
 // storage arrays for all gridfns
 gridfn_storage_         = new fp[N_storage];
@@ -370,8 +395,11 @@ ghosted_gridfn_storage_ = new fp[ghosted_N_storage];
 	p.setup_gridfn_storage(gridfn_pars, ghosted_gridfn_pars);
 	}
 
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      checking that storage is partitioned properly");
+if (print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "      checking that storage is partitioned properly");
+	}
 
 // check to make sure storage for distinct gridfns
 // forms a partition of the overall storage array
@@ -508,9 +536,12 @@ const patch& plast  = ith_patch(N_patches()-1);
 //
 void patch_system::interlink_full_sphere_patch_system
 	(int N_overlap_points,
-	 int interp_handle, int interp_par_table_handle)
+	 int interp_handle, int interp_par_table_handle,
+	 bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   interlinking full sphere patch system");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   interlinking full sphere patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -520,7 +551,8 @@ patch& my = ith_patch(patch_number_of_name("-y"));
 patch& mz = ith_patch(patch_number_of_name("-z"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, mx, N_overlap_points);
@@ -535,7 +567,8 @@ create_interpatch_ghost_zones(mz, mx, N_overlap_points);
 create_interpatch_ghost_zones(mz, my, N_overlap_points);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -584,9 +617,12 @@ assert_all_ghost_zones_fully_setup();
 //
 void patch_system::interlink_plus_z_hemisphere_patch_system
 	(int N_overlap_points,
-	 int interp_handle, int interp_par_table_handle)
+	 int interp_handle, int interp_par_table_handle,
+	 bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +z hemisphere patch system");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   interlinking +z hemisphere patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -595,7 +631,8 @@ patch& mx = ith_patch(patch_number_of_name("-x"));
 patch& my = ith_patch(patch_number_of_name("-y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, mx, N_overlap_points);
@@ -610,7 +647,8 @@ mx.create_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge());
 my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge());
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -647,9 +685,12 @@ assert_all_ghost_zones_fully_setup();
 //
 void patch_system::interlink_plus_xy_quadrant_patch_system
 	(int N_overlap_points,
-	 int interp_handle, int interp_par_table_handle)
+	 int interp_handle, int interp_par_table_handle,
+	 bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xy quadrant patch system");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   interlinking +xy quadrant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -657,7 +698,8 @@ patch& py = ith_patch(patch_number_of_name("+y"));
 patch& mz = ith_patch(patch_number_of_name("-z"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(px, py, N_overlap_points);
@@ -674,7 +716,8 @@ create_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(),
 				     true);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -702,9 +745,12 @@ assert_all_ghost_zones_fully_setup();
 //
 void patch_system::interlink_plus_xz_quadrant_patch_system
 	(int N_overlap_points,
-	 int interp_handle, int interp_par_table_handle)
+	 int interp_handle, int interp_par_table_handle,
+	 bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xz quadrant patch system");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   interlinking +xz quadrant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
@@ -712,7 +758,8 @@ patch& py = ith_patch(patch_number_of_name("+y"));
 patch& my = ith_patch(patch_number_of_name("-y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(pz, my, N_overlap_points);
@@ -729,7 +776,8 @@ create_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(),
 				     false);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -757,16 +805,20 @@ assert_all_ghost_zones_fully_setup();
 //
 void patch_system::interlink_plus_xyz_octant_patch_system
 	(int N_overlap_points,
-	 int interp_handle, int interp_par_table_handle)
+	 int interp_handle, int interp_par_table_handle,
+	 bool print_msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   interlinking +xyz octant patch system");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   interlinking +xyz octant patch system");
 
 patch& pz = ith_patch(patch_number_of_name("+z"));
 patch& px = ith_patch(patch_number_of_name("+x"));
 patch& py = ith_patch(patch_number_of_name("+y"));
 
 // create the ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      creating ghost zones");
 create_interpatch_ghost_zones(pz, px, N_overlap_points);
 create_interpatch_ghost_zones(pz, py, N_overlap_points);
 create_interpatch_ghost_zones(px, py, N_overlap_points);
@@ -780,7 +832,8 @@ create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(),
 				     true);
 
 // finish setting up the interpatch ghost zones
-CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      finishing interpatch setup");
 finish_interpatch_setup(pz, px,
 			N_overlap_points,
 			interp_handle, interp_par_table_handle);
@@ -990,7 +1043,8 @@ error_exit(ERROR_EXIT,
 // design would be to return this via a  patch*&  argument, but the design
 // used here seems slightly cleaner to use in practice.)
 //
-patch& patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
+const patch&
+  patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
 	const
 {
 assert( gpn >= 0 );
@@ -1001,7 +1055,7 @@ assert( gpn < N_grid_points() );
 	// 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 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() );
@@ -1034,7 +1088,8 @@ error_exit(PANIC_EXIT,
 // design would be to return this via a  patch*&  argument, but the design
 // used here seems slightly cleaner to use in practice.)
 //
-patch& patch_system::ghosted_patch_irho_isigma_of_gpn
+const patch&
+  patch_system::ghosted_patch_irho_isigma_of_gpn
 	(int gpn, int& irho, int& isigma)
 	const
 {
@@ -1048,7 +1103,7 @@ assert( gpn < ghosted_N_grid_points() );
 	if (    (gpn >= ghosted_starting_gpn_[pn])
 	     && (gpn <  ghosted_starting_gpn_[pn+1]))
 	   then {
-		patch& p = ith_patch(pn);
+		const 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() );
@@ -1071,6 +1126,28 @@ error_exit(PANIC_EXIT,
 //******************************************************************************
 
 //
+// This function sets a (nominal-grid) gridfn to a constant value.
+//
+void patch_system::set_gridfn_to_constant(fp a, int dst_gfn)
+{
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	patch& p = ith_patch(pn);
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		p.gridfn(dst_gfn, irho,isigma) = a;
+		}
+		}
+	}
+}
+
+//******************************************************************************
+
+//
 // This function copies one (nominal-grid) gridfn to another.
 //
 void patch_system::gridfn_copy(int src_gfn, int dst_gfn)
@@ -1120,6 +1197,7 @@ void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn)
 // This function computes norms of a nominal-grid gridfn.
 //
 void patch_system::gridfn_norms(int src_gfn, jtutil::norm<fp>& norms)
+	const
 {
 if (! is_valid_gfn(src_gfn))
    then error_exit(ERROR_EXIT,
@@ -1130,7 +1208,7 @@ norms.reset();
 
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
-	patch& p = ith_patch(pn);
+	const patch& p = ith_patch(pn);
 		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
 		{
 		for (int isigma = p.min_isigma() ;
@@ -1150,6 +1228,7 @@ norms.reset();
 //
 void patch_system::ghosted_gridfn_norms(int ghosted_src_gfn,
 					jtutil::norm<fp>& norms)
+	const
 {
 if (! is_valid_ghosted_gfn(ghosted_src_gfn))
    then error_exit(ERROR_EXIT,
@@ -1159,7 +1238,7 @@ norms.reset();
 
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
-	patch& p = ith_patch(pn);
+	const patch& p = ith_patch(pn);
 		for (int irho = p.ghosted_min_irho() ;
 		     irho <= p.ghosted_max_irho() ;
 		     ++irho)
@@ -1175,6 +1254,35 @@ norms.reset();
 }
 
 //******************************************************************************
+
+//
+// This function computes an approximation to the (surface) integral of
+// a gridfn over the whole patch system,
+//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
+//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
+// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in
+// $(\rho,sigma)$ coordinates.
+//
+// The integration scheme is selected based on the  method  argument.
+//
+// Bugs:
+// The way the integration coefficients are computed is somewhat inefficient.
+//
+fp patch_system::integrate_gridfn(int src_gfn,
+				  bool area_weighting_flag,
+				  enum patch::integration_method method)
+	const
+{
+fp integral = 0.0;
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	const patch& p = ith_patch(pn);
+	integral += p.integrate_gridfn(src_gfn, area_weighting_flag, method);
+	}
+return integral;
+}
+
+//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -1483,7 +1591,7 @@ global_min_ym_ = +INT_MAX;
 global_max_ym_ = -INT_MAX;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
-	patch& p = ith_patch(pn);
+	const patch& p = ith_patch(pn);
 		// n.b. these loops must use _int_ variables for the loop
 		//      to terminate!
 		for (int want_min = false ; want_min <= true ; ++want_min)
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 5688d48..75ca8cc 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -43,6 +43,7 @@ class	patch_system
 	//
 public:
 	// what patch-system type are supported?
+	// (see "patch_system_info.hh" for detailed descriptions of these)
 	enum patch_system_type {
 			       full_sphere_patch_system,
 			       plus_z_hemisphere_patch_system,
@@ -52,7 +53,8 @@ public:
 			       };
 
 	// maximum number of patches in any patch-system type
-	static const int max_N_patches = 6;
+	static
+	  const int max_N_patches = 6;
 
 	// decode patch system type into N_patches
 	static
@@ -102,7 +104,9 @@ public:
 	int N_patches() const { return N_patches_; }
 
 	// get patches by patch number
-	patch &ith_patch(int pn) const
+	const patch &ith_patch(int pn) const
+		{ return * all_patches_[pn]; }
+	patch &ith_patch(int pn)
 		{ return * all_patches_[pn]; }
 
 	// find a patch by name, return patch number; error_exit() if not found
@@ -241,6 +245,9 @@ private:
 	// ***** gridfn operations *****
 	//
 public:
+	// dst = a
+	void set_gridfn_to_constant(fp a, int dst_gfn);
+
 	// dst = src
 	void gridfn_copy(int src_gfn, int dst_gfn);
 
@@ -248,8 +255,21 @@ public:
 	void add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn);
 
 	// compute norms of gridfn
-	void         gridfn_norms(int         src_gfn, jtutil::norm<fp>& norms);
-	void ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms);
+	void         gridfn_norms(int         src_gfn, jtutil::norm<fp>& norms)
+		const;
+	void ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms)
+		const;
+
+	// integrate a gridfn: computes an approximation to the surface
+	// integral over the whole patch system
+	//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
+	//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
+	// where $dA$ is the area element in $(\rho,sigma)$ coordinates
+	// ... integration method selected by  method  argument
+	fp integrate_gridfn(int src_gfn,
+			    bool area_weighting_flag,
+			    enum patch::integration_method method)
+		const;
 
 
 	//
@@ -353,9 +373,11 @@ public:
 		}
 	// ... 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&
+	  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 patch&
+	  ghosted_patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
 		const;
 
 	// access actual gridfn data arrays
@@ -387,7 +409,8 @@ public:
 		     fp delta_drho_dsigma,
 		     int min_gfn_in, int max_gfn_in,
 		     int ghosted_min_gfn_in, int ghosted_max_gfn_in,
-		     int interp_handle_in, int interp_par_table_handle_in);
+		     int interp_handle_in, int interp_par_table_handle_in,
+		     bool print_msg_flag);
 	~patch_system();
 
 
@@ -401,29 +424,36 @@ private:
 	// does *NOT* set up gridfns
 	void create_patches(const struct patch_info patch_info_in[],
 			    int N_ghost_points, int N_extend_points,
-			    fp delta_drho_dsigma);
+			    fp delta_drho_dsigma,
+			    bool print_msg_flag);
 
 	// setup all gridfns with contiguous-across-patches storage
 	void setup_gridfn_storage
 		(int min_gfn_in, int max_gfn_in,
-		 int ghosted_min_gfn_in, int ghosted_max_gfn_in);
+		 int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+		 bool print_msg_flag);
 
 	// create/interlink ghost zones
 	void interlink_full_sphere_patch_system
 		(int N_overlap_points,
-		 int interp_handle, int interp_par_table_handle);
+		 int interp_handle, int interp_par_table_handle,
+		 bool print_msg_flag);
 	void interlink_plus_z_hemisphere_patch_system
 		(int N_overlap_points,
-		 int interp_handle, int interp_par_table_handle);
+		 int interp_handle, int interp_par_table_handle,
+		 bool print_msg_flag);
 	void interlink_plus_xy_quadrant_patch_system
 		(int N_overlap_points,
-		 int interp_handle, int interp_par_table_handle);
+		 int interp_handle, int interp_par_table_handle,
+		 bool print_msg_flag);
 	void interlink_plus_xz_quadrant_patch_system
 		(int N_overlap_points,
-		 int interp_handle, int interp_par_table_handle);
+		 int interp_handle, int interp_par_table_handle,
+		 bool print_msg_flag);
 	void interlink_plus_xyz_octant_patch_system
 		(int N_overlap_points,
-		 int interp_handle, int interp_par_table_handle);
+		 int interp_handle, int interp_par_table_handle,
+		 bool print_msg_flag);
 
 	// create/interlink a pair of periodic-symmetry ghost zones
 	static