finish support for Jacobian queries (not properly tested yet)

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@45 df1f8a13-aa1d-4dd4-9681-27ded5b42416

3 files changed, 415 insertions, 45 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
index eba21b6..bb3bc2b 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
@@ -40,6 +40,8 @@
 static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c)
 
+/* prototypes for static functions (= local to this file) */
+
 /******************************************************************************/
 /******************************************************************************/
 /******************************************************************************/
@@ -190,13 +192,15 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
   @endvar
 
   @var		input_array_offsets
-  @vdesc	(pointer to) array giving an "offset" for each input array,
-		to use in the subscripting computation:  for input array
-		number a, this computation is
-		  data_pointer[offset + i*istride + j*jstride + k*kstride + ...]
+  @vdesc	If this key is present, the value should be a pointer to
+		an array giving an "offset" for each input array, use in
+		the subscripting computation:  for input array number in,
+		this computation (given for 3D; the generalization to other
+		numbers of dimensions is obvious) is
+		  input_arrays[in][offset + i*istride + j*jstride + k*kstride]
 		where
-		  data_pointer = input_arrays[a]
-		  offset = input_array_offsets[a]
+		  offset = input_array_offsets[in]
+			   or is 0 if input_array_offsets is absent
 		  (istride,jstride,kstride,...) = input_array_stride[]
 		and where (i,j,k,...) run from input_array_min_subscripts[]
 		to input_array_max_subscripts[] inclusive.
@@ -339,6 +343,58 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
   @vio		out
   @endvar
 
+  @var		Jacobian_pointer
+  @vdesc	If this key is present, then the value should be an
+		array of N_output_arrays pointers to (caller-supplied)
+		arrays of CCTK_INTs; the interpolator will store the
+		Jacobian elements in the pointed-to arrays.  For output
+		array number out, the subscripting computation (given
+		for 3D; the generalization to other numbers of dimensions
+		is obvious) is
+		 Jacobian_pointer[out][offset
+				       + pt*Jacobian_interp_point_stride
+				       + mi*stride_i + mj*stride_j + mk*stride_k
+				       + part*Jacobian_part_stride]
+		where
+		  offset = Jacobian_offset[out]
+			   or is 0 if input_array_offsets is absent
+		  pt = the point number
+		  (stride_i,stride_j,stride_k) = Jacobian_m_strides[]
+		  part = 0 for real values,
+			 0 for the real parts of complex values,
+			 1 for the imaginary parts of complex values, or
+			 0 if Jacobian_part_stride is absent
+  @vtype	CCTK_REAL *const Jacobian_pointer[N_output_arrays]
+  @endvar
+
+  @var		Jacobian_offset
+  @vdesc	If this key is present, then the value should be a
+		pointer to an array of N_output_arrays CCTK_INTs giving
+		offsets to use in the Jacobian subscripting computation.
+		If this key is absent, the effect is as if a default
+		array of all 0s had been supplied.
+  @vtype	const CCTK_INT Jacobian_offset[N_output_arrays]
+  @endvar
+
+  @var		Jacobian_interp_point_stride
+  @vdesc	The pt-stride for the Jacobian subscripting computation.
+  @vtype	const CCTK_INT Jacobian_interp_point_stride
+  @endvar
+
+  @var		Jacobian_m_strides
+  @vdesc	An array of N_dims CCTK_INTs giving the m strides along
+		each axis for the Jacobian subscripting computation
+  @vtype	const CCTK_INT Jacobian_m_strides[N_dims]
+  @endvar
+
+  @var		Jacobian_part_stride
+  @vdesc	If this key is present, it gives the part stride for the
+		Jacobian subscripting computation.  If this key is absent,
+		the default value 0 is used (n.b. this is suitable only
+		for real data arrays)
+  @vtype	const CCTK_INT Jacobian_m_strides[N_dims]
+  @endvar
+
   ***** return result *****
 
   @returntype	int
@@ -851,7 +907,7 @@ else	{
 /******************************************************************************/
 
 /*
- * set up for molecule min/max m and position queries
+ * set up for any molecule min/max m and position queries
  */
 
 /* molecule min/max m */
@@ -921,6 +977,167 @@ else	{
 /******************************************************************************/
 
 /*
+ * set up for any Jacobian queries
+ */
+
+  {
+struct Jacobian_info *p_Jacobian_info = NULL;
+struct Jacobian_info Jacobian_info;
+Jacobian_info.Jacobian_pointer = NULL;
+Jacobian_info.Jacobian_offset = NULL;
+
+/* first find out if we are doing Jacobian queries at all */
+status = Util_TableQueryValueInfo(param_table_handle,
+				  NULL, NULL,
+				  "Jacobian_pointer");
+if (status < 0)
+   then {
+	free(operation_codes);
+	free(operand_indices);
+	free(input_array_offsets);
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n"
+"                         error status=%d"
+,
+		   status);
+	 return status;					/*** ERROR RETURN ***/
+	}
+
+if (status)
+   then {
+	/* yes, we're doing Jacobian queries */
+	Jacobian_info.Jacobian_pointer
+	    = (CCTK_REAL **) malloc(N_output_arrays * sizeof(CCTK_REAL *));
+	if (Jacobian_info.Jacobian_pointer == NULL)
+	   then {
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
+		}
+	  {
+	CCTK_POINTER *Jacobian_pointer_temp
+	    = (CCTK_POINTER *) malloc(N_output_arrays * sizeof(CCTK_POINTER));
+	if (Jacobian_pointer_temp == NULL)
+	   then {
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		free(Jacobian_info.Jacobian_pointer);
+		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
+		}
+	status = Util_TableGetPointerArray(param_table_handle,
+					   N_output_arrays,
+					   Jacobian_pointer_temp,
+					   "Jacobian_pointer");
+	if	(status == N_output_arrays)
+	   then {
+		/* type-cast CCTK_POINTER pointers to CCTK_REAL* pointers */
+		/* (which point to where the query results will be stored) */
+		int out;
+			for (out = 0 ; out < N_output_arrays ; ++out)
+			{
+			Jacobian_info.Jacobian_pointer[out]
+				= (CCTK_REAL *) Jacobian_pointer_temp[out];
+			}
+		free(Jacobian_pointer_temp);
+		}
+	else	{
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		free(Jacobian_info.Jacobian_pointer);
+		free(Jacobian_pointer_temp);
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n"
+"                         error status=%d"
+,
+			   status);
+		return status;				/*** ERROR RETURN ***/
+		}
+	  }
+
+	Jacobian_info.Jacobian_offset
+		= (CCTK_INT *) malloc(N_output_arrays * sizeof(CCTK_INT));
+	if (Jacobian_info.Jacobian_offset == NULL)
+	   then {
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		free(Jacobian_info.Jacobian_pointer);
+		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
+		}
+	status = Util_TableGetIntArray(param_table_handle,
+				       N_output_arrays,
+				       Jacobian_info.Jacobian_offset,
+				       "Jacobian_offset");
+	if	(status == N_output_arrays)
+	   then {
+		/* ok ==> no-op here */
+		}
+	else if (status2 == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+	   then {
+		/* default offset = all 0 */
+		LocalInterp_zero_int_array(N_output_arrays,
+					   Jacobian_info.Jacobian_offset);
+		}
+	else	{
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		free(Jacobian_info.Jacobian_pointer);
+		free(Jacobian_info.Jacobian_offset);
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   InterpLocalUniform(): bad \"Jacobian_offset\" table entry!\n"
+"                         error status=%d"
+,
+			   status);
+		return status;				/*** ERROR RETURN ***/
+		}
+
+	status1 = Util_TableGetInt(param_table_handle,
+				   & Jacobian_info.Jacobian_interp_point_stride,
+				   "Jacobian_interp_point_stride");
+	status2 = Util_TableGetIntArray(param_table_handle,
+					N_dims,
+					Jacobian_info.Jacobian_m_strides,
+					"Jacobian_m_strides");
+	status3 = Util_TableGetInt(param_table_handle,
+				   & Jacobian_info.Jacobian_part_stride,
+				   "Jacobian_part_stride");
+	if ((status1 < 0) || (status2 != N_dims) || (status3 < 0))
+	   then {
+		free(operation_codes);
+		free(operand_indices);
+		free(input_array_offsets);
+		free(Jacobian_info.Jacobian_pointer);
+		free(Jacobian_info.Jacobian_offset);
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   InterpLocalUniform(): bad \"Jacobian_interp_point_stride and/or\n"
+"                         \"Jacobian_m_strides\" and/or\n"
+"                         \"Jacobian_part_stride\" table entry/entries!\n"
+"                         error status1=%d status2=%d status3=%d"
+	,
+			   status1, status2, status3);
+		return   (status1 < 0) ? status1	/*** ERROR RETURN ***/
+		       : (status2 < 0) ? status2	/*** ERROR RETURN ***/
+		       :		 status3;	/*** ERROR RETURN ***/
+		}
+
+	p_Jacobian_info = & Jacobian_info;
+	}
+
+/******************************************************************************/
+
+/*
  * decode (N_dims, molecule_family, order, smoothing)
  * and call the appropriate subfunction to do the actual interpolation
  */
@@ -1017,6 +1234,11 @@ if (p_interp_fn == NULL_INTERP_FN_PTR)
 	free(operation_codes);
 	free(operand_indices);
 	free(input_array_offsets);
+	if (p_Jacobian_info != NULL)
+	   then {
+		free(Jacobian_info.Jacobian_pointer);
+		free(Jacobian_info.Jacobian_offset);
+		}
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -1050,7 +1272,7 @@ const int return_code
 			 &molecule_structure_flags,
 			 p_molecule_min_max_m_info,
 			 p_molecule_positions,
-			 NULL);
+			 p_Jacobian_info);
 
 /******************************************************************************/
 
@@ -1075,6 +1297,11 @@ if (return_code == CCTK_ERROR_INTERP_POINT_X_RANGE)
 		free(operation_codes);
 		free(operand_indices);
 		free(input_array_offsets);
+		if (p_Jacobian_info != NULL)
+		   then {
+			free(Jacobian_info.Jacobian_pointer);
+			free(Jacobian_info.Jacobian_offset);
+			}
 		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 			   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -1112,6 +1339,11 @@ if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0))
 	free(operation_codes);
 	free(operand_indices);
 	free(input_array_offsets);
+	if (p_Jacobian_info != NULL)
+	   then {
+		free(Jacobian_info.Jacobian_pointer);
+		free(Jacobian_info.Jacobian_offset);
+		}
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -1149,6 +1381,11 @@ if (p_molecule_min_max_m_info != NULL)
 		free(operation_codes);
 		free(operand_indices);
 		free(input_array_offsets);
+		if (p_Jacobian_info != NULL)
+		   then {
+			free(Jacobian_info.Jacobian_pointer);
+			free(Jacobian_info.Jacobian_offset);
+			}
 		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 			   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -1165,6 +1402,11 @@ if (p_molecule_min_max_m_info != NULL)
 free(operation_codes);
 free(operand_indices);
 free(input_array_offsets);
+if (p_Jacobian_info != NULL)
+   then {
+	free(Jacobian_info.Jacobian_pointer);
+	free(Jacobian_info.Jacobian_offset);
+	}
 
 return return_code;
   }
@@ -1181,4 +1423,5 @@ return return_code;
   }
   }
   }
+  }
 }
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
index 3939c4a..32ee294 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
@@ -125,10 +125,10 @@ struct	error_flags
 
 struct	molecule_structure_flags
 	{
-	int* MSS_is_fn_of_interp_coords;
-	int* MSS_is_fn_of_which_operation;
-	int* MSS_is_fn_of_input_array_values;
-	int* Jacobian_is_fn_of_input_array_values;
+	bool MSS_is_fn_of_interp_coords;
+	bool MSS_is_fn_of_which_operation;
+	bool MSS_is_fn_of_input_array_values;
+	bool Jacobian_is_fn_of_input_array_values;
 	};
 
 struct	molecule_min_max_m_info
@@ -139,8 +139,8 @@ struct	molecule_min_max_m_info
 
 struct	Jacobian_info
 	{
-	CCTK_REAL* const* Jacobian_pointer;
-	const CCTK_INT*   Jacobian_offset;
+	CCTK_REAL** Jacobian_pointer;
+	CCTK_INT*   Jacobian_offset;
 	CCTK_INT Jacobian_interp_point_stride;
 	CCTK_INT Jacobian_m_strides[MAX_N_DIMS];
 	CCTK_INT Jacobian_part_stride;
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index 900f11b..32d0203 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -390,6 +390,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
  * Naming conventions:
  * in, out = 0-origin indices each selecting an input/output array
  * pt = 0-origin index selecting an interpolation point
+ * part = 0-origin index selecting real/imaginary part of a complex number
  */
 
 /* number of real "parts" in a complex number */
@@ -936,11 +937,11 @@ case CCTK_VARIABLE_REAL:
 		= (const CCTK_REAL *) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_real[sub_temp + SUB1(i)]
+    #define DATA(mi)	    input_array_ptr_real[sub_temp + SUB1(mi)]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_real[sub_temp + SUB2(i,j)]
+    #define DATA(mi,mj)	    input_array_ptr_real[sub_temp + SUB2(mi,mj)]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_real[sub_temp + SUB3(i,j,k)]
+    #define DATA(mi,mj,mk)  input_array_ptr_real[sub_temp + SUB3(mi,mj,mk)]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -955,11 +956,11 @@ case CCTK_VARIABLE_REAL4:
 		= (const CCTK_REAL4 *) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_real4[sub_temp + SUB1(i)]
+    #define DATA(mi)	    input_array_ptr_real4[sub_temp + SUB1(mi)]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_real4[sub_temp + SUB2(i,j)]
+    #define DATA(mi,mj)	    input_array_ptr_real4[sub_temp + SUB2(mi,mj)]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_real4[sub_temp + SUB3(i,j,k)]
+    #define DATA(mi,mj,mk)  input_array_ptr_real4[sub_temp + SUB3(mi,mj,mk)]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -975,11 +976,11 @@ case CCTK_VARIABLE_REAL8:
 		= (const CCTK_REAL8 *) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_real8[sub_temp + SUB1(i)]
+    #define DATA(mi)	    input_array_ptr_real8[sub_temp + SUB1(mi)]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_real8[sub_temp + SUB2(i,j)]
+    #define DATA(mi,mj)	    input_array_ptr_real8[sub_temp + SUB2(mi,mj)]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_real8[sub_temp + SUB3(i,j,k)]
+    #define DATA(mi,mj,mk)  input_array_ptr_real8[sub_temp + SUB3(mi,mj,mk)]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -997,11 +998,11 @@ case CCTK_VARIABLE_REAL16:
 		= (const CCTK_REAL16 *) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_real16[sub_temp + SUB1(i)]
+    #define DATA(mi)	    input_array_ptr_real16[sub_temp + SUB1(mi)]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_real16[sub_temp + SUB2(i,j)]
+    #define DATA(mi,mj)	    input_array_ptr_real16[sub_temp + SUB2(mi,mj)]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_real16[sub_temp + SUB3(i,j,k)]
+    #define DATA(mi,mj,mk)  input_array_ptr_real16[sub_temp + SUB3(mi,mj,mk)]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -1016,11 +1017,14 @@ case CCTK_VARIABLE_COMPLEX:
 		= (const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_complex[sub_temp + SUB1(i)][part]
+    #define DATA(mi)	    \
+	    input_array_ptr_complex[sub_temp + SUB1(mi)][part]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_complex[sub_temp + SUB2(i,j)][part]
+    #define DATA(mi,mj)	    \
+	    input_array_ptr_complex[sub_temp + SUB2(mi,mj)][part]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_complex[sub_temp + SUB3(i,j,k)][part]
+    #define DATA(mi,mj,mk)  \
+	    input_array_ptr_complex[sub_temp + SUB3(mi,mj,mk)][part]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -1035,11 +1039,14 @@ case CCTK_VARIABLE_COMPLEX8:
 		= (const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_complex8[sub_temp + SUB1(i)][part]
+    #define DATA(mi)	    \
+	    input_array_ptr_complex8[sub_temp + SUB1(mi)][part]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_complex8[sub_temp + SUB2(i,j)][part]
+    #define DATA(mi,mj)	    \
+	    input_array_ptr_complex8[sub_temp + SUB2(mi,mj)][part]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_complex8[sub_temp + SUB3(i,j,k)][part]
+    #define DATA(mi,mj,mk)  \
+	    input_array_ptr_complex8[sub_temp + SUB3(mi,mj,mk)][part]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -1055,11 +1062,14 @@ case CCTK_VARIABLE_COMPLEX16:
 		= (const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_complex16[sub_temp + SUB1(i)][part]
+    #define DATA(mi)	    \
+	    input_array_ptr_complex16[sub_temp + SUB1(mi)][part]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_complex16[sub_temp + SUB2(i,j)][part]
+    #define DATA(mi,mj)	    \
+	    input_array_ptr_complex16[sub_temp + SUB2(mi,mj)][part]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_complex16[sub_temp + SUB3(i,j,k)][part]
+    #define DATA(mi,mj,mk)  \
+	    input_array_ptr_complex16[sub_temp + SUB3(mi,mj,mk)][part]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -1075,11 +1085,14 @@ case CCTK_VARIABLE_COMPLEX32:
 		= (const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr;
 	#undef DATA
   #if   (N_DIMS == 1)
-    #define DATA(i)	input_array_ptr_complex32[sub_temp + SUB1(i)][part]
+    #define DATA(mi)	    \
+	    input_array_ptr_complex32[sub_temp + SUB1(mi)][part]
   #elif (N_DIMS == 2)
-    #define DATA(i,j)	input_array_ptr_complex32[sub_temp + SUB2(i,j)][part]
+    #define DATA(mi,mj)	    \
+	    input_array_ptr_complex32[sub_temp + SUB2(mi,mj)][part]
   #elif (N_DIMS == 3)
-    #define DATA(i,j,k)	input_array_ptr_complex32[sub_temp + SUB3(i,j,k)][part]
+    #define DATA(mi,mj,mk)  \
+	    input_array_ptr_complex32[sub_temp + SUB3(mi,mj,mk)][part]
   #else
     #error "N_DIMS must be 1, 2, or 3!"
   #endif
@@ -1176,11 +1189,11 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 				}
 
 
-			/*
-			 * ***store*** the result in the output array
-			 */
-			switch	(output_array_type_codes[out])
-				{
+/*
+ * ***store*** the result in the output array
+ */
+switch	(output_array_type_codes[out])
+	{
 
 case CCTK_VARIABLE_REAL:
 	  {
@@ -1263,8 +1276,122 @@ default:
 		   "output datatype %d not supported",
 		   output_array_type_codes[out]);		/*NOTREACHED*/
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
-				/* end of switch (output type code) */
-				}
+	/* end of switch (output type code) */
+	}
+
+/*
+ * handle querying the Jacobian
+ */
+if ((Jacobian_info != NULL) && (Jacobian_info->Jacobian_pointer[out] != NULL))
+   then {
+	/*
+	 * set up the Jacobian storage addressing
+	 */
+	const int J_offset      = Jacobian_info->Jacobian_offset[out];
+	const int J_pt_stride   = Jacobian_info->Jacobian_interp_point_stride;
+	const int J_part_stride = Jacobian_info->Jacobian_part_stride;
+	double *const Jacobian_ptr
+		= & Jacobian_info->Jacobian_pointer[out][J_offset
+							 + pt*J_pt_stride
+							 + part*J_part_stride];
+
+	#if (N_DIMS >= 1)
+	  const int J_mx_stride = Jacobian_info->Jacobian_m_strides[X_AXIS];
+	#endif
+	#if (N_DIMS >= 2)
+	  const int J_my_stride = Jacobian_info->Jacobian_m_strides[Y_AXIS];
+	#endif
+	#if (N_DIMS >= 3)
+	  const int J_mz_stride = Jacobian_info->Jacobian_m_strides[Z_AXIS];
+	#endif
+
+	#if   (N_DIMS == 1)
+	  #define COEFF(mi)	Jacobian_ptr[mi*J_mx_stride]
+	#elif (N_DIMS == 2)
+	  #define COEFF(mi,mj)	Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride]
+	#elif (N_DIMS == 3)
+	  #define COEFF(mi,mj,mk)	\
+		Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride + mk*J_mz_stride]
+	#else
+	  #error "N_DIMS must be 1, 2, or 3!"
+	#endif
+
+	  {
+	fp factor;
+	switch	(operation_codes[out])
+		{
+	#ifdef HAVE_OP_I
+	  case 0:
+		#include COEFF_I_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DX
+	  case 1:
+		factor = inverse_delta_x;
+		#include COEFF_DX_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DY
+	  case 2:
+		factor = inverse_delta_y;
+		#include COEFF_DY_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DZ
+	  case 3:
+		factor = inverse_delta_z;
+		#include COEFF_DZ_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DXX
+	  case 11:
+		factor = inverse_delta_x2;
+		#include COEFF_DXX_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DXY
+	  case 12:
+	  case 21:
+		factor = inverse_delta_x * inverse_delta_y;
+		#include COEFF_DXY_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DXZ
+	  case 13:
+	  case 31:
+		factor = inverse_delta_x * inverse_delta_z;
+		#include COEFF_DXZ_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DYY
+	  case 22:
+		factor = inverse_delta_y2;
+		#include COEFF_DYY_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DYZ
+	  case 23:
+	  case 32:
+		factor = inverse_delta_y * inverse_delta_z;
+		#include COEFF_DYZ_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DZZ
+	  case 33:
+		factor = inverse_delta_z2;
+		#include COEFF_DZZ_VAR_STORE_FILE_NAME
+		break;
+	#endif
+	  default:
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	"Generalized interpolation operation_code %d not supported!",
+			   operation_codes[out]);		/*NOTREACHED*/
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		/* end of switch(operation_codes[out])*/
+		}
+	  }
+	/* end of Jacobian-query code */
+	}
 			  }
 			/* end of  for (part = ...)  loop */
 			}