diff options
Diffstat (limited to 'src/UniformCartesian')
| -rw-r--r-- | src/UniformCartesian/Interpolate.c | 609 | ||||
| -rw-r--r-- | src/UniformCartesian/Interpolate.h | 45 | ||||
| -rw-r--r-- | src/UniformCartesian/Operator.c | 323 | ||||
| -rw-r--r-- | src/UniformCartesian/README | 14 | ||||
| -rw-r--r-- | src/UniformCartesian/Startup.c | 185 | ||||
| -rw-r--r-- | src/UniformCartesian/make.code.defn | 4 |
6 files changed, 0 insertions, 1180 deletions
diff --git a/src/UniformCartesian/Interpolate.c b/src/UniformCartesian/Interpolate.c deleted file mode 100644 index 7bf9113..0000000 --- a/src/UniformCartesian/Interpolate.c +++ /dev/null @@ -1,609 +0,0 @@ -/*@@ - @file Interpolate.c - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - Interpolation of arrays to arbitrary points - - This interpolator is based on the Cactus 3.x Fortran version - written by Paul Walker. It also contains some nice optimization - features from Erik Schnetter. Jonathan Thornburg added some - additional comments in October 2001. - @enddesc - - @history - @date Wed 17 Jan 2001 - @author Thomas Radke - @hdesc Translation from Fortran to C - @date Thu 18 Oct 2001 - @author Jonathan Thornburg - @hdesc Add lots of comments, LOCALINTERP_VERBOSE_DEBUG debugging code - @date 22 Jan 2002 - @author Jonathan Thornburg - @hdesc Move all local-interpolation code from LocalInterp to here - @endhistory - - @version $Header$ - @@*/ - -#include <math.h> -#include <stdlib.h> -#include <string.h> -#include <stdio.h> - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "Interpolate.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Interpolate_c) - -#ifdef LOCALINTERP_VERBOSE_DEBUG - /* if this is >= 0, we print verbose debugging information at this point */ - int LocalInterp_verbose_debug_n = -1; -#endif - -/* the highest order of interpolation we support so far */ -#define MAXORDER 3 - -/* the highest dimension for variables we can deal with (so far) */ -#define MAXDIM 3 - -/******************************************************************************/ - -/*@@ - @routine INTERPOLATE (macro) - @date 18 Oct 2001 - @author code by ???, these comments by Jonathan Thornburg - @desc - This macro does the interpolation of in_array[] to compute - a single value out_array[n] (actually out_array[n]subpart ; - see the comments below for details). The data to be interpolated - must be real numbers of some type, i.e. if the arrays are - complex this macro must be called separately for the real - and imaginary parts. - @enddesc - - @var cctk_type - @vdesc C type of input and output array elements (might be complex) - @endvar - - @var cctk_subtype - @vdesc C type of actual numbers being interpolated (must be real) - @endvar - - @var subpart - @vdesc string to be suffixed to input/output array element to get - to get real number, i.e. empty string if cctk_type is real, - .Re or .Im as appropriate if cctk_type is complex - @endvar - - @var in_array - @vdesc A pointer to array to be interpolated (strictly speaking, to - the array's [0][0]...[0] element); this is typically passed - as a void * pointer so we typecast it as necessary. - @endvar - - @var out_array - @vdesc A 1-dimensional array where the interpolation result should be - stored; this is typically passed as a void * pointer so we - typecast it as necessary. - @endvar - - @var order - @vdesc The order of the interpolation (1=linear, 2=quadratic, 3=cubic, ...) - @endvar - - @var point - @vdesc [MAXDIM] array of integers giving the integer grid coordinates - of the closest grid point to the interpolation point; the - interpolation stencil/molecule is centered at this point. - @endvar - - @var dims - @vdesc [MAXDIM] array of integers giving the dimensions of in_array . - @endvar - - @var n - @vdesc Position in out_array where we should store the interpolation - result. - @endvar - - @var coeff - @vdesc [MAXDIM][MAX_ORDER+1] array of (floating-point) interpolation - coefficients; detailed semantics are that coeff[axis][m] is the - coefficient of y[m] when the 1-dimensional Lagrange interpolation - polynomial passing through the order+1 points - {(0,y[0]), (1,y[1]), ..., (order,y[order])} - is evaluated at the position x=offset[axis]. - @endvar - @@*/ -/* - * The basic idea here is that conceptually we first interpolate the - * (say) 3D gridfn in the x direction at each y and z grid point, - * then interpolate that 2D plane of values in the y direction at - * each z grid point, and finally interpolate that 1D line of values - * in the z direction. The implementation actually interleaves the - * different directions' interpolations so that only 3 scalar temporaries - * are needed. - */ -#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array, \ - order, point, dims, n, coeff) \ - { \ - int ii, jj, kk; \ - const cctk_type *fi; \ - /* for some reason (probably loop unrolling & pipelining) the compiler \ - produces faster code if arrays are used instead of scalars */ \ - cctk_subtype interp_result, fj[1], fk[1]; \ - \ - \ - interp_result = 0; \ - \ - /* NOTE-MAXDIM: support >3D arrays by adding more loops */ \ - for (kk = 0; kk <= order; kk++) \ - { \ - fk[0] = 0; \ - for (jj = 0; jj <= order; jj++) \ - { \ - /* NOTE-MAXDIM: for >3D arrays adapt the index calculation here */ \ - fi = (const cctk_type *) in_array + \ - point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk)); \ - \ - fj[0] = 0; \ - for (ii = 0; ii <= order; ii++) \ - { \ - fj[0] += fi[ii]subpart * coeff[0][ii]; \ - } \ - /* at this point we have just computed */ \ - /* fj[0] = in_array[*][jj][kk] interpolated to x=offset[0] */ \ - \ - fk[0] += fj[0] * coeff[1][jj]; \ - } \ - /* at this point we have just computed */ \ - /* fk[0] = fj[0][*][kk] interpolated to y=offset[1] */ \ - /* = in_array[*][*][kk] interpolated to */ \ - /* x=offset[0], y=offset[1] */ \ - \ - interp_result += fk[0] * coeff[2][kk]; \ - } \ - /* at this point we have just computed */ \ - /* interp_result = fk[0][*] interpolated to z=offset[2] */ \ - /* = in_array[*][*][*] interpolated to */ \ - /* x=offset[0], y=offset[1], z=offset[2] */ \ - \ - /* assign the result */ \ - ((cctk_type *) out_array)[n]subpart = interp_result; \ - } /* end of macro */ - -/* this is needed by some preprocessors to pass into INTERPOLATE - as a dummy macro */ -#define NOTHING - -/******************************************************************************/ - -/*@@ - @routine LocalInterp_Interpolate - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - This routine interpolates a set of input arrays - to a set of output arrays (one-to-one) at arbitrary points - which are given by their coordinates and the underlying - regular, uniform grid. - - Current limitations of this implementation are: - - arrays up to three (MAXDIM) dimensions only can be handled - - interpolation orders up to three (MAXORDER) only are supported - - coordinates must be given as CCTK_REAL types - - input and output array types must be the same - (no type casting of interpolation results supported) - - Despite of these limitations, the code was programmed almost - generically in that it can easily be extended to support - higher-dimensional arrays or more interpolation orders. - Places where the code would need to be changed to do this, - are marked with NOTE-MAXDIM and/or NOTE-MAXORDER comments - as appropriate. - @enddesc - - @var num_points - @vdesc number of points to interpolate at - @vtype int - @vio in - @endvar - @var num_dims - @vdesc dimensionality of the input arrays - @vtype int - @vio in - @endvar - @var num_arrays - @vdesc number of input/output arrays - @vtype int - @vio in - @endvar - @var dims - @vdesc dimensions of the input arrays - @vtype int[ num_dims ] - @vio in - @endvar - @var coord - @vdesc coordinates to interpolate at - @vtype CCTK_REAL coord[ num_dims * num_points ] - @vio in - @endvar - @var origin - @vdesc origin of the underlying grid - @vtype CCTK_REAL origin[ num_dims ] - @vio in - @endvar - @var delta - @vdesc deltas of the underlying grid - @vtype CCTK_REAL delta[ num_dims ] - @vio in - @endvar - @var in_types - @vdesc CCTK variable types of input arrays - @vtype int in_types[ num_arrays ] - @vio in - @endvar - @var in_arrays - @vdesc list of input arrays - @vtype void *in_arrays[ num_arrays ] - @vio in - @endvar - @var out_types - @vdesc CCTK variable types of output arrays - @vtype int out_types[ num_arrays ] - @vio in - @endvar - @var out_arrays - @vdesc list of output arrays - @vtype void *out_arrays[ num_arrays ] - @vio out - @endvar - - @returntype int - @returndesc - 0 - successful interpolation - -1 - in case of any errors - @endreturndesc - @@*/ -int LocalInterp_Interpolate (int order, - int num_points, - int num_dims, - int num_arrays, - const int dims[], - const CCTK_REAL coord[], - const CCTK_REAL origin[], - const CCTK_REAL delta[], - const int in_types[], - const void *const in_arrays[], - const int out_types[], - void *const out_arrays[]) -{ - int retval; - int i, a, n, shift; -#if 0 - int out_of_bounds; -#endif - int max_dims[MAXDIM], point[MAXDIM]; - CCTK_REAL delta_inv[MAXDIM]; - CCTK_REAL below[MAXDIM]; - CCTK_REAL offset[MAXDIM]; - CCTK_REAL coeff[MAXDIM][MAXORDER + 1]; - - /* verify parameters and check against our restrictions */ - retval = -1; - if (num_dims < 1) - { - CCTK_WARN (1, "Number of dimensions must be positive"); - } - else if (num_dims > MAXDIM) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation of %d-dimensional arrays not implemented", - num_dims); - } - else if (order < 1) - { - CCTK_WARN (1, "Inperpolation order must be positive"); - } - else if (order > MAXORDER) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation order %d not implemented", order); - } - else if (num_points < 0) - { - CCTK_WARN (1, "Negative number of points given"); - } - else - { - retval = 0; - } - - /* immediately return in case of errors */ - if (retval) - { - return (retval); - } - - /* also immediately return if there's nothing to do */ - if (num_points == 0) - { - return (retval); - } - - /* avoid divisions by delta later on */ - for (i = 0; i < num_dims; i++) - { - delta_inv[i] = 1.0 / delta[i]; - } - - /* duplicate the dims[] vector into one with MAXDIM-1 elements - (with the remaining elements zeroed out) - so that we can use nested loops over MAXDIM dimensions later on */ - memset (max_dims, 0, sizeof (max_dims)); - memcpy (max_dims, dims, (num_dims - 1) * sizeof (int)); - - /* zero out the coefficients and set the elements with index 'order' to one - so that we can use nested loops over MAXDIM dimensions later on */ - memset (coeff, 0, sizeof (coeff)); - for (i = num_dims; i < MAXDIM; i++) - { - coeff[i][0] = 1; - } - - /* zero out the iterator */ - memset (point, 0, sizeof (point)); - - /* loop over all points to interpolate at */ - for (n = 0; n < num_points; n++) - { -#if 0 - /* reset the out-of-bounds flag */ - out_of_bounds = 0; -#endif - - /* loop over all dimensions */ - for (i = 0; i < num_dims; i++) - { - /* closest grid point for stencil/molecule */ - point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i] - - 0.5 * (order - 1)); - -#if 0 - /* test bounds */ - out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; -#endif - - /* check that the stencil/molecule isn't bigger than the grid */ - if (order+1 > dims[i]) /* stencil/molecule size = order+1 */ - { - return -1; - } - - /* if beyond lower bound shift the grid point to the right */ - shift = point[i]; - if (shift < 0) - { - point[i] -= shift; - } - - /* if beyond upper bound shift the grid point to the left */ - shift = point[i] + order - (dims[i] - 1); - if (shift > 0) - { - point[i] -= shift; - } - - /* physical coordinate of that grid point */ - below[i] = origin[i] + point[i] * delta[i]; - - /* offset from that grid point, in fractions of grid points */ - offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i]; - } - -#ifdef LOCALINTERP_VERBOSE_DEBUG -if (n == LocalInterp_verbose_debug_n) - { - int ii; - printf("out_of_bounds = %d\n", out_of_bounds); - for (ii = 0 ; ii < num_dims ; ++ii) - { - printf("offset[%d] = %g\n", ii, (double) offset[ii]); - } - } -#endif /* LOCALINTERP_VERBOSE_DEBUG */ - -#if 0 - /* check bounds */ - if (out_of_bounds) - { - char *msg; - - - /* put all information into a single message string for output */ - msg = (char *) malloc (100 + num_dims*(10 + 4*20)); - sprintf (msg, - "Interpolation stencil/molecule out of bounds at grid point [%d", - point[0]); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %d", msg, point[i]); - } - sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg, - (double) below[0], (double) (below[0] + offset[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double) below[i], (double) (below[i] + offset[i])); - } - sprintf (msg, "%s]\ngrid is min/max [%f / %f", msg, - (double) origin[0], (double) (origin[0] + (dims[0]-1)*delta[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i])); - } - sprintf (msg, "%s]", msg); - CCTK_WARN (1, msg); - free (msg); - - continue; - } -#endif - - /* - * *** compute the interpolation coefficients according to the order *** - * - * (Thanks to Erik for formulating these so nicely.) - * - * These formulas are "just" the coefficients of the classical - * Lagrange interpolation polynomials along each dimension. - * For example, in 1 dimension the unique quadratic passing - * through the 3 points {(x0,y0), (x1,y1), (x2,y2)} is: - * ( x-x1)( x-x2) ( x-x0)( x-x2) ( x-x0)( x-x1) - * -------------- y0 + -------------- y1 + -------------- y2 - * (x0-x1)(x0-x2) (x1-x0)(x1-x2) (x2-x0)(x2-x1) - * (It's easy to see this: each of the terms is yi if x=xi, or - * zero if x=any other xj.) To get the formulas below, just negate - * each (x-x) factor, and substitute the values xi=i. - */ - /* - * NOTE-MAXORDER: support higher interpolation orders by adding the - * appropriate coefficients in another else branch - */ - if (order == 1) - { - /* first order (linear) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = 1 - offset[i]; - coeff[i][1] = offset[i]; - } - } - else if (order == 2) - { - /* second order (quadratic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) / ( 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) / ( 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2)); - } - } - else if (order == 3) - { - /* third order (cubic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 3 * 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 2 * 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) / - ( 1 * (-1) * (-2)); - coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) / - ((-1) * (-2) * (-3)); - } - } - else - { - CCTK_WARN (0, "Implementation error"); - } - -#ifdef LOCALINTERP_VERBOSE_DEBUG -if (n == LocalInterp_verbose_debug_n) - { - int ii,mm; - for (ii = 0 ; ii < num_dims ; ++ii) - { - for (mm = 0 ; mm <= order ; ++mm) - { - printf("coeff[%d][%d] = %g\n", - ii, mm, (double) coeff[ii][mm]); - } - } - } -#endif /* LOCALINTERP_VERBOSE_DEBUG */ - - /* now loop over all arrays to interpolate at the current point */ - for (a = 0; a < num_arrays; a++) - { - /* sorry, for now input and output arrays must be of same type */ - if (in_types[a] != out_types[a]) - { - CCTK_WARN (1, "Type casting of interpolation results not implemented"); - continue; - } - - /* now do the interpolation according to the array type - we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */ - if (in_types[a] == CCTK_VARIABLE_REAL) - { - INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX) - { - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#ifdef CCTK_REAL4 - else if (in_types[a] == CCTK_VARIABLE_REAL4) - { - INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX8) - { - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL8 - else if (in_types[a] == CCTK_VARIABLE_REAL8) - { - INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX16) - { - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL16 - else if (in_types[a] == CCTK_VARIABLE_REAL16) - { - INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX32) - { - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif - else - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Unsupported variable type %d", (int) in_types[a]); - } - } /* end of loop over all arrays */ - - } /* end of loop over all points to interpolate at */ - - /* we're done */ - return (retval); -} diff --git a/src/UniformCartesian/Interpolate.h b/src/UniformCartesian/Interpolate.h deleted file mode 100644 index 35dac21..0000000 --- a/src/UniformCartesian/Interpolate.h +++ /dev/null @@ -1,45 +0,0 @@ -/*@@ - @file Interpolate.h - @date 22 Jan 2002 - @author Jonathan Thornburg - @desc function prototypes - @enddesc - @version $Header$ - @history - @date 22 Jan 2002 - @author Jonathan Thornburg - @hdesc created by editing down LocalInterp::pughInterpGH.h - (originally dated 4 July 1999, by Thomas Radke) - @endhistory - @@*/ - -/* prototypes of interpolation operator to be registered */ -int LocalInterp_InterpLocal(cGH *GH, - int order, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]); - -/* prototypes of routines used internally */ -int LocalInterp_Interpolate(int order, - int num_points, - int num_dims, - int num_arrays, - const int dims[ /* num_dims */ ], - const CCTK_REAL coord[ /* num_dims*num_points */ ], - const CCTK_REAL origin[ /* num_dims */ ], - const CCTK_REAL delta[ /* num_dims */ ], - const int in_types[ /* num_arrays */ ], - const void *const in_arrays[ /* num_arrays */ ], - const int out_types[ /* num_arrays */ ], - void *const out_arrays[ /* num_arrays */ ]); diff --git a/src/UniformCartesian/Operator.c b/src/UniformCartesian/Operator.c deleted file mode 100644 index 28de748..0000000 --- a/src/UniformCartesian/Operator.c +++ /dev/null @@ -1,323 +0,0 @@ -/*@@ - @file Operator.c - @date Tue Apr 15 18:22:45 1997 - @author Paul Walker - @desc - Definition of interpolation operators for regular uniform grids. - @enddesc - - @history - @date Sun Jul 04 1999 - @author Thomas Radke - @hdesc conversion to Cactus 4.0 (copied from pughGetPoints.c) - @date Wed 31 Jan 2001 - @author Thomas Radke - @hdesc translation of fortran interpolators into C - @date 22 Jan 2002 - @author Jonathan Thornburg - @hdesc Move all local-interpolation code from LocalInterp to here - @endhistory - - @version $Header$ - @@*/ - -#include <stdlib.h> -#include <math.h> -#include <string.h> - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "CactusPUGH/PUGH/src/include/pugh.h" -#include "Interpolate.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Operator_c) - -/* uncomment this to get some debugging output */ -/* #define PUGHINTERP_DEBUG 1 */ - - - -/* prototypes of routines defined in this source file */ -static int LocalInterp_CheckArguments (cGH *GH, - int num_dims, - int num_points, - int num_in_arrays, - int num_out_arrays, - const int interp_coord_array_types[]); - -/******************************************************************************/ - -/*@@ - @routine LocalInterp_InterpLocal - @date Sun Jul 04 1999 - @author Thomas Radke - @desc - The interpolation operator registered with the CCTK - under the name "uniform cartesian". - - Interpolates a list of non-distributed, processor-local - input arrays to a list of output arrays (one-to-one) - at a given number of interpolation points (indicated by - their coordinates). The points are located on a coordinate - system which is assumed to be a uniform cartesian. - @enddesc - - @var GH - @vdesc Pointer to CCTK grid hierarchy - @vtype cGH * - @vio in - @endvar - @var order - @vdesc interpolation order - @vtype int - @vio in - @endvar - @var num_points - @vdesc number of points to be interpolated on this processor - @vtype int - @vio in - @endvar - @var num_dims - @vdesc dimensionality of the underlying grid - @vtype int - @vio in - @endvar - @var num_in_arrays - @vdesc number of input arrays to interpolate from - @vtype int - @vio in - @endvar - @var num_out_arrays - @vdesc number of output arrays to interpolate to - @vtype int - @vio in - @endvar - @var coord_dims - @vdesc dimensions of the underlying grid - @vtype int [num_dims] - @vio in - @endvar - @var coord_arrays - @vdesc list of grid coordinate arrays - @vtype void *const [num_dims] - @vio in - @endvar - @var coord_array_types - @vdesc CCTK data type of coordinate arrays - @vtype int [num_dims] - @vio int - @endvar - @var interp_coord_arrays - @vdesc coordinates of points to interpolate at - @vtype void *const [num_dims] - @vio in - @endvar - @var interp_coord_array_types - @vdesc CCTK data type of coordinate arrays to interpolate at - @vtype int [num_dims] - @vio out - @endvar - @var in_arrays - @vdesc list of input arrays to interpolate on - @vtype void *const [num_in_arrays] - @vio in - @endvar - @var in_array_types - @vdesc CCTK data types of input arrays - @vtype int [num_in_arrays] - @vio in - @endvar - @var out_arrays - @vdesc list of output arrays to interpolate to - @vtype void *const [num_out_arrays] - @vio out - @endvar - @var out_array_types - @vdesc CCTK data types of output arrays - @vtype int [num_out_arrays] - @vio in - @endvar - - @returntype int - @returndesc - 0 - successful interpolation - -1 - in case of any errors - @endreturndesc - @@*/ -int LocalInterp_InterpLocal (cGH *GH, - int order, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]) -{ - int dim, point, retval; - CCTK_REAL *coords, *origin, *delta; - const CCTK_REAL *const *data; - - /* check arguments */ - retval = LocalInterp_CheckArguments (GH, num_dims, num_points, - num_in_arrays, num_out_arrays, - interp_coord_array_types); - if (retval <= 0) - { - return (retval); - } - for (dim = 0; dim < num_dims; dim++) - { - if (coord_array_types[dim] != CCTK_VARIABLE_REAL) - { - CCTK_WARN (1, "Coordinates should be of type CCTK_REAL"); - return (-1); - } - if (interp_coord_array_types[dim] != CCTK_VARIABLE_REAL) - { - CCTK_WARN (1, "Interpolation coordinates should be of type CCTK_REAL"); - return (-1); - } - } - - /* get the grid spacings - this assumes a cartesian grid */ - origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL)); - delta = origin + num_dims; - data = (const CCTK_REAL *const *) coord_arrays; - for (dim = 0; dim < num_dims; dim++) - { - origin[dim] = data[dim][0]; - delta[dim] = data[dim][1] - data[dim][0]; - } - - /* sort the individual interpolation coordinate arrays into a single one */ - coords = (CCTK_REAL *) malloc (num_dims * num_points * sizeof (CCTK_REAL)); - data = (const CCTK_REAL *const *) interp_coord_arrays; - for (point = 0; point < num_points; point++) - { - for (dim = 0; dim < num_dims; dim++) - { - *coords++ = data[dim][point]; - } - } - coords -= num_dims * num_points; - - /* call the interpolator function */ - retval = LocalInterp_Interpolate (order, - num_points, num_dims, num_out_arrays, - coord_dims, coords, origin, delta, - in_array_types, in_arrays, - out_array_types, out_arrays); - - /* free allocated resources */ - free (coords); - free (origin); - - return (retval); -} - -/**************************************************************************/ -/* local routines */ -/**************************************************************************/ - -/*@@ - @routine LocalInterp_CheckArguments - @date Thu 25 Jan 2001 - @author Thomas Radke - @desc - Checks the interpolation arguments passed in via - the flesh's general interpolation calling interface - - This routine also verifies that the parameters meet - the limitations of LocalInterp's interpolation operators. - @enddesc - - @var GH - @vdesc Pointer to CCTK grid hierarchy - @vtype cGH * - @vio in - @endvar - @var num_dims - @vdesc dimensionality of the underlying grid - @vtype int - @vio in - @endvar - @var num_points - @vdesc number of points to interpolate at - @vtype int - @vio in - @endvar - @var num_in_arrays - @vdesc number of passed input arrays - @vtype int - @vio in - @endvar - @var num_out_arrays - @vdesc number of passed input arrays - @vtype int - @vio in - @endvar - @var interp_coord_array_types - @vdesc types of passed coordinates to interpolate at - @vtype int [num_dims] - @vio in - @endvar - - @returntype int - @returndesc - +1 for success - 0 for success but nothing to do - -1 for failure (wrong parameters passed or limitations not met) - @endreturndesc - @@*/ -static int LocalInterp_CheckArguments (cGH *GH, - int num_dims, - int num_points, - int num_in_arrays, - int num_out_arrays, - const int interp_coord_array_types[]) -{ - int i; - - - /* check for invalid arguments */ - if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0) - { - return (-1); - } - - /* check if there's anything to do at all */ - /* NOTE: num_points can be 0 in a collective call */ - if (num_dims == 0 || num_in_arrays == 0 || num_out_arrays == 0) - { - return (0); - } - - /* for now we can only deal with coordinates of type CCTK_REAL */ - for (i = 0; i < num_dims; i++) - { - if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL) - { - CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL"); - return (-1); - } - } - - /* LocalInterp's interpolation operators compute one output array - per input array */ - if (num_in_arrays != num_out_arrays) - { - CCTK_WARN (1, "Number of input arrays must match number of output arrays"); - return (-1); - } - - return (1); -} diff --git a/src/UniformCartesian/README b/src/UniformCartesian/README deleted file mode 100644 index 2342312..0000000 --- a/src/UniformCartesian/README +++ /dev/null @@ -1,14 +0,0 @@ -This directory contains the interpolation functions registered for - CCTK_InterpLocal() -under the names - "first-order uniform cartesian" - "second-order uniform cartesian" - "third-order uniform cartesian" - -Implementation Notes -==================== -The interpolation operators registered for different orders are mapped -via wrappers (in "Startup.c") onto a single routine (in "Operator.c"), -just passing the order as an additional argument. - -The actual core interpolation routine is located in "Interpolate.c". diff --git a/src/UniformCartesian/Startup.c b/src/UniformCartesian/Startup.c deleted file mode 100644 index b25dc9f..0000000 --- a/src/UniformCartesian/Startup.c +++ /dev/null @@ -1,185 +0,0 @@ -/*@@ - @file Startup.c - @date Sun Jul 04 1999 - @author Thomas Radke - @desc - Startup routines for LocalInterp/UniformCartesian - @enddesc - - @history - @date 22 Jan 2002 - @author Jonathan Thornburg - @hdesc Move all local-interpolation code from LocalInterp to here - @endhistory - - @version $Header$ - @@*/ - -#include <stdlib.h> - -#include "cctk.h" -#include "cctk_Interp.h" -#include "Interpolate.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Startup_c) - - -/* prototypes of externally-visible routines defined in this source file */ -void LocalInterp_UC_Startup(void); - -/* prototypes of static routines defined in this source file */ -static int LocalInterp_InterpLocal_1stOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]); -static int LocalInterp_InterpLocal_2ndOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]); -static int LocalInterp_InterpLocal_3rdOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]); - -/******************************************************************************/ - -/*@@ - @routine LocalInterp_InterpLocal_NthOrder - @date Wed 14 Feb 2001 - @author Thomas Radke - @desc - Wrappers for the different interpolation operators - registered for first/second/third order interpolation. - These wrappers just call the common interpolation routine - passing all arguments plus the interpolation order. - @enddesc - - @returntype int - @returndesc - the return code of the common interpolation routine - @endreturndesc - @@*/ -static int LocalInterp_InterpLocal_1stOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]) -{ - return (LocalInterp_InterpLocal (GH, 1, num_points, num_dims, - num_in_arrays, num_out_arrays, - coord_dims, coord_arrays, coord_array_types, - interp_coord_arrays, interp_coord_array_types, - in_arrays, in_array_types, - out_arrays, out_array_types)); -} - - -static int LocalInterp_InterpLocal_2ndOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]) -{ - return (LocalInterp_InterpLocal (GH, 2, num_points, num_dims, - num_in_arrays, num_out_arrays, - coord_dims, coord_arrays, coord_array_types, - interp_coord_arrays, interp_coord_array_types, - in_arrays, in_array_types, - out_arrays, out_array_types)); -} - - -static int LocalInterp_InterpLocal_3rdOrder(cGH *GH, - int num_points, - int num_dims, - int num_in_arrays, - int num_out_arrays, - const int coord_dims[], - const void *const coord_arrays[], - const int coord_array_types[], - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const void *const in_arrays[], - const int in_array_types[], - void *const out_arrays[], - const int out_array_types[]) -{ - return (LocalInterp_InterpLocal (GH, 3, num_points, num_dims, - num_in_arrays, num_out_arrays, - coord_dims, coord_arrays, coord_array_types, - interp_coord_arrays, interp_coord_array_types, - in_arrays, in_array_types, - out_arrays, out_array_types)); -} - -/******************************************************************************/ - -/*@@ - @routine LocalInterp_Startup - @date Sun Jul 04 1999 - @author Thomas Radke - @desc - The startup registration routine for LocalInterp. - Registers the interpolation operators with the flesh. - @enddesc - @calls CCTK_InterpRegisterOperatorLocal - @@*/ -void LocalInterp_UC_Startup(void) -{ - CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_1stOrder, - "first-order uniform cartesian"); - CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_2ndOrder, - "second-order uniform cartesian"); - CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_3rdOrder, - "third-order uniform cartesian"); -} diff --git a/src/UniformCartesian/make.code.defn b/src/UniformCartesian/make.code.defn deleted file mode 100644 index 4a506e8..0000000 --- a/src/UniformCartesian/make.code.defn +++ /dev/null @@ -1,4 +0,0 @@ -# $Header$ - -# Source files in this directory -SRCS = Startup.c Operator.c Interpolate.c |