grid_reader.c 15.53 KiB
#include <stdio.h>
#include <errno.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <stddef.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <unistd.h>
#include <ctype.h>
#include <math.h>
#include "grid_reader.h"
typedef struct
{
size_t size;
size_t nx, ny, nz;
double* data;
double* origin;
double** basis;
} grid;
void remove_new_line(char* str)
{
while(*str != '\n')
++str;
*str = '\0';
}
void remove_white_space(char* str)
{
size_t len = strlen(str);
char* ptr = str;
//int count = 0;
while(*ptr != '\0')
{
if(*ptr == ' ' || *ptr == '\t')
ptr++;
else
break;
}
memmove(str, ptr, sizeof(char)*(len-(size_t)(ptr-str)+1));
len = strlen(str); //new length
ptr = str+len-1;
while(*ptr == ' ')
ptr--;
//truncated the trailing white space
*(ptr+1) = '\0';
}
void* gridformat_reader(const char* name)
{
printf("Reading from file %s\n", name);
char buffer[512];
char delim[]=" \t";
int set_count = 0;
FILE* infile = fopen(name, "r");
grid* grids = (grid*)malloc(sizeof(grid));
grids->nx = 0;
grids->ny = 0;
grids->nz = 0;
grids->size = 0;
grids->origin = (double*)malloc(sizeof(double)*3);
grids->basis = (double**)malloc(sizeof(double*)*3);
for(size_t i = 0; i < 3; ++i)
grids->basis[i] = (double*)malloc(sizeof(double)*3);
grids->data = NULL;
if(infile == NULL)
{
perror("Error: ");
printf("Failure opening file %s\n", name);
return NULL;
}
size_t count = 0;
size_t count_data = 0;
while(fgets(buffer, 512, infile) != NULL)
{
char tmp[512];
char str[512];
strcpy(tmp, buffer);
remove_new_line(tmp);
remove_white_space(tmp); //strip the possible newline and puts '\0'
//skip comment
if(buffer[0] == '#')
continue;
strcpy(str, tmp);
if(strlen(str) == 0)
continue;
char* rest = strtok(tmp, " ");
if(rest == NULL)
continue;
if(!strcmp(rest, "object") && set_count == 0)
{
strtok(str," ");
strtok(NULL, " ");
strtok(NULL, " ");
strtok(NULL, " ");
strtok(NULL, " ");
rest = strtok(NULL, " ");
grids->nx = (size_t)atoi(rest);
rest = strtok(NULL, " ");
grids->ny = (size_t)atoi(rest);
rest = strtok(NULL, " ");
grids->nz = (size_t)atoi(rest);
grids->size = grids->nx*grids->ny*grids->nz;
grids->data = (double*)malloc(sizeof(double)*(grids->size));
//grids->data = (double***)malloc(sizeof(double**)*(grids->nx));
if(grids->data == NULL)
{
perror("Error in malloc:");
//printf("Error message %s\n", strerror(errno));
return NULL;
}
/*
for(size_t i = 0; i < grids->nx; ++i)
{
grids->data[i] = (double**)malloc(sizeof(double*)*(grids->ny));
if(grids->data[i] == NULL)
{
perror("Error in mallc:");
for(size_t j = 0; j < i; ++j)
free(grids->data[j]);
free(grids->data);
return NULL;
}
}
for(size_t i = 0; i < grids->nx; ++i)
{
for(size_t j = 0; j < grids->ny; ++j)
{
grids->data[i][j] = (double*)malloc(sizeof(double)*(grids->nz));
if(grids->data[i][j] == NULL)
{
perror("Error in malloc:");
for(size_t i1= 0; i1 < i; ++i1)
for(size_t j1 = 0; j1 < grids->ny; ++j1)
free(grids->data[i1][j1]);
for(size_t j1 = 0; j1 < j; ++j1)
free(grids->data[i][j1]);
for(size_t i1 = 0; i1 < grids->nx; ++i1)
free(grids->data[i1]);
free(grids->data);
return NULL;
}
}
}*/
set_count = 1;
}
else if(!strcmp(rest, "origin"))
{
strtok(str," ");
grids->origin[0] = (double)atof(strtok(NULL, " "));
grids->origin[1] = (double)atof(strtok(NULL, " "));
grids->origin[2] = (double)atof(strtok(NULL, " "));
}
else if(!strcmp(rest, "delta"))
{
strtok(str," ");
grids->basis[count][0] = (double)atof(strtok(NULL, " "));
grids->basis[count][1] = (double)atof(strtok(NULL, " "));
grids->basis[count][2] = (double)atof(strtok(NULL, " "));
++count;
}
else if(isdigit(str[0]) || isdigit(str[1]))
{
rest = strtok(str, delim);
grids->data[count_data++] = (double)atof(rest);
rest = strtok(NULL, delim);
if(rest != NULL)
grids->data[count_data++] = (double)atof(rest);
else
continue;
rest = strtok(NULL, delim);
if(rest != NULL)
grids->data[count_data++] = (double)atof(rest);
else
continue;
}
}
if(count_data != grids->size)
{
printf("Invalid file, data and count is not consistent.\n");
grid_destroy(grids);
//free(grids);
return NULL;
} printf("Sucessfully reading file %s\n", name);
if(fclose(infile) != 0)
perror("Error in closing file:");
return grids;
}
void grid_destroy(void* src)
{
free(((grid*)src)->data);
for(size_t i = 0; i < 3; ++i)
free(((grid*)src)->basis[i]);
free(((grid*)src)->origin);
free(((grid*)src)->basis);
free(src);
}
int gridformat_writer(char* name, void* src)
{
ssize_t bytes_write;
printf("Starting writing files %s\n", name);
assert(src);
grid* g = (grid*)src;
int fildes = open(name, O_RDWR | O_CREAT | O_TRUNC, 0664);
if(fildes == -1)
{
perror("Error in opeing file for writing");
//exit(errno);
//return -1;
return errno;
}
bytes_write = write(fildes, "#OpenDx file generated from grids reader/writer program\n",56);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < 56)
{
return 9999;
}
char buff[1024];
ssize_t len;
sprintf(buff, "object 1 class gridpositions counts %lu %lu %lu\n", g->nx, g->ny, g->nz);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
sprintf(buff, "origin %lf %lf %lf\n", g->origin[0], g->origin[1], g->origin[2]);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
for(size_t i = 0; i < 3; i++)
{
sprintf(buff, "delta %lf %lf %lf\n", g->basis[i][0], g->basis[i][1], g->basis[i][2]);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
}
sprintf(buff, "object 2 class gridconnections counts %lu %lu %lu\n", g->nx, g->ny, g->nz);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
sprintf(buff, "object 3 class array type float rank 0 items %lu data follows\n", g->size);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
size_t lines = g->size/3;
for(size_t l = 0; l < lines; ++l)
{
sprintf(buff, "%.15lf %.15lf %.15lf\n", g->data[3*l], g->data[3*l+1], g->data[3*l+2]);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
}
if(g->size%3 == 1)
{
sprintf(buff, "%.15lf\n", g->data[g->size-1]);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
}
else if(g->size%3 == 2) {
sprintf(buff, "%.15lf %.15lf\n", g->data[g->size-2], g->data[g->size-1]);
len = strlen(buff);
bytes_write = write(fildes, buff, len);
if(bytes_write < -1)
{
return errno;
}
else if(bytes_write < len)
{
return 9999;
}
}
if(close(fildes) == -1)
{
int err = errno;
perror("Error in closing files");
//exit(errno);
return err;
}
return 0;
}
void set_grid_data(void* src, compute_type compute, void* kernel)
{
grid* g = (grid*)src;
for(size_t i = 0; i < g->nx; ++i)
{
for(size_t j = 0; j < g->ny; ++j)
{
for(size_t k = 0; k < g->nz; ++k)
{
size_t idx = k + g->nz * j + g->nz * g->ny * i;
g->data[idx] = (*compute)(i, j, k, g->origin, g->basis, kernel);
}
}
}
}
static void* create_new_grid(size_t nx, size_t ny, size_t nz, double* origin, double** basis)
{
void* g = malloc(sizeof(grid));
((grid*)g)->nx = nx;
((grid*)g)->ny = ny;
((grid*)g)->nz = nz;
((grid*)g)->size = nx*ny*nz;
((grid*)g)->data = (double*)malloc(sizeof(double)*(((grid*)g)->size));
((grid*)g)->origin = (double*)malloc(sizeof(double)*3);
for(size_t i = 0; i < 3; ++i)
((grid*)g)->origin[i] = origin[i];
((grid*)g)->basis = (double**)malloc(sizeof(double*)*3);
for(size_t i = 0; i < 3; ++i)
((grid*)g)->basis[i] = (double*)malloc(sizeof(double)*3);
for(size_t i = 0; i < 3; ++i)
for(size_t j = 0; j < 3; ++j)
((grid*)g)->basis[i][j] = basis[i][j];
return g;
}
static size_t Index(size_t nx, size_t ny, size_t nz, int x, int y, int z)
{
x = (x < 0) ? 0 : x;
y = (y < 0) ? 0 : y;
z = (z < 0) ? 0 : z;
x = ((size_t)x >= nx) ? (int)nx-1 : x;
y = ((size_t)y >= ny) ? (int)ny-1 : y; z = ((size_t)z >= nz) ? (int)nz-1 : z;
return z + nz * (y + ny * x);
}
/*
static void inverse_3by3(double** m, double inv[3][3])
{
double det = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])-
m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])+
m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
inv[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1])/det;
inv[0][1] =-(m[0][1]*m[2][2]-m[0][2]*m[2][1])/det;
inv[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1])/det;
inv[1][0] =-(m[1][0]*m[2][2]-m[1][2]*m[2][0])/det;
inv[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0])/det;
inv[1][2] =-(m[0][0]*m[1][2]-m[0][2]*m[1][0])/det;
inv[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0])/det;
inv[2][1] =-(m[0][0]*m[2][1]-m[0][1]*m[2][0])/det;
inv[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0])/det;
}*/
//resample the grid data such that each dimensition is even number by using linear interpolation
//I currently implement for diagonal basis set
void* resample_grid_data(void* src)
{
grid* old_g = (grid*)src;
grid* new_g;
size_t nx, ny, nz;
double** b = (double**)malloc(sizeof(double*)*3);
for(size_t i = 0; i < 3; ++i)
{
b[i] = (double*)malloc(sizeof(double)*3);
memset(b[i],0,sizeof(double)*3);
}
nx = (old_g->nx/2)*2;
ny = (old_g->ny/2)*2;
nz = (old_g->nz/2)*2;
b[0][0] = old_g->basis[0][0]*(old_g->nx/(double)nx);
b[1][1] = old_g->basis[1][1]*(old_g->ny/(double)ny);
b[2][2] = old_g->basis[2][2]*(old_g->nz/(double)nz);
new_g = create_new_grid(nx,ny,nz,old_g->origin, b);
for(size_t i = 0; i < 3; ++i)
free(b[i]);
free(b);
for(size_t i = 0; i < new_g->size; ++i)
{
size_t ix = i / new_g->nz / new_g->ny;
size_t iy = i / new_g->nz % new_g->ny;
size_t iz = i % new_g->nz;
double posx = new_g->basis[0][0]*ix;
double posy = new_g->basis[1][1]*iy;
double posz = new_g->basis[2][2]*iz;
int x0 = floor(posx/old_g->basis[0][0]);
int y0 = floor(posy/old_g->basis[1][1]);
int z0 = floor(posz/old_g->basis[2][2]);
double posx0 = old_g->basis[0][0] * x0;
double posy0 = old_g->basis[1][1] * y0;
double posz0 = old_g->basis[2][2] * z0;
int x1 = x0 + 1; int y1 = y0 + 1;
int z1 = z0 + 1;
double posx1 = old_g->basis[0][0] * x1;
double posy1 = old_g->basis[1][1] * y1;
double posz1 = old_g->basis[2][2] * z1;
double xd = (posx-posx0)/(posx1-posx0);
double yd = (posy-posy0)/(posy1-posy0);
double zd = (posz-posz0)/(posz1-posz0);
double c0 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z0)] +
xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z0)])*(1.-yd) +
((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z0)] +
xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z0)])*yd;
double c1 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z1)] +
xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z1)])*(1.-yd) +
((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z1)] +
xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z1)])*yd;
new_g->data[i] = c0*(1.-zd) + c1*zd;
}
return new_g;
}
double* get_grid_origin(void* src)
{
return ((grid*)src)->origin;
}
double** get_grid_basis(void* src)
{
return ((grid*)src)->basis;
}
double* get_grid_data(void* src)
{
return ((grid*)src)->data;
}
size_t get_nx(void* src)
{
return ((grid*)src)->nx;
}
size_t get_ny(void* src)
{
return ((grid*)src)->ny;
}
size_t get_nz(void* src)
{
return ((grid*)src)->nz;
}
size_t get_size(void* src)
{
return ((grid*)src)->size;
}
void* copy_grid_format(void* ptr)
{
grid* src = (grid*)ptr;
grid* dest = (grid*)malloc(sizeof(grid));
dest->nx = src->nx;
dest->ny = src->ny;
dest->nz = src->nz;
dest->size = src->size;
dest->data = (double*)malloc(sizeof(double)*dest->size);
dest->origin = (double*)malloc(sizeof(double)*3);
dest->basis = (double**)malloc(sizeof(double*)*3); for(size_t i = 0; i < 3; ++i)
dest->basis[i] = (double*)malloc(sizeof(double)*3);
for(size_t i = 0; i < 3; ++i)
{
for(size_t j = 0; j < 3; ++j)
{
dest->basis[i][j] = src->basis[i][j];
}
dest->origin[i] = src->origin[i];
}
return dest;
}
void sum_grids(void* dest, int num, void** src)
{
grid** in = (grid**)src;
grid* out = (grid*)dest;
for(size_t idx = 0; idx < out->size; idx++)
{
out->data[idx] = 0.;
for(size_t i = 0; i < (size_t)num; ++i)
out->data[idx] += in[i]->data[idx];
}
}