Raw File
/*
 * 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 <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <assert.h>
#include <math.h>
#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_redu<nodex & nodex_redu>0)
	{
	nodex = nodex_redu;
	}
	if(nodey_redu<nodey & nodey_redu>0)
	{
	nodey = nodey_redu;
	}
	if(nodez_redu<nodez & nodez_redu>0)
	{
	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;n<nodex*nodey;n++)
				{
				//Coordinates						
				coordinate.x = botm_coord->data[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;n<nodex*nodey;n++)
				{
				//Coordinates						
				coordinate.x = surf_coord->data[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;i<nodes*caps;i++)
			{
			fprintf(file, "%f %f %f\n",coordinates[i].x,coordinates[i].y,coordinates[i].z);
			}
		}
		else
		{
		//fwrite(&coordinates,sizeof(coordinates_t),nodes*caps,file);
			
		for(i=0;i<nodes*caps;i++){
			
			fwrite(&coordinates[i].x,sizeof(float),1,file);
			fwrite(&coordinates[i].y,sizeof(float),1,file);			
			fwrite(&coordinates[i].z,sizeof(float),1,file);			
			}
		}
		
	int cells = ((nodex-1)*(nodey-1)*(nodez-1))*caps;
	//Write Cells
	fprintf(file, "CELLS %d %d\n",cells,cells*9);
	for(i=0;i<cells;i++)
		{
		fprintf(file, "8 %d %d %d %d %d %d %d %d\n",
					connectivity[i].c1,
					connectivity[i].c2,
					connectivity[i].c3,
					connectivity[i].c4,
					connectivity[i].c5,
					connectivity[i].c6,
					connectivity[i].c7,
					connectivity[i].c8	
					);
		}
		
	//Write Cell Types Hexahedron
	
		fprintf(file,"CELL_TYPES %d\n",cells);
		int j=0;
		for(i=0;i<cells;i++)
		{
		fprintf(file,"12 ");
		j++;
		if(j==8)				//nicer formating
			{
			fprintf(file,"\n");
			j=0;
			}
	
		}
		
		
		
	//Write Scalar Temperature
	fprintf(file,"POINT_DATA %d\n",nodes*caps);
	fprintf(file,"SCALARS Temperature_scalars float 1\n");
	fprintf(file,"LOOKUP_TABLE default\n");
	
	
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file,"%f ", temperature[i]);
		
		if((i+1)%nodex==0)
			{
			fprintf(file,"\n");
			
			}			
		}

	
	//Write Scalar Viscosity
	fprintf(file,"SCALARS Viscosity_scalars float 1\n");
	fprintf(file,"LOOKUP_TABLE default\n");
	
	
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file,"%f ", viscosity[i]);
		
		if((i+1)%nodex==0)
			{
			fprintf(file,"\n");
			
			}			
		}

	
	//Write Velocity Vectors
	fprintf(file,"Vectors Velocity_vectors float\n");
	
	
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file,"%f %f %f \n", velocity[i].x, velocity[i].y, velocity[i].z);
		}
		
		

	fclose(file);
		
}

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)
{
	FILE *file;
    char filename[100];
    
	
	snprintf(filename, (size_t)99, "%s.%d.vtk", filename_prefix, timestep);
    fprintf(stderr, "Writing %s\n", filename);
	file = fopen(filename, "w");
    int i;
	int nodes = nodex*nodey;

	//Write Header	
	fprintf(file,"# vtk DataFile Version 2.0\n");
	fprintf(file,"CitcomS Output %s Timestep:%d NX:%d NY:%d NZ:%d\n",
				filename_prefix,timestep,nodex,nodey,nodez);
	
	fprintf(file, "ASCII\n");
	fprintf(file, "DATASET UNSTRUCTURED_GRID\n");
	fprintf(file, "POINTS %d float\n",nodes*caps);
	
	//Write Coordinates
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file, "%f %f %f\n",coordinates[i].x,coordinates[i].y,coordinates[i].z);
		}

	//Write Cells
	int cells = ((nodex-1)*(nodey-1))*caps;
	fprintf(file, "CELLS %d %d\n",cells,cells*5);
	for(i=0;i<(((nodex-1)*(nodey-1))*caps);i++)
		{
		fprintf(file, "4 %d %d %d %d\n",
					connectivity[i].c1,
					connectivity[i].c2,
					connectivity[i].c3,
					connectivity[i].c4
					);
		}
	
	//Write Cell Types Pixel
	
		fprintf(file,"CELL_TYPES %d\n",cells);
		int j=0;
		for(i=0;i<cells;i++)
		{
		fprintf(file,"8 ");
		j++;
		if(j==8)				//nicer formating
			{
			fprintf(file,"\n");
			j=0;
			}
	
		}
	
		
		
	//Write Scalar Temperature
	fprintf(file,"POINT_DATA %d\n",nodes*caps);
	fprintf(file,"SCALARS Temperature_scalars float 1\n");
	fprintf(file,"LOOKUP_TABLE default\n");
	
	
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file,"%f ", heatflux[i]);
		
		if((i+1)%nodex==0)
			{
			fprintf(file,"\n");
			}			
		}

	//Write Velocity Vectors
	fprintf(file,"Vectors Velocity_vectors float\n");
	
	
	for(i=0;i<nodes*caps;i++)
		{
		fprintf(file,"%f %f %f \n", velocity[i].x, velocity[i].y, velocity[i].z);
		}

	fclose(file);
}

static float calc_mean(field_t *topo,int nodex,int nodey)
{
int i;
int nodes = (nodex*nodey);
double f;
for(i=0;i<nodes;i++)
	{
	f = f+topo->data[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");
}
back to top