aboutsummaryrefslogtreecommitdiff
path: root/src/Write3D.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Write3D.c')
-rw-r--r--src/Write3D.c292
1 files changed, 292 insertions, 0 deletions
diff --git a/src/Write3D.c b/src/Write3D.c
new file mode 100644
index 0000000..f07ef55
--- /dev/null
+++ b/src/Write3D.c
@@ -0,0 +1,292 @@
+/*@@
+ @file Write3D.c
+ @date Sat 12 June 2004
+ @author Thomas Radke
+ @desc
+ Three-dimensional output of variables in SDF file format.
+ @enddesc
+ @version $Id$
+ @@*/
+
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "util_Table.h"
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "CactusBase/IOUtil/src/ioutil_AdvertisedFiles.h"
+#include "ioSDFGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusIO_IOSDF_Write3D_c)
+
+
+/********************************************************************
+ ******************** Internal Routines ************************
+ ********************************************************************/
+static char *OpenFile (const cGH *GH, const char *fullname, const char *alias);
+
+/*@@
+ @routine IOSDF_Write3D
+ @date Sat 12 June 2004
+ @author Thomas Radke
+ @desc
+ Writes the 3D volume of a variable into a gnuplot SDF file.
+ @enddesc
+ @calls Hyperslab_GlobalMappingByIndex
+ Hyperslab_FreeMapping
+ Hyperslab_GetList
+ OpenFile
+
+ @var GH
+ @vdesc Pointer to CCTK GH
+ @vtype const cGH *
+ @vio in
+ @endvar
+ @var vindex
+ @vdesc index of variable to output
+ @vtype int
+ @vio in
+ @endvar
+ @var alias
+ @vdesc alias name of variable to output
+ @vtype const char *
+ @vio in
+ @endvar
+
+ @returntype int
+ @returndesc
+ 0 for success, or<BR>
+ -1 if variable has no storage assigned<BR>
+ -2 if output file couldn't be opened<BR>
+ -3 if hyperslab for coordinates and/or variable couldn't be
+ extracted
+ @endreturndesc
+@@*/
+int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
+{
+ int i, total_hsize;
+ int myproc, gindex, have_coords;
+ int num_requested_hslabs, num_returned_hslabs;
+ cGroup gdata;
+ CCTK_INT coord_system_handle, coord_handles[3];
+ char *fullname, *groupname, *filename;
+ double *hdata[4];
+ int extent_int[3];
+ double offset[3];
+ int mapping;
+ CCTK_INT vindices[4], extent[3], downsample[3], hsize[3];
+ const double dtime = GH->cctk_time;
+ const CCTK_INT origin[] = {0, 0, 0},
+ direction[] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
+ DECLARE_CCTK_PARAMETERS
+
+
+ /* get the variable group information */
+ fullname = CCTK_FullName (vindex);
+ gindex = CCTK_GroupIndexFromVarI (vindex);
+ CCTK_GroupData (gindex, &gdata);
+
+ /* check if variable has storage assigned */
+ if (! CCTK_QueryGroupStorageI (GH, gindex))
+ {
+ CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "No IOSDF 3D output for '%s' (no storage)", fullname);
+ free (fullname);
+ return (-1);
+ }
+
+ /* get the coordinate system associated with this grid variable */
+ vindices[0] = vindex;
+ groupname = CCTK_GroupName (gindex);
+ coord_system_handle = Coord_GroupSystem (GH, groupname);
+ free (groupname);
+
+ have_coords = coord_system_handle >= 0 &&
+ Util_TableGetIntArray (coord_system_handle, 3,
+ coord_handles, "COORDINATES") >= 0;
+ if (have_coords)
+ {
+ /* get the coordinate functions and coordinate physical minimum */
+ for (i = 1; i <= 3; i++)
+ {
+ vindices[i] = -1;
+ Util_TableGetInt (coord_handles[i-1], &vindices[i], "GAINDEX");
+ have_coords &= vindices[i] >= 0;
+ }
+ }
+ num_requested_hslabs = have_coords ? 4 : 1;
+
+ /* What processor are we on? */
+ myproc = CCTK_MyProc (GH);
+
+ /* Open the file on processor 0 */
+ filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL;
+
+ /* get the total number of grid points to check for zero-sized variables */
+ /* set the extent vector (copy from 'int' to 'CCTK_INT') */
+ CCTK_GroupgshVI (GH, 3, extent_int, vindex);
+ for (i = 0, total_hsize = 1; i < 3; i++)
+ {
+ total_hsize *= extent_int[i];
+ extent[i] = extent_int[i];
+ }
+ if (total_hsize <= 0)
+ {
+ free (fullname);
+ return (0);
+ }
+
+ downsample[0] = out_downsample_x;
+ downsample[1] = out_downsample_y;
+ downsample[2] = out_downsample_z;
+
+ /* get the hyperslab mapping */
+ mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 3,
+ direction, origin, extent,
+ downsample,
+ -1, /* table handle */
+ NULL /* conversion fn */,
+ hsize);
+ if (mapping < 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "IOSDF_Write3D: Failed to define hyperslab mapping for "
+ "variable '%s'", fullname);
+ free (fullname);
+ return (-1);
+ }
+ for (i = 0, total_hsize = 1; i < 3; i++)
+ {
+ extent_int[i] = hsize[i];
+ total_hsize *= hsize[i];
+ }
+ if (total_hsize <= 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "IOSDF_Write3D: selected hyperslab has zero size for "
+ "variable '%s'", fullname);
+ Hyperslab_FreeMapping (mapping);
+ free (fullname);
+ return (-1);
+ }
+
+ hdata[0] = hdata[1] = hdata[2] = hdata[3];
+ if (myproc == 0)
+ {
+ /* allocate hyperslab buffers */
+ hdata[0] = malloc (total_hsize * sizeof (double));
+ hdata[1] = have_coords ? malloc (3 * total_hsize * sizeof (double)) : NULL;
+ hdata[2] = hdata[1] + 1*total_hsize;
+ hdata[3] = hdata[1] + 2*total_hsize;
+ }
+
+ /* get the hyperslabs */
+ num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs,
+ NULL, vindices, NULL, NULL,
+ (void **) hdata, NULL);
+
+ /* And dump the data to file */
+ if (myproc == 0)
+ {
+ if (num_returned_hslabs == num_requested_hslabs)
+ {
+ if (have_coords)
+ {
+ /* get the staggering offset for the coordinates */
+ for (i = 0; i < 3; i++)
+ {
+ offset[i] = 0.5 * GH->cctk_delta_space[i] *
+ CCTK_StaggerDirIndex (i, gdata.stagtype);
+ }
+ for (i = 0; i < total_hsize; i++)
+ {
+ hdata[1][i] += offset[0];
+ hdata[2][i] += offset[1];
+ hdata[3][i] += offset[2];
+ }
+ }
+
+ i = have_coords ? gft_out_full (filename, dtime, extent_int, "x|y|z", 3,
+ hdata[1], hdata[0]) :
+ gft_out_brief (filename, dtime, extent_int, 3,hdata[0]);
+ if (i <= 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Error writing to 3D IOSDF output file '%s'", filename);
+ }
+ }
+ else
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "IOSDF_Write3D: Failed to extract hyperslab for "
+ "variable '%s'", fullname);
+ }
+
+ /* clean up */
+ free (hdata[0]);
+ free (hdata[1]);
+ } /* end of outputting the data by processor 0 */
+
+ free (fullname);
+
+ return (0);
+}
+
+
+/********************************************************************
+ ******************** Internal Routines ************************
+ ********************************************************************/
+/*@@
+ @routine OpenFile
+ @date Sat 12 June 2004
+ @author Thomas Radke
+ @desc
+ Opens an SDF file for a given alias name.
+ If this is the first time through, it will advertise it
+ to IOUtil.
+ @enddesc
+
+ @returntype char *
+ @returndesc
+ the full filename of the SDF output file
+ @endreturndesc
+ @@*/
+static char *OpenFile (const cGH *GH, const char *fullname, const char *alias)
+{
+ char *filename;
+ ioSDFGH *myGH;
+ ioAdvertisedFileDesc advertised_file;
+ DECLARE_CCTK_PARAMETERS
+
+
+ /* get handle for IOSDF GH extensions */
+ myGH = CCTK_GHExtension (GH, "IOSDF");
+
+ /* see if we are the first time through */
+ filename = GetNamedData (myGH->fileList_3D, alias);
+ if (! filename)
+ {
+ filename = malloc (strlen (myGH->out3D_dir) + strlen (alias) + 9);
+
+ /* open/create the file */
+ sprintf (filename, "%s%s_3D.sdf", myGH->out3D_dir, alias);
+
+ /* advertise the file for downloading and write file info */
+ advertised_file.slice = "";
+ advertised_file.thorn = CCTK_THORNSTRING;
+ advertised_file.varname = fullname;
+ advertised_file.description = "Full-dimensional variable contents";
+ advertised_file.mimetype = "application/sdf";
+
+ IOUtil_AdvertiseFile (GH, filename, &advertised_file);
+
+ /* store filename in database */
+ StoreNamedData (&myGH->fileList_3D, alias, filename);
+ }
+
+ return (filename);
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 21:53:36 +0000