change $Id:$ --> $Header:$

rename src/util/ --> src/patch/
add Jacobian code for computing surface integrals
add support for choosing ASCII or HDF5 data files (HDF5 not implemented yet)


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

51 files changed, 312 insertions, 184 deletions

diff --git a/src/README b/src/README
index 5ff17b4..0ee1180 100644
--- a/src/README
+++ b/src/README
@@ -4,11 +4,12 @@ See ../doc/ for further documentation.
 The main subdirectories under this directory are:
 driver/		high-level driver routines to solve the H(h) = 0
 		equations and interface to the rest of Cactus
+elliptic/	code to solve elliptic equations on the multipatch $S^2$
 gr/		relativity code; all knowledge of the actual apparent
 		horizon equation lives in the code in this directory
 gr.cg/		Maple-generated code to compute the H(h) function
 		and its Jacobian coefficients
-util/		basic multipatch infrastructure for storing and
+patch/		basic multipatch infrastructure for storing and
 		finite differencing gridfns in angular coordinates
 		on the surface of a 2-sphere
 jtutil/		low-level utility library in C++ namespace jtutil::
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index e7ce466..7267088 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -1,5 +1,5 @@
 // Newton.cc -- solve H(h) = 0 via Newton's method
-// $Id$
+// $Header$
 //
 // Newton_solve - driver to solve H(h) = 0 via Newton's method
 //
@@ -21,14 +21,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -67,7 +67,7 @@ bool Newton_solve(patch_system& ps,
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "Newton iteration %d", iteration);
 
-	if (solver_info.output_h_and_H_at_each_Newton_iteration)
+	if (IO_info.output_h_and_H_at_each_Newton_iteration)
 	   then output_gridfn(ps, gfns::gfn__h,
 			      IO_info, IO_info.h_base_file_name,
 			      hn, verbose_info.print_algorithm_details,
@@ -86,7 +86,7 @@ bool Newton_solve(patch_system& ps,
 			   iteration,
 			   H_norms.rms_norm(), H_norms.infinity_norm());
 
-	if (solver_info.output_h_and_H_at_each_Newton_iteration)
+	if (IO_info.output_h_and_H_at_each_Newton_iteration)
 	   then output_gridfn(ps, gfns::gfn__H,
 			      IO_info, IO_info.H_base_file_name,
 			      hn, verbose_info.print_algorithm_details,
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index f01a64c..38efffc 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -38,6 +38,16 @@ enum	Jacobian_method
 	Jacobian_method__symbolic_diff // no comma
 	};
 
+//
+// this enum holds the (a) decoded  file_format  parameter,
+// i.e. it specifies what format of input/output file(s) we should use
+//
+enum	file_format
+	{
+	file_format__ASCII,
+	file_format__HDF5 // no comma
+	};
+
 //******************************************************************************
 
 //
@@ -56,7 +66,6 @@ struct	Jacobian_info
 struct	solver_info
 	{
 	int max_Newton_iterations;
-	bool output_h_and_H_at_each_Newton_iteration;
 	fp max_Delta_h_over_h;
 	fp H_norm_for_convergence;
 	fp Delta_h_norm_for_convergence;
@@ -68,6 +77,12 @@ struct	solver_info
 //
 struct	IO_info
 	{
+	bool output_initial_guess;
+	bool output_h_and_H_at_each_Newton_iteration;
+	bool output_h, output_H;
+	enum file_format file_format;
+	const char* ASCII_file_name_extension;
+	const char* HDF5_file_name_extension;
 	const char* h_base_file_name;
 	const char* H_base_file_name;
 	const char* Jacobian_base_file_name;
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 3f3279e..0739a3c 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -4,6 +4,9 @@
 // <<<access to persistent data>>>
 // <<<prototypes for functions local to this file>>>
 // AHFinderDirect_find_horizons - top-level driver to find apparent horizons
+// find_horizon - find a horizon
+// BH_diagnostics - compute BH diagnostics for a horizon
+// surface_integral - compute surface integral of a gridfn over the horizon
 //
 
 #include <stdio.h>
@@ -24,14 +27,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -64,6 +67,8 @@ bool find_horizon(enum method method,
 void BH_diagnostics(enum patch::integration_method surface_integral_method,
 		    const struct verbose_info& verbose_info,
 		    struct AH_info& AH_info);
+fp surface_integral(const patch_system& ps, int src_gfn,
+		    enum patch::integration_method method);
 	  };
 
 //******************************************************************************
@@ -105,7 +110,7 @@ if (state.timer_handle >= 0)
 
 		if (verbose_info.print_physics_details)
 		   then CCTK_VInfo(CCTK_THORNSTRING,
-				   "AH found: A=%g m=%g at (%g,%g,%g)",
+				   "AH found: A=%.10g m=%.10g at (%g,%g,%g)",
 				   double(AH_info.area), double(AH_info.mass),
 				   double(AH_info.centroid_x),
 				   double(AH_info.centroid_y),
@@ -172,9 +177,10 @@ case method__horizon_function:
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "   H(h) rms-norm %.2e, infinity-norm %.2e",
 			   H_norms.rms_norm(), H_norms.infinity_norm());
-	output_gridfn(ps, gfns::gfn__H,
-		      IO_info, IO_info.H_base_file_name,
-		      hn, true);
+	if (IO_info.output_H)
+	   then output_gridfn(ps, gfns::gfn__H,
+			      IO_info, IO_info.H_base_file_name,
+			      hn, true);
 	return true;					// *** NORMAL RETURN ***
 	  }
 
@@ -249,12 +255,14 @@ case method__Newton_solve:
 	if (! status)
 	   then return false;				// *** ERROR RETURN ***
 
-	output_gridfn(ps, gfns::gfn__h,
-		      IO_info, IO_info.h_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-	output_gridfn(ps, gfns::gfn__H,
-		      IO_info, IO_info.H_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
+	if (IO_info.output_h)
+	   then output_gridfn(ps, gfns::gfn__h,
+			      IO_info, IO_info.h_base_file_name,
+			      hn, verbose_info.print_algorithm_details);
+	if (IO_info.output_H)
+	   then output_gridfn(ps, gfns::gfn__H,
+			      IO_info, IO_info.H_base_file_name,
+			      hn, verbose_info.print_algorithm_details);
 	return true;					// *** NORMAL RETURN ***
 	  }
 
@@ -292,18 +300,14 @@ const patch_system& ps = * AH_info.ps_ptr;
 //
 // compute raw surface integrals
 //
-fp integral_one = ps.integrate_gridfn(gfns::gfn__one,
-				      true,	// area weighting
-				      surface_integral_method);
-fp integral_x = ps.integrate_gridfn(gfns::gfn__global_x,
-				    true,		// area weighting
-				    surface_integral_method);
-fp integral_y = ps.integrate_gridfn(gfns::gfn__global_y,
-				    true,		// area weighting
-				    surface_integral_method);
-fp integral_z = ps.integrate_gridfn(gfns::gfn__global_z,
-				    true,		// area weighting
-				    surface_integral_method);
+fp integral_one = surface_integral(ps, gfns::gfn__one,
+				   surface_integral_method);
+fp integral_x = surface_integral(ps, gfns::gfn__global_x,
+				 surface_integral_method);
+fp integral_y = surface_integral(ps, gfns::gfn__global_y,
+				 surface_integral_method);
+fp integral_z = surface_integral(ps, gfns::gfn__global_z,
+				 surface_integral_method);
 
 //
 // adjust integrals to take into account patch system not covering full sphere
@@ -353,3 +357,23 @@ AH_info.centroid_y = integral_y / integral_one;
 AH_info.centroid_z = integral_z / integral_one;
 }
 	  }
+
+//******************************************************************************
+
+//
+// This function computes the surface integral of a gridfn over the
+// horizon.
+//
+namespace {
+fp surface_integral(const patch_system& ps, int src_gfn,
+		    enum patch::integration_method method)
+{
+return ps.integrate_gridfn
+	   (src_gfn,
+	    gfns::gfn__h,
+	    gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+				gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+						    gfns::gfn__g_dd_33,
+	    method);
+}
+	  }
diff --git a/src/driver/horizon_Jacobian.cc b/src/driver/horizon_Jacobian.cc
index 4cc751d..ec50d78 100644
--- a/src/driver/horizon_Jacobian.cc
+++ b/src/driver/horizon_Jacobian.cc
@@ -1,5 +1,5 @@
 // horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
-// $Id$
+// $Header$
 //
 // <<<prototypes for functions local to this file>>>
 //
@@ -28,14 +28,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
index 94d81d2..03aef25 100644
--- a/src/driver/initial_guess.cc
+++ b/src/driver/initial_guess.cc
@@ -26,14 +26,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -140,7 +140,7 @@ IO_info.time_iteration = cctk_iteration;
 			   initial_guess_method);	/*NOTREACHED*/
 
 	// write initial guess back to the data file?
-	if (output_initial_guess != 0)
+	if (IO_info.output_initial_guess)
 	   then output_gridfn(ps, gfns::gfn__h,
 			      IO_info, IO_info.h_base_file_name,
 			      hn, verbose_info.print_algorithm_highlights);
diff --git a/src/driver/io.cc b/src/driver/io.cc
index 6feead1..30f15b5 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -27,14 +27,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -78,9 +78,26 @@ if (print_msg_flag)
 			   "   reading \"%s\"", file_name);
 	}
 
-ps.read_ghosted_gridfn(unknown_gfn,
-		       file_name,
-		       false);			// no ghost zones in data file
+switch	(IO_info.file_format)
+	{
+case file_format__ASCII:
+	ps.read_ghosted_gridfn(unknown_gfn,
+			       file_name,
+			       false);		// no ghost zones in data file
+	break;
+
+case file_format__HDF5:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"input_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/
+
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   input_gridfn(): unknown IO_info.file_format=(int)%d!\n"
+"                   (this should never happen!)"
+		   ,
+		   int(IO_info.file_format));			/*NOTREACHED*/
+	}
 }
 
 //******************************************************************************
@@ -99,32 +116,49 @@ void output_gridfn(patch_system& ps, int unknown_gfn,
 {
 const char* file_name = io_file_name(IO_info, base_file_name,
 				     hn, AHF_iteration);
-if (print_msg_flag)
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   writing \"%s\"", file_name);
 
-switch	(unknown_gfn)
+switch	(IO_info.file_format)
 	{
-case gfns::gfn__h:
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   writing h to \"%s\"", file_name);
-	ps.print_ghosted_gridfn_with_xyz(unknown_gfn,
-					 true, gfns::gfn__h,
-					 file_name,
-					 false);	// no ghost zones
-	break;
-case gfns::gfn__H:
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   writing H(h) to \"%s\"", file_name);
-	ps.print_gridfn(unknown_gfn,
-			file_name);
+case file_format__ASCII:
+	switch	(unknown_gfn)
+		{
+	case gfns::gfn__h:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "   writing h to \"%s\"", file_name);
+		ps.print_ghosted_gridfn_with_xyz(unknown_gfn,
+						 true, gfns::gfn__h,
+						 file_name,
+						 false);      // no ghost zones
+		break;
+	case gfns::gfn__H:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "   writing H to \"%s\"", file_name);
+		ps.print_gridfn(unknown_gfn,
+				file_name);
+		break;
+	default:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "   writing \"%s\"", file_name);
+		ps.print_gridfn(unknown_gfn,
+				file_name);
+		break;
+		}
 	break;
+
+case file_format__HDF5:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"output_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/
+
 default:
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   writing \"%s\"", file_name);
-	ps.print_gridfn(unknown_gfn,
-			file_name);
-	break;
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   output_gridfn(): unknown IO_info.file_format=(int)%d!\n"
+"                    (this should never happen!)"
+		   ,
+		   int(IO_info.file_format));			/*NOTREACHED*/
 	}
 }
 
@@ -134,6 +168,9 @@ default:
 // This function prints one or two Jacobian matrices (and their difference
 // in the latter case) to a named output file.
 //
+// Bugs:
+// - The file format is hardwired to ASCII.
+//
 // Arguments:
 // A null Jacobian pointer means to skip that Jacobian.
 //
@@ -277,13 +314,35 @@ const char* io_file_name(struct IO_info& IO_info, const char base_file_name[],
 const int N_file_name_buffer = 200;
 static char file_name_buffer[N_file_name_buffer];
 
+const char* file_name_extension;
+switch	(IO_info.file_format)
+	{
+case file_format__ASCII:
+	file_name_extension = IO_info.ASCII_file_name_extension;
+	break;
+case file_format__HDF5:
+	file_name_extension = IO_info.HDF5_file_name_extension;
+	break;
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   io_file_name(): unknown IO_info.file_format=(int)%d!\n"
+"                   (this should never happen!)"
+		   ,
+		   int(IO_info.file_format));			/*NOTREACHED*/
+	}
+
 if (AHF_iteration == 0)
    then snprintf(file_name_buffer, N_file_name_buffer,
-		 "%s.t%d.ah%d.dat",
-		 base_file_name, IO_info.time_iteration, hn);
+		 "%s.t%d.ah%d.%s",
+		 base_file_name,
+		 IO_info.time_iteration, hn,
+		 file_name_extension);
    else snprintf(file_name_buffer, N_file_name_buffer,
-		 "%s.t%d.ah%d.it%d.dat",
-		 base_file_name, IO_info.time_iteration, hn, AHF_iteration);
+		 "%s.t%d.ah%d.it%d.%s",
+		 base_file_name,
+		 IO_info.time_iteration, hn, AHF_iteration,
+		 file_name_extension);
 
 return file_name_buffer;
 }
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 86ca722..87fcd86 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -8,6 +8,7 @@
 ///
 /// decode_method - decode the method parameter
 /// decode_verbose_level - decode the verbose_level parameter
+/// decode_file_format - decode the file_format parameter
 /// decode_Jacobian_method - decode the Jacobian_method parameter
 /// decode_geometry_method - decode the geometry_method parameter
 ///
@@ -32,14 +33,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -59,6 +60,8 @@ enum method
   decode_method(const char method_string[]);
 enum verbose_level
   decode_verbose_level(const char verbose_level_string[]);
+enum file_format
+  decode_file_format(const char file_format_string[]);
 enum Jacobian_method
   decode_Jacobian_method(const char Jacobian_method_string[]);
 enum geometry_method
@@ -122,6 +125,14 @@ state.surface_integral_method
 const struct verbose_info& verbose_info = state.verbose_info;
 
 struct IO_info& IO_info = state.IO_info;
+IO_info.output_initial_guess     = (output_initial_guess != 0);
+IO_info.output_h_and_H_at_each_Newton_iteration
+			      = (output_h_and_H_at_each_Newton_iteration != 0);
+IO_info.output_h = (output_h != 0);
+IO_info.output_H = (output_H_of_h != 0);
+IO_info.file_format = decode_file_format(file_format);
+IO_info.ASCII_file_name_extension = ASCII_file_name_extension;
+IO_info.HDF5_file_name_extension  = HDF5_file_name_extension;
 IO_info.h_base_file_name        = h_base_file_name;
 IO_info.H_base_file_name        = H_of_h_base_file_name;
 IO_info.Jacobian_base_file_name = Jacobian_base_file_name;
@@ -135,8 +146,6 @@ Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
 
 struct solver_info& solver_info = state.solver_info;
 solver_info.max_Newton_iterations = max_Newton_iterations;
-solver_info.output_h_and_H_at_each_Newton_iteration
-	= (output_h_and_H_at_each_Newton_iteration != 0);
 solver_info.max_Delta_h_over_h           = max_Delta_h_over_h;
 solver_info.H_norm_for_convergence       = H_norm_for_convergence;
 solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
@@ -339,6 +348,26 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 //******************************************************************************
 
 //
+// This function decodes the  file_format  parameter (string) into an
+// internal enum for future use.
+//
+namespace {
+enum file_format
+  decode_file_format(const char file_format_string[])
+{
+if	(STRING_EQUAL(file_format_string, "ASCII"))
+   then return file_format__ASCII;
+else if (STRING_EQUAL(file_format_string, "HDF5"))
+   then return file_format__HDF5;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"decode_file_format(): unknown file_format_string=\"%s\"!",
+		   file_format_string);				/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+
+//
 // This function decodes the  Jacobian_method  parameter (string) into
 // an internal enum for future use.
 //
diff --git a/src/driver/state.cc b/src/driver/state.cc
index b1f78c7..4c6c1b8 100644
--- a/src/driver/state.cc
+++ b/src/driver/state.cc
@@ -19,14 +19,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc
index 4a70b9e..353c625 100644
--- a/src/elliptic/Jacobian.cc
+++ b/src/elliptic/Jacobian.cc
@@ -1,5 +1,5 @@
 // Jacobian.cc -- data structures for the Jacobian matrix
-// $Id$
+// $Header$
 
 //
 // decode_Jacobian_type
@@ -38,14 +38,14 @@ using std::FILE;
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "Jacobian.hh"
 // FIXME: Cactus's CCTK_FCALL() isn't expanded in .h files (this is a bug),
@@ -54,7 +54,7 @@ using jtutil::error_exit;
 
 //***** begin "lapack.h" contents ******
 /* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */
-/* $Id$ */
+/* $Header$ */
 
 /*
  * prerequisites:
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index a62387c..2fcdcc6 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -1,5 +1,5 @@
 // Jacobian.hh -- data structures for the Jacobian matrix
-// $Id$
+// $Header$
 
 //
 // Jacobian -- abstract base class to describe a Jacobian matrix
diff --git a/src/elliptic/README b/src/elliptic/README
index c1e009d..e20a5ba 100644
--- a/src/elliptic/README
+++ b/src/elliptic/README
@@ -1,5 +1,5 @@
 This directory contains code for solving elliptic systems on our
-($S^2$ topology) multipatch system in ../util/.
+($S^2$ topology) multipatch system in ../patch/.
 
 The main files are as follows:
 
diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc
index 628e641..3d87606 100644
--- a/src/gr/Schwarzschild_EF.cc
+++ b/src/gr/Schwarzschild_EF.cc
@@ -1,5 +1,5 @@
 // Schwarzschild_EF.cc -- set up Schwarzschild/EF geometry variables
-// $Id$
+// $Header$
 //
 // <<<prototypes for functions local to this file>>>
 // Schwarzschild_EF_geometry - top-level driver
@@ -27,14 +27,14 @@
 #include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
@@ -117,7 +117,7 @@ const fp Delta_xyz = gi.geometry__Schwarzschild_EF__Delta_xyz;
 if (msg_flag)
    then {
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "      setting up Schwarzschild/EF geometry");
+		   "      setting up exact Schwarzschild/EF geometry");
 	CCTK_VInfo(CCTK_THORNSTRING,
 		   "                 posn=(%g,%g,%g) mass=%g",
 		   double(x_posn),double(y_posn),double(z_posn), double(mass));
diff --git a/src/gr/auxiliary.maple b/src/gr/auxiliary.maple
index 99069e7..9ae4838 100644
--- a/src/gr/auxiliary.maple
+++ b/src/gr/auxiliary.maple
@@ -1,5 +1,5 @@
 # auxiliary.maple -- Maple code to compute ADM auxiliary quantities
-# $Id$
+# $Header$
 
 #
 # auxiliary - compute and optionally generate C code for auxiliary quantities
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index b63ff2a..02bf370 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -1,5 +1,5 @@
 // cg.hh -- interface between machine-generated and other code
-// $Id$
+// $Header$
 
 //
 // This header file defines the "virtual machine" used by machine-generated
diff --git a/src/gr/curvature.maple b/src/gr/curvature.maple
index 2652a5c..3d58939 100644
--- a/src/gr/curvature.maple
+++ b/src/gr/curvature.maple
@@ -1,5 +1,5 @@
 # curvature.maple - evaluate various "derivatives of the metric" quantities
-# $Id$
+# $Header$
 #
 # curvature - overall driver
 ## inverse_metric_gradient - compute partial_d_g_uu
diff --git a/src/gr/doit.maple b/src/gr/doit.maple
index 3721186..5055e3b 100644
--- a/src/gr/doit.maple
+++ b/src/gr/doit.maple
@@ -1,5 +1,5 @@
 # top-level Maple file to read/run all code in this directory
-# $Id$
+# $Header$
 
 read "../maple/setup.mm";
 
diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh
index a343974..5404f46 100644
--- a/src/gr/gfns.hh
+++ b/src/gr/gfns.hh
@@ -1,5 +1,5 @@
 // gfns.hh -- define gfns of all gridfns
-// $Id$
+// $Header$
 
 namespace gfns
 	  {
diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
index 5a7ebc1..7579e4d 100644
--- a/src/gr/gr_gfas.minc
+++ b/src/gr/gr_gfas.minc
@@ -1,5 +1,5 @@
 # Maple include file listing all gfa functional-dependence forms
-# $Id$
+# $Header$
 
 g_dd,
 K_dd,
diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple
index f7c191c..e0213a7 100644
--- a/src/gr/horizon.maple
+++ b/src/gr/horizon.maple
@@ -1,5 +1,5 @@
 # horizon.maple - evaluate $H$ = LHS of horizon function
-# $Id$
+# $Header$
 #
 # horizon - overall driver
 ## non_unit_normal - compute s_d
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 3bbed78..17a312d 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -1,5 +1,5 @@
 // horizon_function.cc -- evaluate LHS function H(h)
-// $Id$
+// $Header$
 //
 // <<<prototypes for functions local to this file>>>
 // horizon_function - top-level driver
@@ -27,14 +27,14 @@ using jtutil::error_exit;
 using jtutil::pow2;
 using jtutil::pow4;
 
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
 
 #include "../elliptic/Jacobian.hh"
 
diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple
index e833338..ff39681 100644
--- a/src/gr/setup_gr_gfas.maple
+++ b/src/gr/setup_gr_gfas.maple
@@ -1,5 +1,5 @@
 # Maple code to set up all gridfn arrays
-# $Id$
+# $Header$
 
 #
 # setup_gr_gfas - setup all the GR gfas
diff --git a/src/include/config.hh b/src/include/config.hh
index 7afbc32..b5a61bd 100644
--- a/src/include/config.hh
+++ b/src/include/config.hh
@@ -1,5 +1,5 @@
 // config.hh -- compile-time configuration for AHFinderDirect
-// $Id$
+// $Header$
 
 //
 // prerequisites:
diff --git a/src/jtutil/array.cc b/src/jtutil/array.cc
index e48ae27..b660ced 100644
--- a/src/jtutil/array.cc
+++ b/src/jtutil/array.cc
@@ -1,5 +1,5 @@
 // array.cc -- array template classes
-// $Id$
+// $Header$
 //
 // jtutil::array1d::array1d - 1D array template constructor
 // jtutil::array1d::~array1d - 1D array template destructor
diff --git a/src/jtutil/array.hh b/src/jtutil/array.hh
index 9f4715f..51227a6 100644
--- a/src/jtutil/array.hh
+++ b/src/jtutil/array.hh
@@ -1,5 +1,5 @@
 // array.hh -- array template classes
-// $Id$
+// $Header$
 //
 // array1d - 1D array template
 // array2d - 2D array template
diff --git a/src/jtutil/cpm_map.cc b/src/jtutil/cpm_map.cc
index 7b1a897..db04d75 100644
--- a/src/jtutil/cpm_map.cc
+++ b/src/jtutil/cpm_map.cc
@@ -1,5 +1,5 @@
 // cpm_map.cc -- "integer +/-" mapping  i --> j = const +/- i
-// $Id$
+// $Header$
 //
 // jtutil::cpm_map::cpm_map  # mirror map, specified by fixed point
 // jtutil::cpm_map::cpm_map  # generic map specified by sample point & sign
diff --git a/src/jtutil/cpm_map.hh b/src/jtutil/cpm_map.hh
index dacd715..3afce11 100644
--- a/src/jtutil/cpm_map.hh
+++ b/src/jtutil/cpm_map.hh
@@ -1,5 +1,5 @@
 // cpm_map.hh -- "integer +/-" mapping  i --> j = const +/- i
-// $Id$
+// $Header$
 //
 // cpm_map - "integer +/-" mapping
 //
diff --git a/src/jtutil/fuzzy.cc b/src/jtutil/fuzzy.cc
index 9652618..e6b598d 100644
--- a/src/jtutil/fuzzy.cc
+++ b/src/jtutil/fuzzy.cc
@@ -1,5 +1,5 @@
 // fuzzy.cc -- template for fuzzy comparisons et al on floating point values
-// $Id$
+// $Header$
 //
 // jtutil::fuzzy::tolerance - comparison tolerance
 // jtutil::fuzzy::EQ
diff --git a/src/jtutil/linear_map.cc b/src/jtutil/linear_map.cc
index 55ec8af..514f636 100644
--- a/src/jtutil/linear_map.cc
+++ b/src/jtutil/linear_map.cc
@@ -1,5 +1,5 @@
 // linear_map.cc -- linear mapping from integers <--> floating point values
-// $Id$
+// $Header$
 //
 // jtutil::linear_map::linear_map - basic ctor
 // jtutil::linear_map::linear_map - ctor for subrange of existing linear_map
diff --git a/src/jtutil/linear_map.hh b/src/jtutil/linear_map.hh
index 59cdbdf..d37ec75 100644
--- a/src/jtutil/linear_map.hh
+++ b/src/jtutil/linear_map.hh
@@ -1,5 +1,5 @@
 // linear_map.hh -- linear mapping from integers <--> floating point values
-// $Id$
+// $Header$
 //
 // jtutil::linear_map - linear mapping from integers <--> floating point values
 //
diff --git a/src/jtutil/norm.cc b/src/jtutil/norm.cc
index 2e20327..1e6c1b5 100644
--- a/src/jtutil/norm.cc
+++ b/src/jtutil/norm.cc
@@ -1,5 +1,5 @@
 // norm.cc -- compute norms
-// $Id$
+// $Header$
 
 //
 // *** class jtutil::norm::
diff --git a/src/jtutil/round.cc b/src/jtutil/round.cc
index 06309b6..c0b094d 100644
--- a/src/jtutil/round.cc
+++ b/src/jtutil/round.cc
@@ -1,5 +1,5 @@
 // round.hh -- template for rounding floating point values
-// $Id$
+// $Header$
 //
 // *** Implementation Notes ***
 // jtutil::round::to_integer
diff --git a/src/jtutil/test_array.cc b/src/jtutil/test_array.cc
index 1db587a..1f4a001 100644
--- a/src/jtutil/test_array.cc
+++ b/src/jtutil/test_array.cc
@@ -1,5 +1,5 @@
 // test_array.cc -- test driver for array.hh classes
-// $Id$
+// $Header$
 
 #include <assert.h>
 #include <stdio.h>
diff --git a/src/jtutil/test_array2.cc b/src/jtutil/test_array2.cc
index 2e99ae1..6410c52 100644
--- a/src/jtutil/test_array2.cc
+++ b/src/jtutil/test_array2.cc
@@ -1,5 +1,5 @@
 // test_array2.cc -- test driver for array.hh classes
-// $Id$
+// $Header$
 
 #include <assert.h>
 #include <math.h>
diff --git a/src/jtutil/test_cpm_map.cc b/src/jtutil/test_cpm_map.cc
index e9098e3..8b5b182 100644
--- a/src/jtutil/test_cpm_map.cc
+++ b/src/jtutil/test_cpm_map.cc
@@ -1,5 +1,5 @@
 // test_cpm_map.cc -- test driver for cpm_map class
-// $Id$
+// $Header$
 
 #include <assert.h>
 #include <stdio.h>
diff --git a/src/jtutil/test_fuzzy.cc b/src/jtutil/test_fuzzy.cc
index 8620a07..3384d23 100644
--- a/src/jtutil/test_fuzzy.cc
+++ b/src/jtutil/test_fuzzy.cc
@@ -1,5 +1,5 @@
 // test_fuzzy.cc -- test driver for fuzzy<float> and fuzzy<double>
-// $Id$
+// $Header$
 //
 // main
 // check<fp>
diff --git a/src/jtutil/test_linear_map.cc b/src/jtutil/test_linear_map.cc
index a35f40f..e756a7a 100644
--- a/src/jtutil/test_linear_map.cc
+++ b/src/jtutil/test_linear_map.cc
@@ -1,5 +1,5 @@
 // test_linear_map.cc -- test driver for linear_map class
-// $Id$
+// $Header$
 
 #include <assert.h>
 #include <stdio.h>
diff --git a/src/jtutil/test_modulo.cc b/src/jtutil/test_modulo.cc
index a4b2198..322ddd7 100644
--- a/src/jtutil/test_modulo.cc
+++ b/src/jtutil/test_modulo.cc
@@ -1,5 +1,5 @@
 // test_modulo -- quick-n-dirty test of modulo_reduce()
-// $Id$
+// $Header$
 
 #include <stdio.h>
 #include <assert.h>
diff --git a/src/jtutil/test_norm.cc b/src/jtutil/test_norm.cc
index fb24344..6f89030 100644
--- a/src/jtutil/test_norm.cc
+++ b/src/jtutil/test_norm.cc
@@ -1,5 +1,5 @@
 // test_norm -- quick-n-dirty test of norm template class
-// $Id$
+// $Header$
 
 #include <cstdio>
 #include <cassert>
diff --git a/src/jtutil/test_round.cc b/src/jtutil/test_round.cc
index 630cc7d..faf4942 100644
--- a/src/jtutil/test_round.cc
+++ b/src/jtutil/test_round.cc
@@ -1,5 +1,5 @@
 // test_round.cc -- test driver for round<fp> and round<double>
-// $Id$
+// $Header$
 
 #include <stdio.h>
 #include <assert.h>
diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh
index 3d91422..fa46f70 100644
--- a/src/jtutil/util.hh
+++ b/src/jtutil/util.hh
@@ -1,5 +1,5 @@
 // util.hh -- stuff for my C++ utility library
-// $Id$
+// $Header$
 
 //
 // prerequisites:
diff --git a/src/make.code.defn b/src/make.code.defn
index ce8faf8..f9e6261 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -5,4 +5,4 @@
 SRCS = 
 
 # Subdirectories containing source files
-SUBDIRS = jtutil util elliptic gr driver
+SUBDIRS = jtutil patch elliptic gr driver
diff --git a/src/maple/Diff.maple b/src/maple/Diff.maple
index 820ac86..10a66ee 100644
--- a/src/maple/Diff.maple
+++ b/src/maple/Diff.maple
@@ -1,5 +1,5 @@
 # Diff.maple -- Diff() and friends
-# $Id$
+# $Header$
 
 #
 # simplify/Diff - simplify(...) function for expressions containing Diff()
diff --git a/src/maple/codegen2.maple b/src/maple/codegen2.maple
index 3f459c4..eaa881b 100644
--- a/src/maple/codegen2.maple
+++ b/src/maple/codegen2.maple
@@ -1,5 +1,5 @@
 # codegen2.maple -- generate C code from Maple expressions
-# $Id$
+# $Header$
 #
 # codegen2 - generate C code from Maple expressions
 #
diff --git a/src/maple/coeffs.maple b/src/maple/coeffs.maple
index e2b90e3..059d8a7 100644
--- a/src/maple/coeffs.maple
+++ b/src/maple/coeffs.maple
@@ -1,5 +1,5 @@
 # coeffs.maple - set up coordinate coefficients
-# $Id$
+# $Header$
 #
 # setup_coeff_gfas - create the coordinate-transformation coefficient gfas
 #
diff --git a/src/maple/coords.maple b/src/maple/coords.maple
index d76993c..abba0e9 100644
--- a/src/maple/coords.maple
+++ b/src/maple/coords.maple
@@ -1,5 +1,5 @@
 # coords.maple -- set up and covert between coordinates, gridfn arrays, etc
-# $Id$
+# $Header$
 
 #
 # setup_coords - set up coordinates
diff --git a/src/maple/coords.minc b/src/maple/coords.minc
index 0924b94..b72a25e 100644
--- a/src/maple/coords.minc
+++ b/src/maple/coords.minc
@@ -1,5 +1,5 @@
 # coords.minc -- global variables for coordinates
-# $Id$
+# $Header$
 
 # Kronecker delta
 delta,
diff --git a/src/maple/gfa.maple b/src/maple/gfa.maple
index 167b659..b36315c 100644
--- a/src/maple/gfa.maple
+++ b/src/maple/gfa.maple
@@ -1,5 +1,5 @@
 # gfa.maple -- low-level gridfn array "database"
-# $Id$
+# $Header$
 #
 # ***** gridfn arrays *****
 # ***** functional dependences *****
diff --git a/src/maple/gfa.minc b/src/maple/gfa.minc
index daf2cbd..a3adb02 100644
--- a/src/maple/gfa.minc
+++ b/src/maple/gfa.minc
@@ -1,5 +1,5 @@
 # gfa.minc -- global variables for gridfn array system and other misc stuff
-# $Id$
+# $Header$
 
 # arguments to make_gfa
 inert, none,
diff --git a/src/maple/setup.maple b/src/maple/setup.maple
index 6835cb6..931d433 100644
--- a/src/maple/setup.maple
+++ b/src/maple/setup.maple
@@ -1,5 +1,5 @@
 # setup.maple -- top-level Maple input file to set up PDE compiler system
-# $Id$
+# $Header$
 
 #
 # Note not all files here actually *use* the preprocessor, but for
diff --git a/src/maple/util.maple b/src/maple/util.maple
index a099d7a..1ed171d 100644
--- a/src/maple/util.maple
+++ b/src/maple/util.maple
@@ -1,5 +1,5 @@
 # util.maple -- Maple utility functions
-# $Id$
+# $Header$
 
 #
 # msum - multiple-index sum()