* finish Jacobian support in src/util

* test driver for this in src/util/test_patch_system.cc
  --> IT WORKS!!!!!
* also fix comments on #include prerequisites to say ../jtutil/ not jt/


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

1 files changed, 109 insertions, 120 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 0c7181f..c19dd42 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -4,14 +4,13 @@
 //
 // patch_interp::patch_interp
 // patch_interp::~patch_interp
-// patch_interp::compute_[min,max]_ipar
 // patch_interp::interpolate
 //
-// patch_interp::query_interpolator
 // patch_interp::verify_Jacobian_sparsity_pattern_ok
 // patch_interp::molecule_minmax_ipar_m
 // patch_interp::molecule_posn
 // patch_interp::Jacobian
+// patch_interp::query_interpolator
 //
 
 #include <stdio.h>
@@ -56,22 +55,28 @@ 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,
 			   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()),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
-	  min_ipar_(compute_min_ipar()),
-	  max_ipar_(compute_max_ipar()),
+	  min_ipar_(ok_to_use_min_par_corner
+		    ? my_edge_in.min_ipar_with_corners()
+		    : my_edge_in.min_ipar_without_corners()),
+	  max_ipar_(ok_to_use_max_par_corner
+		    ? my_edge_in.max_ipar_with_corners()
+		    : my_edge_in.max_ipar_without_corners()),
 	  min_parindex_array_(min_parindex_array_in),
 	  max_parindex_array_(max_parindex_array_in),
 	  interp_par_(interp_par_in),
 	  interp_handle_(interp_handle_in),
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
-	  // note interpolator coordinate origin = par of min_ipar() point
+	  // note interpolator coordinate origin = min_ipar_
 	  // cf. comments in "patch_interp.hh" on gridfn array subscripting
-	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())),
+	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar_)),
 	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
 	  gridfn_type_codes_      (min_gfn_, max_gfn_),
 	  gridfn_data_ptrs_       (min_gfn_, max_gfn_),
@@ -132,35 +137,6 @@ Util_TableDestroy(interp_par_table_handle_);
 //*****************************************************************************
 
 //
-// These functions compute the [min,max] range of ipar in this patch
-// from which we may use gridfn data in interpolation.  This computation
-// depends on our adjacent ghost zones' types as described in the comments
-// to our constructor, but needs only my_patch_ and my_edge_ members of
-// this object to be valid.
-//
-int patch_interp::compute_min_ipar() const
-{
-const ghost_zone& min_par_adjacent_ghost_zone
-	= my_patch()
-	  .ghost_zone_on_edge(my_edge().min_par_adjacent_edge());
-return min_par_adjacent_ghost_zone.is_symmetry()
-       ? my_edge().min_ipar_with_corners()
-       : my_edge().min_ipar_without_corners();
-}
-
-int patch_interp::compute_max_ipar() const
-{
-const ghost_zone& max_par_adjacent_ghost_zone
-	= my_patch()
-	  .ghost_zone_on_edge(my_edge().max_par_adjacent_edge());
-return max_par_adjacent_ghost_zone.is_symmetry()
-       ? my_edge().max_ipar_with_corners()
-       : my_edge().max_ipar_without_corners();
-}
-
-//*****************************************************************************
-
-//
 // This function interpolates the specified range of ghosted gridfns
 // at all the coordinates specified when we were constructed.  (It does
 // a separate interpolator call for each iperp.)  This function stores
@@ -259,89 +235,6 @@ const CCTK_INT *gridfn_type_codes_ptr
 //******************************************************************************
 
 //
-// This function queries the interpolator with whatever arguments are
-// in the parameter table.  It specifies arguments such that no actual
-// interpolation is done.
-//
-// In particular, the following interpolator arguments are set up properly
-// based on  iperp :
-//	N_dims
-//	operator_handle, param_table_handle,
-//	coord_origin, coord_delta,
-//	N_interp_points, interp_coords_type_code, interp_coords,
-//	input_array_dims	# specifies the correct par size of the
-//				# patch interpolation region at this iperp
-// The following arguments are set to specify a single input array, but
-// with a NULL data pointer so no actual data is used:
-//	N_input_arrays
-//	input_array_type_codes, input_arrays
-// The following arguments are set to specify a single output array, but
-// with a NULL data pointer so no actual data is stored:
-//	N_output_arrays
-//	output_array_type_codes, output_arrays
-//
-// Arguments:
-// iperp = (in) Specifies where in the patch interpolation region the
-//		interpolator query should be done.
-//
-// Results:
-// If the interpolator returns an error status, this function does an
-//  error_exit()  (and does not return to its caller).  Otherwise, this
-// function returns the interpolator status code.
-//
-int patch_interp::query_interpolator(const char function_name[], int iperp)
-	const
-{
-int status;
-const int N_dims = 1;
-
-const int min_parindex = min_parindex_array_(iperp);
-const int max_parindex = max_parindex_array_(iperp);
-const int N_interp_points = jtutil::how_many_in_range(min_parindex,
-						      max_parindex);
-const int interp_coords_type_code = CCTK_VARIABLE_REAL;
-const void* const interp_coords[N_dims]
-	= { static_cast<const void *>(& interp_par_(iperp, min_parindex)) };
-
-const int N_input_arrays = 1;
-const CCTK_INT input_array_dims[N_dims]
-	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
-const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
-const void* const input_arrays[N_input_arrays] = { NULL };
-
-const int N_output_arrays = 1;
-const CCTK_INT output_array_type_codes[N_output_arrays]
-	= { CCTK_VARIABLE_REAL };
-void* const output_arrays[N_output_arrays] = { NULL };
-
-status = CCTK_InterpLocalUniform(N_dims,
-				 interp_handle_, interp_par_table_handle_,
-				 &gridfn_coord_origin_, &gridfn_coord_delta_,
-				 N_interp_points,
-				    interp_coords_type_code,
-				    interp_coords,
-				 N_input_arrays,
-				    input_array_dims,
-				    input_array_type_codes,
-				    input_arrays,
-				 N_output_arrays,
-				    output_array_type_codes,
-				    output_arrays);
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::query_interpolator():\n"
-"        on behalf of patch_interp::%s()\n"
-"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
-,
-		   function_name,
-		   status, iperp, min_iperp(), max_iperp());	/*NOTREACHED*/
-
-return status;
-}
-
-//******************************************************************************
-
-//
 // This function verifies that our interpolator has a Jacobian sparsity
 // pattern which we grok: at present this means that the molecules are
 // fixed-sized hypercubes, with size/shape independent of the interpolation
@@ -505,6 +398,7 @@ int status;
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
 	const int min_parindex = min_parindex_array_(iperp);
+	const int max_parindex = max_parindex_array_(iperp);
 
 	// set up the molecule-position query in the parameter table
 	CCTK_POINTER molecule_posn_ptrs[N_dims]
@@ -523,9 +417,21 @@ int status;
 			   status);				/*NOTREACHED*/
 
 	// query the interpolator to get the molecule positions
-	// for all parindex at this iperp (the interpolator will
-	// store them directly in the  posn_buffer(iperp,*)  array)
+	// for all parindex at this iperp; the interpolator will
+	// store the  parindex-min_parindex  values (cf comments
+	// on array subscripting at the start of "patch_interp.hh")
+	// directly in the  posn_buffer(iperp,*)  array)
 	query_interpolator("molecule_positions", iperp);
+
+	// convert the molecule positions from  parindex-min_ipar
+	// to  parindex  values (again, cf comments on array subscripting
+	// at the start of "patch_interp.hh")
+		for (int parindex = min_parindex ;
+		     parindex <= max_parindex ;
+		     ++parindex)
+		{
+		posn_buffer(iperp, parindex) += min_ipar();
+		}
 	}
 
 //
@@ -658,3 +564,86 @@ if ((status < 0) || (status1 < 0) || (status2 < 0))
 ,
 		   status, status1, status2);			/*NOTREACHED*/
 }
+
+//******************************************************************************
+
+//
+// This function queries the interpolator with whatever arguments are
+// in the parameter table.  It specifies arguments such that no actual
+// interpolation is done.
+//
+// In particular, the following interpolator arguments are set up properly
+// based on  iperp :
+//	N_dims
+//	operator_handle, param_table_handle,
+//	coord_origin, coord_delta,
+//	N_interp_points, interp_coords_type_code, interp_coords,
+//	input_array_dims	# specifies the correct par size of the
+//				# patch interpolation region at this iperp
+// The following arguments are set to specify a single input array, but
+// with a NULL data pointer so no actual data is used:
+//	N_input_arrays
+//	input_array_type_codes, input_arrays
+// The following arguments are set to specify a single output array, but
+// with a NULL data pointer so no actual data is stored:
+//	N_output_arrays
+//	output_array_type_codes, output_arrays
+//
+// Arguments:
+// iperp = (in) Specifies where in the patch interpolation region the
+//		interpolator query should be done.
+//
+// Results:
+// If the interpolator returns an error status, this function does an
+//  error_exit()  (and does not return to its caller).  Otherwise, this
+// function returns the interpolator status code.
+//
+int patch_interp::query_interpolator(const char function_name[], int iperp)
+	const
+{
+int status;
+const int N_dims = 1;
+
+const int min_parindex = min_parindex_array_(iperp);
+const int max_parindex = max_parindex_array_(iperp);
+const int N_interp_points = jtutil::how_many_in_range(min_parindex,
+						      max_parindex);
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const void* const interp_coords[N_dims]
+	= { static_cast<const void *>(& interp_par_(iperp, min_parindex)) };
+
+const int N_input_arrays = 1;
+const CCTK_INT input_array_dims[N_dims]
+	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
+const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
+const void* const input_arrays[N_input_arrays] = { NULL };
+
+const int N_output_arrays = 1;
+const CCTK_INT output_array_type_codes[N_output_arrays]
+	= { CCTK_VARIABLE_REAL };
+void* const output_arrays[N_output_arrays] = { NULL };
+
+status = CCTK_InterpLocalUniform(N_dims,
+				 interp_handle_, interp_par_table_handle_,
+				 &gridfn_coord_origin_, &gridfn_coord_delta_,
+				 N_interp_points,
+				    interp_coords_type_code,
+				    interp_coords,
+				 N_input_arrays,
+				    input_array_dims,
+				    input_array_type_codes,
+				    input_arrays,
+				 N_output_arrays,
+				    output_array_type_codes,
+				    output_arrays);
+if (status < 0)
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::query_interpolator():\n"
+"        on behalf of patch_interp::%s()\n"
+"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
+,
+		   function_name,
+		   status, iperp, min_iperp(), max_iperp());	/*NOTREACHED*/
+
+return status;
+}