initial checkin of new LocalInterp thorn

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

1 files changed, 225 insertions, 0 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/molecule_posn.c b/src/GeneralizedPolynomial-Uniform/molecule_posn.c
new file mode 100644
index 0000000..182a3e4
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/molecule_posn.c
@@ -0,0 +1,225 @@
+ /*@@
+   @file      molecule_posn.c
+   @date      23 Oct 2001
+   @author    Jonathan Thornburg <jthorn@aei.mpg.de>
+   @desc      Worker function to compute molecule positions for interpolator.
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <limits.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#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 <jthorn@aei.mpg.de>
+ @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	out_of_range_tolerance
+ @vdesc	Specifies how out-of-range interpolation points should
+	be handled:
+	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
+ @endvar
+
+ @var	molecule_size
+ @vdesc	The size (number of points) of the molecule.
+ @vtype	int molecule_size
+ @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
+ @vdesc	A pointer to an int where this function should store the minimum
+	molecule coordinate m of the molecule; or NULL to skip storing this.
+ @vtype	int *molecule_m_min
+ @vio	pointer to out
+ @endvar
+
+ @var	p_molecule_m_max
+ @vdesc	A pointer to an int where this function should store the maximum
+	molecule coordinate m of the molecule; or NULL to skip storing this.
+ @vtype	int *molecle_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 <limits.h>)
+	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,
+			      fp out_of_range_tolerance,
+			      int molecule_size,
+			      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 >= 0.0)
+   then {
+	const fp fp_effective_grid_i_min
+		= ((fp) grid_i_min) - out_of_range_tolerance;
+	const fp fp_effective_grid_i_max
+		= ((fp) grid_i_max) + out_of_range_tolerance;
+	if (fp_i < fp_effective_grid_i_min)
+	   then return INT_MIN;				/*** ERROR RETURN ***/
+	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 == -1.0)
+   then {
+	if (centered_i_min < grid_i_min)
+	   then return INT_MIN;				/*** ERROR RETURN ***/
+	if (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;
+  }
+  }
+}