From d5ec2d9d47162065b3cca9dda8a7b8a08f453d1f Mon Sep 17 00:00:00 2001 From: schnetter Date: Tue, 20 Oct 2009 19:11:36 +0000 Subject: Support interpolating to obtain the image values. This makes it possible to use IOJpeg together with mesh refinement. There is a new keyword parameter "gridpoints" with possible values "hyperslab" and "interpolate". Hyperslabbing works only on uniform grids, whereas interpolation also works with AMR. Additional parameters select where in the domain the points should be interpolated. git-svn-id: http://svn.cactuscode.org/arrangements/CactusIO/IOJpeg/trunk@130 eff87b29-5268-4891-90a3-a07138403961 --- src/ChooseOutput.c | 12 +- src/Output.c | 2 +- src/Startup.c | 6 +- src/Write.c | 341 +++++++++++++++++++++++++++++++++++++++++------------ 4 files changed, 274 insertions(+), 87 deletions(-) (limited to 'src') diff --git a/src/ChooseOutput.c b/src/ChooseOutput.c index 7366fa7..04497c5 100644 --- a/src/ChooseOutput.c +++ b/src/ChooseOutput.c @@ -13,6 +13,7 @@ #include #include "cctk.h" +#include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "CactusBase/IOUtil/src/ioutil_Utils.h" #include "ioJpegGH.h" @@ -41,8 +42,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_ChooseOutput_c) /******************************************************************** ******************** External Routines ************************ ********************************************************************/ -void IOJpeg_ChooseOutput (const cGH *GH); - /*@@ @routine IOJpeg_ChooseOutput @@ -60,13 +59,14 @@ void IOJpeg_ChooseOutput (const cGH *GH); @vio in @endvar @@*/ -void IOJpeg_ChooseOutput (const cGH *GH) +void IOJpeg_ChooseOutput (CCTK_ARGUMENTS) { int i, maxdim; ioJpegGH *myGH; int origin_index[3]; CCTK_REAL origin_phys[3]; - DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; GET_SLICE (out2D_xyplane_z, out_xyplane_z, origin_index[0], origin_phys[0]); @@ -74,7 +74,7 @@ void IOJpeg_ChooseOutput (const cGH *GH) GET_SLICE (out2D_yzplane_x, out_yzplane_x, origin_index[2], origin_phys[2]); maxdim = CCTK_MaxDim (); - myGH = (ioJpegGH *) CCTK_GHExtension (GH, "IOJpeg"); + myGH = (ioJpegGH *) CCTK_GHExtension (cctkGH, "IOJpeg"); myGH->sp2xyz = (int **) malloc (3 * sizeof (int *)); for (i = 0; i < maxdim; i++) @@ -83,7 +83,7 @@ void IOJpeg_ChooseOutput (const cGH *GH) if (i > 0 && i < 3) { - IOUtil_2DPlanes (GH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]); + IOUtil_2DPlanes (cctkGH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]); } } } diff --git a/src/Output.c b/src/Output.c index 312e4ca..ea72af8 100644 --- a/src/Output.c +++ b/src/Output.c @@ -377,7 +377,7 @@ static void SetOutputFlag (int vindex, const char *optstring, void *arg) IOUtil_ParseOutputFrequency (CCTK_THORNSTRING, "IOJpeg::out_vars", myGH->stop_on_parse_errors, vindex, optstring, - &myGH->out_every[vindex]); + &myGH->out_every[vindex], NULL); } } } diff --git a/src/Startup.c b/src/Startup.c index 27566d8..c90472b 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -27,8 +27,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_Startup_c) /******************************************************************** ******************** External Routines ************************ ********************************************************************/ -void IOJpeg_Startup (void); - /******************************************************************** ******************** Internal Routines ************************ @@ -49,13 +47,15 @@ static void *SetupGH (tFleshConfig *config, int conv_level, cGH *GH); @calls CCTK_RegisterGHExtension CCTK_RegisterGHExtensionSetupGH @@*/ -void IOJpeg_Startup (void) +int IOJpeg_Startup (void) { if (CCTK_MaxDim () >= 2) { CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("IOJpeg"), SetupGH); } + + return 0; } diff --git a/src/Write.c b/src/Write.c index 3fa5c7d..4ada69b 100644 --- a/src/Write.c +++ b/src/Write.c @@ -9,6 +9,9 @@ @@*/ +#include +#include +#include #include #include @@ -118,103 +121,287 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias) min = colormap_min; max = colormap_max; } - else + else if (CCTK_Equals (colormap, "auto")) { i = CCTK_LocalArrayReductionHandle ("minimum"); - i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min); + CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min); i = CCTK_LocalArrayReductionHandle ("maximum"); - i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max); + CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max); + } + else if (CCTK_Equals (colormap, "auto-old")) + { + i = CCTK_ReductionHandle ("minimum"); + CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &min, 1, vindex); + i = CCTK_ReductionHandle ("maximum"); + CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &max, 1, vindex); + } + else + { + CCTK_WARN (CCTK_WARN_ABORT, "Setting for parameter IOpjeg::colormap not recognised"); } /* get the number of slices to output */ /* in general: maxdir = gdata.dim * (gdata.dim - 1) / 2; */ maxdir = gdata.dim == 2 ? 1 : 3; - /* now do the actual I/O looping over all directions */ - for (dir = 0; dir < maxdir; dir++) + if (CCTK_EQUALS(gridpoints, "hyperslab")) { - /* get the directions to span the hyperslab */ - if (dir == 0) + + /* now do the actual I/O looping over all directions */ + for (dir = 0; dir < maxdir; dir++) { - dir_i = 0; dir_j = 1; /* xy */ - } - else if (dir == 1) - { - dir_i = 0; dir_j = 2; /* xz */ - } - else - { - dir_i = 1; dir_j = 2; /* yz */ - } + /* get the directions to span the hyperslab */ + if (dir == 0) + { + dir_i = 0; dir_j = 1; /* xy */ + } + else if (dir == 1) + { + dir_i = 0; dir_j = 2; /* xz */ + } + else + { + dir_i = 1; dir_j = 2; /* yz */ + } + + /* set the origin using the slice center from IOUtil */ + memset (origin, 0, sizeof (origin)); + if (gdata.grouptype == CCTK_GF) + { + origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir]; + } + + /* set the direction vector */ + memset (direction, 0, sizeof (direction)); + direction[dir_i] = direction[gdata.dim + dir_j] = 1; + + mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2, + direction, origin, extent, + NULL, -1, NULL, hsize); + if (mapping < 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Failed to define hyperslab mapping for variable '%s'", + fullname); + continue; + } + { + // Sanity check + // (if this fails, we requested an insane number of grid points) + int max = INT_MAX; + int d; + for (d=0; d<2; ++d) { + assert (hsize[d] >= 0 && hsize[d] <= max); + if (hsize[d] > 0) max /= hsize[d]; + } + assert (CCTK_VarTypeSize (gdata.vartype) <= max); + } + total_hsize = hsize[0] * hsize[1] * CCTK_VarTypeSize (gdata.vartype); + assert (total_hsize >= 0); + if (total_hsize <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Selected hyperslab has zero size for variable '%s' " + "direction %d", fullname, dir); + Hyperslab_FreeMapping (mapping); + continue; + } + + /* allocate hyperslab buffer */ + hdata = myproc == 0 ? malloc (total_hsize) : NULL; + + /* get the hyperslab */ + i = Hyperslab_Get (GH, mapping, 0, vindex, 0, htype, hdata); + if (i) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Failed to extract hyperslab for variable '%s'", fullname); + } + + /* release the mapping structure */ + Hyperslab_FreeMapping (mapping); + + /* and dump the data to file */ + if (myproc == 0) + { + if (i == 0) + { + /* refine the 2D slice if requested */ + outdata = refinement_factor > 1 ? RefineData (hsize, hdata) : hdata; + if (! outdata) + { + CCTK_WARN (1, "Failed to refine 2D array"); + outdata = hdata; + } + WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata); + if (outdata != hdata) + { + free (outdata); + } + } + + /* clean up */ + free (hdata); + } + } /* end of looping through xyz directions */ - /* set the origin using the slice center from IOUtil */ - memset (origin, 0, sizeof (origin)); - if (gdata.grouptype == CCTK_GF) + } + else if (CCTK_EQUALS(gridpoints, "interpolate")) + { + + int interpolator; + int options_table; + int coord_handle; + + CCTK_REAL * restrict coordsx; + CCTK_REAL * restrict coordsy; + CCTK_REAL * restrict coordsz; + CCTK_POINTER_TO_CONST coords[3]; + CCTK_INT * restrict inputs; + CCTK_INT * restrict output_types; + CCTK_POINTER * restrict outputs; + CCTK_INT * restrict operation_codes; + CCTK_INT * restrict time_deriv_order; + + const int nvars = 1; + int npoints; + + int n; + int i, j, k; + int d; + int ierr; + + CCTK_REAL *outdata; + + + + interpolator = CCTK_InterpHandle (interpolator_name); + assert (interpolator >= 0); + + options_table = Util_TableCreateFromString (interpolator_options); + assert (options_table >= 0); + + coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates); + assert (coord_handle >= 0); + + + { - origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir]; + ierr = Util_TableSetInt (options_table, 1, "want_global_mode"); + assert (! ierr); } - - /* set the direction vector */ - memset (direction, 0, sizeof (direction)); - direction[dir_i] = direction[gdata.dim + dir_j] = 1; - - mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2, - direction, origin, extent, - NULL, -1, NULL, hsize); - if (mapping < 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Failed to define hyperslab mapping for variable '%s'", - fullname); - continue; + + + + /* 2D Arrays */ + npoints = myproc==0 ? array2d_npoints_i * array2d_npoints_j : 0; + + coordsx = malloc (npoints * sizeof * coordsx); + assert (npoints==0 || coordsx); + coordsy = malloc (npoints * sizeof * coordsy); + assert (npoints==0 || coordsy); + coordsz = malloc (npoints * sizeof * coordsz); + assert (npoints==0 || coordsz); + coords[0] = coordsx; + coords[1] = coordsy; + coords[2] = coordsz; + + if (myproc==0) { + n = 0; + for (j=0; j 1 ? RefineData (hsize, hdata) : hdata; - if (! outdata) - { - CCTK_WARN (1, "Failed to refine 2D array"); - outdata = hdata; - } - WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata); - if (outdata != hdata) - { - free (outdata); - } + + outputs = malloc (nvars * sizeof * outputs); + assert (outputs); + + outdata = malloc (npoints * sizeof * outdata); + assert (npoints==0 || outdata); + + for (n=0; n