/* * CitcomS.py by Eh Tan, Eun-seo Choi, and Pururav Thoutireddy. * Copyright (C) 2002-2005, California Institute of Technology. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include #include #include #include #include #include "hdf5.h" #define true 1 #define false 0 #define STRING_BUFFER 256 typedef struct vtk_pixel_t { int c1; int c2; int c3; int c4; } vtk_pixel_t; typedef struct hexahedron_t { int c1; int c2; int c3; int c4; int c5; int c6; int c7; int c8; } hexahedron_t; typedef struct coordinates_t { float x; float y; float z; } coordinates_t; typedef struct cap_t { int id; char name[8]; hid_t group; } cap_t; typedef struct field_t { const char *name; int rank; hsize_t *dims; hsize_t *maxdims; hsize_t *offset; hsize_t *count; int n; float *data; } field_t; static cap_t *open_cap(hid_t file_id, int capid); static herr_t close_cap(cap_t *cap); static field_t *open_field(cap_t *cap, const char *name); static herr_t read_field(cap_t *cap, field_t *field, int timestep); static herr_t close_field(field_t *field); static herr_t get_attribute_str(hid_t obj_id, const char *attr_name, char **data); static herr_t get_attribute_int(hid_t input, const char *name, int *val); static herr_t get_attribute_float(hid_t input, const char *name, float *val); static herr_t get_attribute(hid_t obj_id, const char *attr_name, hid_t mem_type_id, void *data); static herr_t get_attribute_mem(hid_t obj_id, const char *attr_name, hid_t mem_type_id, void *data); static herr_t get_attribute_disk(hid_t loc_id, const char *attr_name, void *attr_out); static herr_t get_attribute_info(hid_t obj_id, const char *attr_name, hsize_t *dims, H5T_class_t *type_class, size_t *type_size, hid_t *type_id); static coordinates_t rtf_to_xyz(coordinates_t coord); static coordinates_t velocity_to_cart(coordinates_t velocity, coordinates_t coord); static float calc_mean(field_t *topo,int nodey,int nodez); static void write_vtk_shell(coordinates_t coordinates[], hexahedron_t connectivity[], float temperature[], float viscosity[], coordinates_t velocity[], int nodex, int nodey, int nodez, int timestep, float radius_inner,int caps,int ascii); static void write_vtk_surface(coordinates_t coordinates[], vtk_pixel_t connectivity[], float heatflux[], coordinates_t velocity[], int timestep,int nodex,int nodey,int nodez,int caps, char* filename_prefix); void print_help(void); int main(int argc, char *argv[]) { char *datafile; hid_t h5file; hid_t input; herr_t status; int caps; int capid; cap_t *cap; int step; int n, i, j, k; int nodex, nodey, nodez; float radius_inner; float radius_outer; int nodex_redu=0; int nodey_redu=0; int nodez_redu=0; int initial=0; int timestep=-1; int current_t=0; int timesteps=0; int cell_counter=0; int cell_counter_surface=0; int *steps; char *endptr; field_t *coord; field_t *velocity; field_t *temperature; field_t *viscosity; //Bottom field_t *botm_coord; field_t *botm_hflux; field_t *botm_velo; field_t *botm_topo; //Surface field_t *surf_coord; field_t *surf_hflux; field_t *surf_velo; field_t *surf_topo; /////////////////////////////////////////////// int bottom=false; int surface=false; int topo=false; int ascii=false; /************************************************************************ * Parse command-line parameters. * ************************************************************************/ /* * HDF5 file must be specified as first argument. */ if (argc < 2) { fprintf(stderr, "Usage: run with -h to get help\n", argv[0]); return EXIT_FAILURE; } char c; char *hdf_filepath; char *output_filepath; extern char *optarg; extern int optind, optopt; int errflg=0; while ((c = getopt(argc, argv, ":p:o:i:t:x:y:z:bscah?")) != -1) { switch(c) { case 'p': hdf_filepath = optarg; break; case 'o': output_filepath = optarg; printf("Got output filepath\n"); break; case 'i': initial = atoi(optarg); printf("Initial: %d\n",initial); break; case 't': timesteps = atoi(optarg); printf("Timesteps: %d\n", timesteps); //Inclusive timesteps++; break; case 'x': nodex_redu=atoi(optarg); break; case 'y': nodey_redu=atoi(optarg); break; case 'z': nodez_redu=atoi(optarg); break; case ':': /* missing operand */ fprintf(stderr, "Option -%c requires an operand\n", optopt); errflg++; break; //////////////////// case 'b': bottom=true; printf("Create Bottom"); break; case 's': surface=true; printf("Create Surface"); break; case 'c': topo=true; printf("Create Topography"); break; case 'a': ascii=true; break; case '?': errflg++; print_help(); break; } } for ( ; optind < argc; optind++) { if (access(argv[optind], R_OK)) { printf("Geht\n"); } } printf("Opening HDF file...\n"); h5file = H5Fopen(hdf_filepath, H5F_ACC_RDONLY, H5P_DEFAULT); if (h5file < 0) { fprintf(stderr, "Could not open HDF5 file \"%s\"\n", argv[1]); return EXIT_FAILURE; } /************************************************************************ * Get mesh parameters. * ************************************************************************/ /* Read input group */ input = H5Gopen(h5file, "input"); if (input < 0) { fprintf(stderr, "Could not open /input group in \"%s\"\n", argv[1]); status = H5Fclose(h5file); return EXIT_FAILURE; } status = get_attribute_str(input, "datafile", &datafile); status = get_attribute_int(input, "nproc_surf", &caps); status = get_attribute_int(input, "nodex", &nodex); status = get_attribute_int(input, "nodey", &nodey); status = get_attribute_int(input, "nodez", &nodez); status = get_attribute_float(input,"radius_inner",&radius_inner); status = get_attribute_float(input,"radius_outer",&radius_outer); //Bound input params against hdf if(nodex_redu0) { nodex = nodex_redu; } if(nodey_redu0) { nodey = nodey_redu; } if(nodez_redu0) { nodez = nodez_redu; } printf("Nodex: %d\n",nodex); printf("Nodey: %d\n",nodey); printf("Nodez: %d\n",nodez); printf("Caps: %d\n",caps); /* Release input group */ status = H5Gclose(input); /************************************************************************ * Create fields using cap00 datasets as a template. * ************************************************************************/ cap = open_cap(h5file, 0); coord = open_field(cap, "coord"); velocity = open_field(cap, "velocity"); temperature = open_field(cap, "temperature"); viscosity = open_field(cap, "viscosity"); /*Create fields bottom and surface*/ botm_coord = open_field(cap,"botm/coord"); botm_hflux = open_field(cap,"botm/heatflux"); botm_velo = open_field(cap,"botm/velocity"); botm_topo = open_field(cap,"botm/topography"); surf_coord = open_field(cap,"surf/coord"); surf_hflux = open_field(cap,"surf/heatflux"); surf_velo = open_field(cap,"surf/velocity"); surf_topo = open_field(cap,"surf/topography"); status = close_cap(cap); /************************************************************************ * Output requested data. * ************************************************************************/ int iterations=0; /* Iterate over timesteps */ for(current_t = initial; current_t < timesteps; current_t++) { printf("\nProcessing timestep: %d\n",current_t); coordinates_t ordered_coordinates[((nodex*nodey*nodez)*caps)]; coordinates_t ordered_velocity[((nodex*nodey*nodez)*caps)*3]; float ordered_temperature[(nodex*nodey*nodez)*caps]; float ordered_viscosity[(nodex*nodey*nodez)*caps]; hexahedron_t connectivity[((nodex-1)*(nodey-1)*(nodez-1))*caps]; coordinates_t ordered_botm_coords[(nodex*nodey*caps)]; float ordered_botm_hflux[(nodex*nodey*caps)]; coordinates_t ordered_botm_velocity[(nodex*nodey*caps)]; coordinates_t ordered_surf_coords[(nodex*nodey*caps)]; float ordered_surf_hflux[(nodex*nodey*caps)]; coordinates_t ordered_surf_velocity[(nodex*nodey*caps)]; vtk_pixel_t connectivity_surface[(nodex*nodey*caps)]; //Holds single coordinate coordinates_t coordinate; //Holds single vector coordinates_t velocity_vector; /* Iterate over caps */ for(capid = 0; capid < caps; capid++) { cap = open_cap(h5file, capid); printf("Processing cap %d of %d\n",capid+1,caps); //snprintf(filename, (size_t)99, "%s.cap%02d.%d", datafile, capid, step); //fprintf(stderr, "Writing %s\n", filename); //file = fopen(filename, "w"); //fprintf(file, "%d x %d x %d\n", nodex, nodey, nodez); /* Read data from HDF5 file. */ read_field(cap, coord, 0); read_field(cap, velocity, current_t); read_field(cap, temperature, current_t); read_field(cap, viscosity, current_t); /*Counts iterations*/ n = 0; //Number of nodes per cap int nodes=nodex*nodey*nodez; /* Traverse data in Citcom order */ for(j = 0; j < nodey; j++) { for(i = 0; i < nodex; i++) { for(k = 0; k < nodez; k++) { //Coordinates coordinate.x = coord->data[3*n+0]; coordinate.y = coord->data[3*n+1]; coordinate.z = coord->data[3*n+2]; coordinate = rtf_to_xyz(coordinate); ordered_coordinates[n+(capid*nodes)].x = coordinate.x; ordered_coordinates[n+(capid*nodes)].y = coordinate.y; ordered_coordinates[n+(capid*nodes)].z = coordinate.z; //Velocity velocity_vector.x = velocity->data[3*n+0]; velocity_vector.y = velocity->data[3*n+1]; velocity_vector.z = velocity->data[3*n+2]; velocity_vector = velocity_to_cart(velocity_vector,coordinate); ordered_velocity[n+(capid*nodes)].x = velocity_vector.x; ordered_velocity[n+(capid*nodes)].y = velocity_vector.y; ordered_velocity[n+(capid*nodes)].z = velocity_vector.z; //Temperature ordered_temperature[n+(capid*nodes)] = temperature->data[n]; //Viscosity ordered_viscosity[n+(capid*nodes)] = viscosity->data[n]; n++; } } } //Create connectivity if(iterations==0) { //For 3d Data int i=1; //Counts X Direction int j=1; //Counts Y Direction int k=1; //Counts Z Direction for(n=0; n<((nodex*nodey*nodez)-(nodex*nodez));n++) { if ((i%nodez)==0) //X-Values { j++; //Count Y-Values } if ((j%nodex)==0) { k++; //Count Z-Values } if (((i%nodez) != 0) && ((j%nodex) != 0)) //Check if Box can be created { //Create Connectivity connectivity[cell_counter].c1 = n+(capid*(nodes)); connectivity[cell_counter].c2 = connectivity[cell_counter].c1+nodez; connectivity[cell_counter].c3 = connectivity[cell_counter].c2+nodez*nodex; connectivity[cell_counter].c4 = connectivity[cell_counter].c1+nodez*nodex; connectivity[cell_counter].c5 = connectivity[cell_counter].c1+1; connectivity[cell_counter].c6 = connectivity[cell_counter].c5+nodez; connectivity[cell_counter].c7 = connectivity[cell_counter].c6+nodez*nodex; connectivity[cell_counter].c8 = connectivity[cell_counter].c5+nodez*nodex; cell_counter++; } i++; } } //Bottom and Surface if(bottom==true){ /*Read Bottom data from HDF5 file.*/ read_field(cap,botm_coord,0); read_field(cap,botm_hflux,current_t); read_field(cap,botm_velo,current_t); read_field(cap,botm_topo,current_t); float botm_mean=0.0; if(topo=true) { botm_mean = calc_mean(botm_topo,nodex,nodey); } for(n=0;ndata[2*n+0]; coordinate.y = botm_coord->data[2*n+1]; if(topo==true) { coordinate.z = radius_inner+(botm_topo->data[n]-botm_mean)* (pow(10.0,21.0)/(pow(6371000,2)/pow(10,-6))/3300*10)/1000; //printf("Z: %f\n",coordinate.z); } else { coordinate.z = radius_inner; } coordinate = rtf_to_xyz(coordinate); ordered_botm_coords[n+(capid*nodex*nodey)].x = coordinate.x; ordered_botm_coords[n+(capid*nodex*nodey)].y = coordinate.y; ordered_botm_coords[n+(capid*nodex*nodey)].z = coordinate.z; ordered_botm_hflux[n+((capid)*nodex*nodey)] = botm_hflux->data[n]; velocity_vector.x = botm_velo->data[3*n+0]; velocity_vector.y = botm_velo->data[3*n+1]; velocity_vector.z = botm_velo->data[3*n+2]; velocity_vector = velocity_to_cart(velocity_vector,coordinate); ordered_botm_velocity[n+(capid*nodex*nodey)].x = velocity_vector.x; ordered_botm_velocity[n+(capid*nodex*nodey)].y = velocity_vector.y; ordered_botm_velocity[n+(capid*nodex*nodey)].z = velocity_vector.z; } } if(surface==true) { /*Read Surface data from HDF5 file.*/ read_field(cap,surf_coord,0); read_field(cap,surf_hflux,current_t); read_field(cap,surf_velo,current_t); read_field(cap,surf_topo,current_t); float surf_mean=0.0; if(topo=true) { surf_mean = calc_mean(surf_topo,nodex,nodey); } for(n=0;ndata[2*n+0]; coordinate.y = surf_coord->data[2*n+1]; if(topo==true) { coordinate.z = radius_outer+(surf_topo->data[n]-surf_mean)* (pow(10.0,21.0)/(pow(6371000,2)/pow(10,-6))/3300*10)/1000; //printf("Z: %f\n",coordinate.z); } else { coordinate.z = radius_outer; } coordinate = rtf_to_xyz(coordinate); ordered_surf_coords[n+(capid*nodex*nodey)].x = coordinate.x; ordered_surf_coords[n+(capid*nodex*nodey)].y = coordinate.y; ordered_surf_coords[n+(capid*nodex*nodey)].z = coordinate.z; ordered_surf_hflux[n+((capid)*nodex*nodey)] = botm_hflux->data[n]; velocity_vector.x = botm_velo->data[3*n+0]; velocity_vector.y = botm_velo->data[3*n+1]; velocity_vector.z = botm_velo->data[3*n+2]; velocity_vector = velocity_to_cart(velocity_vector,coordinate); ordered_surf_velocity[n+(capid*nodex*nodey)].x = velocity_vector.x; ordered_surf_velocity[n+(capid*nodex*nodey)].y = velocity_vector.y; ordered_surf_velocity[n+(capid*nodex*nodey)].z = velocity_vector.z; } } //Create connectivity information 2d if(iterations==0){ if(surface==true | bottom==true) { for(n=0;n<(nodex*nodey)-nodey;n++) { if ((n+1)%nodey!=0){ connectivity_surface[cell_counter_surface].c1 = n+(capid*((nodex)*(nodey))); connectivity_surface[cell_counter_surface].c2 = connectivity_surface[cell_counter_surface].c1+1; connectivity_surface[cell_counter_surface].c3 = connectivity_surface[cell_counter_surface].c1+nodey; connectivity_surface[cell_counter_surface].c4 = connectivity_surface[cell_counter_surface].c3+1; cell_counter_surface++; } } } } close_cap(cap); } iterations++; //Write data to file write_vtk_shell(ordered_coordinates, connectivity, ordered_temperature, ordered_viscosity, ordered_velocity, nodex, nodey, nodez, current_t, radius_inner,caps,ascii); if(bottom==true){ write_vtk_surface(ordered_botm_coords,connectivity_surface,ordered_botm_hflux, ordered_botm_velocity,current_t,nodex,nodey,nodez,caps,"bottom"); } if(surface==true){ write_vtk_surface(ordered_surf_coords,connectivity_surface,ordered_surf_hflux, ordered_surf_velocity,current_t,nodex,nodey,nodez,caps,"surface"); } }//end timesteps loop /* Release resources. */ status = close_field(coord); status = close_field(velocity); status = close_field(temperature); status = close_field(viscosity); status = H5Fclose(h5file); return EXIT_SUCCESS; } static cap_t *open_cap(hid_t file_id, int capid) { cap_t *cap; cap = (cap_t *)malloc(sizeof(cap_t)); cap->id = capid; snprintf(cap->name, (size_t)7, "cap%02d", capid); cap->group = H5Gopen(file_id, cap->name); if (cap->group < 0) { free(cap); return NULL; } return cap; } static herr_t close_cap(cap_t *cap) { herr_t status; if (cap != NULL) { cap->id = -1; cap->name[0] = '\0'; status = H5Gclose(cap->group); free(cap); } return 0; } static field_t *open_field(cap_t *cap, const char *name) { hid_t dataset; hid_t dataspace; herr_t status; int d; int rank; field_t *field; if (cap == NULL) return NULL; /* Allocate field and initialize. */ field = (field_t *)malloc(sizeof(field_t)); field->name = name; field->rank = 0; field->dims = NULL; field->maxdims = NULL; field->n = 0; dataset = H5Dopen(cap->group, name); if(dataset < 0) { free(field); return NULL; } dataspace = H5Dget_space(dataset); if (dataspace < 0) { free(field); return NULL; } /* Calculate shape of field. */ rank = H5Sget_simple_extent_ndims(dataspace); field->rank = rank; field->dims = (hsize_t *)malloc(rank * sizeof(hsize_t)); field->maxdims = (hsize_t *)malloc(rank * sizeof(hsize_t)); status = H5Sget_simple_extent_dims(dataspace, field->dims, field->maxdims); /* DEBUG printf("Field %s shape (", name); for(d = 0; d < rank; d++) printf("%d,", (int)(field->dims[d])); printf(")\n"); // */ /* Allocate memory for hyperslab selection parameters. */ field->offset = (hsize_t *)malloc(rank * sizeof(hsize_t)); field->count = (hsize_t *)malloc(rank * sizeof(hsize_t)); /* Allocate enough memory for a single time-slice buffer. */ field->n = 1; if (field->maxdims[0] == H5S_UNLIMITED) for(d = 1; d < rank; d++) field->n *= field->dims[d]; else for(d = 0; d < rank; d++) field->n *= field->dims[d]; field->data = (float *)malloc(field->n * sizeof(float)); /* Release resources. */ status = H5Sclose(dataspace); status = H5Dclose(dataset); return field; } static herr_t read_field(cap_t *cap, field_t *field, int timestep) { hid_t dataset; hid_t filespace; hid_t memspace; herr_t status; int d; if (cap == NULL || field == NULL) return -1; dataset = H5Dopen(cap->group, field->name); if (dataset < 0) return -1; for(d = 0; d < field->rank; d++) { field->offset[d] = 0; field->count[d] = field->dims[d]; } if (field->maxdims[0] == H5S_UNLIMITED) { field->offset[0] = timestep; field->count[0] = 1; } /* DEBUG printf("Reading step %d on field %s with offset (", timestep, field->name); for(d = 0; d < field->rank; d++) printf("%d,", (int)(field->offset[d])); printf(") and count ("); for(d = 0; d < field->rank; d++) printf("%d,", (int)(field->count[d])); printf(")\n"); // */ filespace = H5Dget_space(dataset); status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, field->offset, NULL, field->count, NULL); memspace = H5Screate_simple(field->rank, field->count, NULL); status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, filespace, H5P_DEFAULT, field->data); status = H5Sclose(filespace); status = H5Sclose(memspace); status = H5Dclose(dataset); return 0; } static herr_t close_field(field_t *field) { if (field != NULL) { free(field->dims); free(field->maxdims); free(field->offset); free(field->count); free(field->data); free(field); } return 0; } static herr_t get_attribute_str(hid_t obj_id, const char *attr_name, char **data) { hid_t attr_id; hid_t attr_type; size_t type_size; herr_t status; *data = NULL; attr_id = H5Aopen_name(obj_id, attr_name); if (attr_id < 0) return -1; attr_type = H5Aget_type(attr_id); if (attr_type < 0) goto out; /* Get the size */ type_size = H5Tget_size(attr_type); if (type_size < 0) goto out; /* malloc enough space for the string, plus 1 for trailing '\0' */ *data = (char *)malloc(type_size + 1); status = H5Aread(attr_id, attr_type, *data); if (status < 0) goto out; /* Set the last character to '\0' in case we are dealing with * null padded or space padded strings */ (*data)[type_size] = '\0'; status = H5Tclose(attr_type); if (status < 0) goto out; status = H5Aclose(attr_id); if (status < 0) return -1; return 0; out: H5Tclose(attr_type); H5Aclose(attr_id); if (*data) free(*data); return -1; } static herr_t get_attribute_int(hid_t input, const char *name, int *val) { hid_t attr_id; hid_t type_id; H5T_class_t type_class; size_t type_size; herr_t status; char *strval; attr_id = H5Aopen_name(input, name); type_id = H5Aget_type(attr_id); type_class = H5Tget_class(type_id); type_size = H5Tget_size(type_id); H5Tclose(type_id); H5Aclose(attr_id); switch(type_class) { case H5T_STRING: status = get_attribute_str(input, name, &strval); if (status < 0) return -1; *val = atoi(strval); free(strval); return 0; case H5T_INTEGER: status = get_attribute(input, name, H5T_NATIVE_INT, val); if (status < 0) return -1; return 0; } return -1; } static herr_t get_attribute_float(hid_t input, const char *name, float *val) { hid_t attr_id; hid_t type_id; H5T_class_t type_class; size_t type_size; herr_t status; char *strval; attr_id = H5Aopen_name(input, name); type_id = H5Aget_type(attr_id); type_class = H5Tget_class(type_id); type_size = H5Tget_size(type_id); H5Tclose(type_id); H5Aclose(attr_id); switch(type_class) { case H5T_STRING: status = get_attribute_str(input, name, &strval); if (status < 0) return -1; *val = atof(strval); free(strval); return 0; case H5T_FLOAT: status = get_attribute(input, name, H5T_NATIVE_FLOAT, val); if (status < 0) return -1; return 0; } return -1; } static herr_t get_attribute(hid_t obj_id, const char *attr_name, hid_t mem_type_id, void *data) { herr_t status; status = get_attribute_mem(obj_id, attr_name, mem_type_id, data); if (status < 0) return -1; return 0; } static herr_t get_attribute_mem(hid_t obj_id, const char *attr_name, hid_t mem_type_id, void *data) { hid_t attr_id; herr_t status; attr_id = H5Aopen_name(obj_id, attr_name); if (attr_id < 0) return -1; status = H5Aread(attr_id, mem_type_id, data); if (status < 0) { H5Aclose(attr_id); return -1; } status = H5Aclose(attr_id); if (status < 0) return -1; return 0; } static herr_t get_attribute_disk(hid_t loc_id, const char *attr_name, void *attr_out) { hid_t attr_id; hid_t attr_type; herr_t status; attr_id = H5Aopen_name(loc_id, attr_name); if (attr_id < 0) return -1; attr_type = H5Aget_type(attr_id); if (attr_type < 0) goto out; status = H5Aread(attr_id, attr_type, attr_out); if (status < 0) goto out; status = H5Tclose(attr_type); if (status < 0) goto out; status = H5Aclose(attr_id); if (status < 0) return -1; return 0; out: H5Tclose(attr_type); H5Aclose(attr_id); return -1; } static herr_t get_attribute_info(hid_t obj_id, const char *attr_name, hsize_t *dims, H5T_class_t *type_class, size_t *type_size, hid_t *type_id) { hid_t attr_id; hid_t space_id; herr_t status; int rank; /* Open the attribute. */ attr_id = H5Aopen_name(obj_id, attr_name); if (attr_id < 0) return -1; /* Get an identifier for the datatype. */ *type_id = H5Aget_type(attr_id); /* Get the class. */ *type_class = H5Tget_class(*type_id); /* Get the size. */ *type_size = H5Tget_size(*type_id); /* Get the dataspace handle */ space_id = H5Aget_space(attr_id); if (space_id < 0) goto out; /* Get dimensions */ rank = H5Sget_simple_extent_dims(space_id, dims, NULL); if (rank < 0) goto out; /* Terminate access to the dataspace */ status = H5Sclose(space_id); if (status < 0) goto out; /* End access to the attribute */ status = H5Aclose(attr_id); if (status < 0) goto out; return 0; out: H5Tclose(*type_id); H5Aclose(attr_id); return -1; } static coordinates_t rtf_to_xyz(coordinates_t coord) { coordinates_t output; output.x = coord.z * sin(coord.x) * cos(coord.y); output.y = coord.z * sin(coord.x) * sin(coord.y); output.z = coord.z * cos(coord.x); return output; } static coordinates_t velocity_to_cart(coordinates_t velocity, coordinates_t coord) { coordinates_t output; output.x = velocity.z*sin(coord.x)*cos(coord.y)+velocity.x*cos(coord.x)*cos(coord.y)-velocity.y*sin(coord.y); output.y = velocity.z*sin(coord.x)*sin(coord.y)+velocity.x*cos(coord.x)*sin(coord.y)+velocity.y*cos(coord.y); output.z = velocity.z*cos(coord.x)-velocity.x*sin(coord.x); return output; } static void write_vtk_shell(coordinates_t coordinates[], hexahedron_t connectivity[], float temperature[], float viscosity[], coordinates_t velocity[], int nodex, int nodey, int nodez, int timestep, float radius_inner,int caps,int ascii) { FILE *file; char filename[100]; char header[STRING_BUFFER]; char puffer; int i; snprintf(filename, (size_t)99, "%s.%d.vtk", "datafile", timestep); fprintf(stderr, "Writing %s\n", filename); file = fopen(filename, "w"); int nodes = nodex*nodey*nodez; char *type_info; if(ascii==true) { type_info = "Ascii"; } else { type_info = "Binary"; } //Write Header sprintf(header,"# vtk DataFile Version 2.0\n\ CitcomS Output Timestep:%d NX:%d NY:%d NZ:%d Radius_Inner:%f\n\ %s\n\ DATASET UNSTRUCTURED_GRID\n",timestep,nodex,nodey,nodez,radius_inner,ascii ? "ASCII" : "BINARY"); fprintf(file,header); fprintf(file, "POINTS %d float\n",nodes*caps); //printf("Float: %d\n",sizeof(float)); //printf("Field: %d\n",sizeof(coordinates_t)); //Write Coordinates if(ascii==true) { for(i=0;idata[i]; } return (float) f/nodes; } void print_help(void){ printf("Converts CitcomS HDF5 file to VTK\n"); printf("-p, --path [path to hdf] \n\t Specify input file.\n"); printf("-o, --output [output filename] \n\t Specify the path to the folder for output files.\n"); printf("-i, --initial [initial timestep] \n\t Specify initial timestep to export. If not \n \ \t specified script starts exporting from timestep 0.\n"); printf("-t, --timestep [max timestep] \n\t Specify to which timestep you want to export. If not\n \ \t specified export all all timestep starting from intial timestep.\n"); printf("-x, --nx_reduce [nx] \n\t Set new nx to reduce output grid.\n"); printf("-y, --ny_reduce [ny] \n\t Set new ny to reduce output grid.\n"); printf("-z, --nz_reduce [nz] \n\t Set new nz to reduce output grid.\n"); printf("-b, --bottom \n\t Set to export Bottom information to Vtk file.\n"); printf("-s, --surface \n\t Set to export Surface information to Vtk file.\n"); printf("-c, --createtopo \n\t Set to create topography information in bottom and surface Vtk file.\n"); printf("-a, --ascii \n\t Create Vtk ASCII encoded files instead if binary.\n"); printf("-h, --help, -? \n\t Print this help.\n"); }