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

1 files changed, 158 insertions, 31 deletions

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 */
 			}