diff options
| -rw-r--r-- | src/InterpGridArrays.c | 29 |
1 files changed, 24 insertions, 5 deletions
diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c index 39de534..3c95211 100644 --- a/src/InterpGridArrays.c +++ b/src/InterpGridArrays.c @@ -203,7 +203,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH, CCTK_REAL *range_min, *range_max; CCTK_REAL *origin_global; CCTK_REAL *interp_coords_proc, *coords, *bbox_local; - const CCTK_REAL **bbox_interp_coords; + CCTK_REAL **bbox_interp_coords; #endif @@ -566,17 +566,29 @@ int PUGHInterp_InterpGridArrays (const cGH *GH, This is done as a query, simply by passing all arguments as for the real call to the local interpolator, except for the output_arrays[] argument which contains NULL pointers. */ - bbox_interp_coords = malloc (N_dims * sizeof (void *)); - for (i = 0; i < N_dims; i++) + bbox_interp_coords = malloc (N_dims * sizeof (CCTK_REAL *)); + bbox_interp_coords[0] = malloc (2 * N_dims * N_dims * sizeof (CCTK_REAL)); + for (i = 1; i < N_dims; i++) + { + bbox_interp_coords[i] = bbox_interp_coords[i-1] + 2*N_dims; + } + for (j = 0; j < 2 * N_dims; j++) { - bbox_interp_coords[i] = &bbox_local[2 * i]; + for (i = 0; i < N_dims; i++) + { + bbox_interp_coords[i][j] = 0.5 * (bbox_local[2*i+0] + bbox_local[2*i+1]); + } + if (! GH->cctk_bbox[j]) + { + bbox_interp_coords[j/2][j] = bbox_local[j]; + } } output_arrays_local = N_output_arrays > 0 ? calloc (N_output_arrays, sizeof (void *)) : NULL; retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle, param_table_handle, - origin_local, delta, 2, + origin_local, delta, 2 * N_dims, CCTK_VARIABLE_REAL, (const void *const *) bbox_interp_coords, N_input_arrays, input_array_dims, @@ -612,6 +624,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH, { free (output_arrays_local); } + free (bbox_interp_coords[0]); free (bbox_interp_coords); free (range_min); @@ -627,6 +640,12 @@ int PUGHInterp_InterpGridArrays (const cGH *GH, free (input_array_dims); free (type_desc); + /* if the query call returned 'point-outside' error + we turn this into a 'ghost-size-too-small' error here */ + if (retval == CCTK_ERROR_INTERP_POINT_OUTSIDE) + { + retval = CCTK_ERROR_INTERP_GHOST_SIZE_TOO_SMALL; + } return (retval); } |