diff options
Diffstat (limited to 'src/CoordinateStuff.c')
| -rw-r--r-- | src/CoordinateStuff.c | 165 |
1 files changed, 165 insertions, 0 deletions
diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c new file mode 100644 index 0000000..59e99a1 --- /dev/null +++ b/src/CoordinateStuff.c @@ -0,0 +1,165 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "cctk.h" + + /*@@ + @routine IndexFloor + @date Fri Jan 7 10:34:29 2000 + @author Gerd Lanfermann + @desc + For a given physical, global coordinate, IndexFloor returns + the closest lower grid index *locally* in the direction d, + which is owned by the grid path. (Ghostzone are not owned!) + Other return values: + -1 : gridindex not on this gridpatch + -3 : coordinate is outside global grid + + This routine and the Fortran wrapper are Fortran/C index aware. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +int IndexFloor(cGH *GH, CCTK_REAL coord_value, int d) +{ + int ierr,index_low,idummy,ugs,lgs; + char *message; + CCTK_REAL cmin,cmax; + + if (d==0) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x"); + else if (d==1) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y"); + else if (d==2) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"z"); + else + CCTK_WARN(0,"IndexFloor: dimension is not valid. Fortran:1..3 C:0..2"); + + if ((coord_value<cmin)||(coord_value>cmax)) { + message = (char*)malloc(128*sizeof(char)); + sprintf(message, + "IndexFloor: coordinate value outside grid [%f,%f]: %f \n", + cmin, cmax, coord_value); + CCTK_WARN(2,message); + free(message); + return(-3); + } + + idummy = floor((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]); + index_low = idummy-GH->cctk_lbnd[d]; + + if (GH->cctk_bbox[2*d]==1) + lgs = 0; + else + lgs = GH->cctk_nghostzones[d]; + + if (GH->cctk_bbox[2*d+1]==1) + ugs = 0; + else + ugs = GH->cctk_nghostzones[d]; + + if (index_low<lgs) + index_low = -1; + + if (index_low>=GH->cctk_lsh[d]-ugs) + index_low = -1; + + return(index_low); +} + + + /*@@ + @routine IndexCeil + @date Fri Jan 7 10:34:29 2000 + @author Gerd Lanfermann + @desc + For a given physical, global coordinate, IndexCeil returns + the closest upper grid index *locally* in the direction d, + which is owned by the grid path. (Ghostzone are not owned!) + Other return values: + -1 : grid index not on this gridpatch + -3 : coordinate is outside global grid + + This routine and the Fortran wrapper are Fortran/C index aware. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +int IndexCeil(cGH *GH, CCTK_REAL coord_value, int d) +{ + int ierr,index_up,idummy,ugs,lgs; + char *message; + CCTK_REAL cmin,cmax; + + if (d==0) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x"); + else if (d==1) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y"); + else if (d==2) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"z"); + else + CCTK_WARN(0,"IndexCeil: dimension is not valid. Fortran:1..3 C:0..2"); + + + if ((coord_value<cmin)||(coord_value>cmax)) { + message = (char*) malloc(256*sizeof(char)); + sprintf(message, + "IndexCeil: coordinate value outside grid [%f,%f]: %f \n", + cmin,cmax, coord_value); + CCTK_WARN(2,message); + free(message); + return(-3); + } + + idummy = ceil((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]); + index_up = idummy-GH->cctk_lbnd[d]; + + if (GH->cctk_bbox[2*d]==1) + lgs = 0; + else + lgs = GH->cctk_nghostzones[d]; + + if (GH->cctk_bbox[2*d+1]==1) + ugs = 0; + else + ugs = GH->cctk_nghostzones[d]; + + if (index_up<lgs) + index_up = -1; + + if (index_up>=GH->cctk_lsh[d]-ugs) + index_up = -1; + + return(index_up); +} + +void FMODIFIER FORTRAN_NAME(IndexFloor) + (cGH *GH, CCTK_REAL *coord_value, int *index_low, int *fdim) { + + /* note the increase/decrease of indeces for fortran compatability */ + *index_low = IndexFloor(GH, *coord_value, (*fdim)-1); + if ((*index_low)>=0) (*index_low)++; +} + +void FMODIFIER FORTRAN_NAME(IndexCeil) + (cGH *GH, CCTK_REAL *coord_value, int *index_up, int *fdim) { + /* note the increase/decrease of indeces for fortran compatability */ + *index_up = IndexCeil(GH, *coord_value, (*fdim)-1); + if ((*index_up)>=0) (*index_up)++; +} + + |