/*@@ @file molecule_posn.c @date 23 Oct 2001 @author Jonathan Thornburg @desc Worker function to compute molecule positions for interpolator. @enddesc @version $Id$ @@*/ #include #include #include #include #ifndef LOCALINTERP_STANDALONE_BUILD #include "cctk.h" #endif #include "InterpLocalUniform.h" /* the rcs ID and its dummy function to use it */ static const char *rcsid = "$Header$"; #ifndef LOCALINTERP_STANDALONE_BUILD CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_molecule_posn_c) #endif /******************************************************************************/ /******************************************************************************/ /******************************************************************************/ /*@@ @routine LocalInterp_molecule_posn @date 23 Oct 2001 @author Jonathan Thornburg @desc Given a uniformly-spaced grid, this function computes molecule positions for an interpolation or similar operation: (1) It computes the closest grid point to the input coordinate (2) It does the slightly tricky odd/even computation to decide where to center the molecule (3) It checks for the molecule falling off the edge of the grid, and off-centers it as appropriate. For (1), we assume that grid points have floating-point coordinates (we refer to these as "x") which are an arbitrary linear function of the integer grid coordinates (we refer to these as "i"), x = grid_origin + i*grid_delta. For (2), suppose we have a data array indexed by [i], from which we wish to select an N-point molecule [i_lo, i_hi] centered as close as possible to the point x . The problem is just how to choose the centering. The following diagram illustrates the issues: N i_lo i_hi [i-3] [i-2] [i-1] [i] [i+1] [i+2] [i+3] [i+4] - ---- ---- ----- ----- ----- ----- ----- ----- ----- ----- 2 i i+1 *xxxxxx* 3 i-1 i+1 * xxx*xxx * 4 i-1 i+2 * *xxxxxx* * 5 i-2 i+2 * * xxx*xxx * * 6 i-2 i+3 * * *xxxxxx* * * 7 i-3 i+3 * * * xxx*xxx * * * 8 i-3 i+4 * * * *xxxxxx* * * * The diagram shows the range of x values relative to the grid, for which we'll choose each centering. Thus if N is even, the centering is based on which grid zone contains x (take the floor of the floating-point i coordinate), while if N is odd, the centering is based on which grid point is closest to x (round the floating- -point i coordinate to the nearest integer). It's also convenient to introduce the integer molecule coordinate m; this is the integer grid coordinate i relative to that of the molecule center, i.e. the above diagram has columns labelled with [i+m]. @enddesc @var grid_origin @vdesc The floating-point coordinate x of the grid point i=0. @vtype fp grid_origin @endvar @var grid_delta @vdesc The grid spacing (in the floating-point coordinate x). @vtype fp grid_delta @endvar @var grid_i_min @vdesc The minimum integer coordinate i of the grid. @vtype int grid_i_min @endvar @var grid_i_max @vdesc The maximum integer coordinate i of the grid. @vtype int grid_i_max @endvar @var molecule_size @vdesc The size (number of points) of the molecule. @vtype int molecule_size @endvar @var out_of_range_tolerance_min, out_of_range_tolerance_max @vdesc Specifies how out-of-range interpolation points should be handled for the {minimum,maximum} ends of the grid respectively. If out_of_range_tolerance >= 0.0, then an interpolation point is considered to be "out of range" if and only if the interpolation point is > out_of_range_tolerance * grid_delta outside the grid in any coordinate. If out_of_range_tolerance == -1.0, then an interpolation point is considered to be "out of range" if and only if a centered molecule would require data from outside the grid. Other values of out_of_range_tolerance are illegal. @vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max @endvar @var x @vdesc The floating-point coordinate x at which we would like to center the molecule. @vtype fp x @endvar @var x_rel @vdesc A pointer to a floating-point value where this function should store the x coordinate relative to the molecule center, measured in units of the grid spacing; or NULL to skip storing this. @vtype fp *x_rel @vio pointer to out @endvar @var molecule_m_min, molecule_m_max @vdesc A pointer to an int where this function should store the {minimum,maximum} molecule coordinate m of the molecule; or NULL to skip storing this. @vtype int *molecule_m_min, *molecule_m_max @vio pointer to out @endvar @returntype int @returndesc This function returns the integer coordinate i_center at which the molecule should be centered, or one of the error codes (from ) INT_MIN if x is out-of-range on the min end of the grid (i.e. x < the minimum allowable x) INT_MAX if x is out-of-range on the max end of the grid (i.e. x > the maximum allowable x) @endreturndesc @@*/ int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, fp out_of_range_tolerance_min, fp out_of_range_tolerance_max, fp x, fp *x_rel, int *molecule_m_min, int *molecule_m_max) { /* min/max molecule coordinates */ const int m_max = (molecule_size >> 1); /* a positive number */ const int m_min = m_max - molecule_size + 1; /* a negative number */ /* integer coordinate i of point x, as a floating-point number */ const fp fp_i = (x - grid_origin) / grid_delta; /* is point x out-of-range? */ if (out_of_range_tolerance_min >= 0.0) then { const fp fp_effective_grid_i_min = ((fp) grid_i_min) - out_of_range_tolerance_min; if (fp_i < fp_effective_grid_i_min) then return INT_MIN; /*** ERROR RETURN ***/ } if (out_of_range_tolerance_max >= 0.0) then { const fp fp_effective_grid_i_max = ((fp) grid_i_max) + out_of_range_tolerance_max; if (fp_i > fp_effective_grid_i_max) then return INT_MAX; /*** ERROR RETURN ***/ } /* integer coordinate i where we will center the molecule */ /* (see diagram in header comment for explanation) */ /* ... as a floating-point number */ { const fp fp_i_center = is_even(molecule_size) ? floor(fp_i) : jt_round(fp_i); /* ... as an integer */ int i_center = (int) fp_i_center; /* min/max integer coordinates in molecule assuming it's fully centered */ const int centered_i_min = i_center + m_min; const int centered_i_max = i_center + m_max; /* check if off-centered molecules are forbidden? */ if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min)) then return INT_MIN; /*** ERROR RETURN ***/ if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max)) then return INT_MAX; /*** ERROR RETURN ***/ /* off-center as needed if we're close to the edge of the grid */ { int shift = 0; if (centered_i_min < grid_i_min) then shift = grid_i_min - centered_i_min; /* a positive number */ if (centered_i_max > grid_i_max) then shift = grid_i_max - centered_i_max; /* a negative number */ i_center += shift; /* store the results */ if (x_rel != NULL) then *x_rel = fp_i - ((fp) i_center); if (molecule_m_min != NULL) then *molecule_m_min = m_min; if (molecule_m_max != NULL) then *molecule_m_max = m_max; return i_center; } } }