diff --git a/examples/Makefile.am b/examples/Makefile.am
index ce5007d37733daefea82bf46a4754954a5389921..361de521c051eea0613e524505648f68ef03021c 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -21,6 +21,8 @@ AUTOMAKE_OPTIONS=gnu
# Add the source directory and debug to CFLAGS
AM_CFLAGS = -g -Wall -Werror -I../src $(OPENMP_CFLAGS) -DCPU_TPS=2.67e9
+AM_LDFLAGS = $(HDF5_LDFLAGS)
+
# Set-up the library
bin_PROGRAMS = test
diff --git a/examples/test.c b/examples/test.c
index 0498a9d024108cceac4d4e637750a7803cdbf0fa..5bcec2a09b20d51a9d3497f98c8592a67d03f5c5 100644
--- a/examples/test.c
+++ b/examples/test.c
@@ -628,108 +628,176 @@ int main ( int argc , char *argv[] ) {
int nr_threads = 1, nr_queues = -1, runs = 1;
int data[2];
double dim[3] = { 1.0 , 1.0 , 1.0 }, shift[3] = { 0.0 , 0.0 , 0.0 };
- double r_min = 0.01, r_max = 0.1, h_max = -1.0 , scaling = 1.0, rho = 0.0;
+ double /*r_min = 0.01, r_max = 0.1,*/ h_max = -1.0 , scaling = 1.0, rho = 0.0;
struct part *parts = NULL, *p;
struct space s;
struct engine e;
+ char ICfileName[200];
ticks tic;
/* Init the space. */
bzero( &s , sizeof(struct space) );
-
- /* Parse the options. */
- while ( ( c = getopt( argc , argv , "a:b:p:d:N:c:h:v:m:s:t:q:r:i:m:z:" ) ) != -1 )
- switch ( c ) {
- case 'N':
- if ( sscanf( optarg , "%d" , &N ) != 1 )
- error( "Error parsing number of particles." );
- if ( posix_memalign( (void *)&parts , 32 , N * sizeof(struct part) ) != 0 )
- error( "Call to posix_memalign failed." );
- for ( k = 0 ; k < N ; k++ ) {
- parts[k].x[0] = ((double)rand()) / RAND_MAX * dim[0];
- parts[k].x[1] = ((double)rand()) / RAND_MAX * dim[1];
- parts[k].x[2] = ((double)rand()) / RAND_MAX * dim[2];
- parts[k].id = k;
- parts[k].h = r_min + ((r_max - r_min)*rand())/RAND_MAX;
- parts[k].mass = 1.0f;
- parts[k].u = 1.0f;
- }
- printf( "main: allocated memory for %i parts.\n" , N ); fflush(stdout);
- break;
- case 'a':
- if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
- error( "Error parsing cutoff scaling." );
- printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout);
- for ( k = 0 ; k < N ; k++ )
- parts[k].h *= scaling;
- break;
- case 'b':
- if ( sscanf( optarg , "%lf %lf %lf" , &dim[0] , &dim[1] , &dim[2] ) != 3 )
- error( "Error parsing box dimensions." );
- break;
- case 'c':
- printf( "main: reading parts from %s...\n" , optarg ); fflush(stdout);
- if ( parts == NULL && posix_memalign( (void *)&parts , 16 , N * sizeof(struct part) ) != 0 )
- error( "Call to calloc failed." );
- read_coords( optarg , parts , N );
- break;
- case 'd':
- printf( "main: reading dt from %s...\n" , optarg ); fflush(stdout);
- read_dt( optarg , parts , N );
- break;
- case 'h':
- printf( "main: reading cutoffs from %s...\n" , optarg ); fflush(stdout);
- read_cutoffs( optarg , parts , N );
- break;
- case 'i':
- printf( "main: reading ids from %s...\n" , optarg ); fflush(stdout);
- read_id( optarg , parts , N );
- break;
- case 'm':
- if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
- error( "Error parsing h_max." );
- printf( "main: maximum h set to %e.\n" , h_max );
- break;
- case 'p':
- if ( sscanf( optarg , "%d" , &periodic ) != 1 )
- error( "Error parsing periodicity." );
- printf( "main: periodicity switched %s.\n" , periodic ? "on" : "off" );
- break;
- case 'q':
- if ( sscanf( optarg , "%d" , &nr_queues ) != 1 )
- error( "Error parsing number of queues." );
- break;
- case 'r':
- if ( sscanf( optarg , "%d" , &runs ) != 1 )
- error( "Error parsing number of runs." );
- break;
- case 's':
- if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 )
- error( "Error parsing shift." );
- for ( k = 0 ; k < N ; k++ ) {
- parts[k].x[0] += shift[0];
- parts[k].x[1] += shift[1];
- parts[k].x[2] += shift[2];
- }
- printf( "main: shifted parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] );
- break;
- case 't':
- if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
- error( "Error parsing number of threads." );
- omp_set_num_threads( nr_threads );
- break;
- case 'z':
- if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
- error( "Error parsing split size." );
- printf( "main: split size set to %i.\n" , space_splitsize );
- break;
- case '?':
- error( "Unknown option." );
- break;
- }
+
+ /* Parse the options */
+ while ( ( c = getopt( argc , argv , "f:a:m:s:t:q:r:z:" ) ) != -1 )
+ switch( c )
+ {
+ case 'f':
+ if( !strcpy(ICfileName, optarg))
+ error("Error parsing IC file name.");
+ printf("main: IC to be read from file '%s'\n", ICfileName);
+ break;
+ case 'a':
+ if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
+ error( "Error parsing cutoff scaling." );
+ printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout);
+ break;
+ case 'm':
+ if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
+ error( "Error parsing h_max." );
+ printf( "main: maximum h set to %e.\n" , h_max ); fflush(stdout);
+ break;
+ case 'q':
+ if ( sscanf( optarg , "%d" , &nr_queues ) != 1 )
+ error( "Error parsing number of queues." );
+ break;
+ case 'r':
+ if ( sscanf( optarg , "%d" , &runs ) != 1 )
+ error( "Error parsing number of runs." );
+ break;
+ case 's':
+ if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 )
+ error( "Error parsing shift." );
+ printf( "main: will shift parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] );
+ break;
+ case 't':
+ if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
+ error( "Error parsing number of threads." );
+ omp_set_num_threads( nr_threads );
+ break;
+ case 'z':
+ if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
+ error( "Error parsing split size." );
+ printf( "main: split size set to %i.\n" , space_splitsize );
+ break;
+ case '?':
+ error( "Unknown option." );
+ break;
+
+ }
+
+ /* while ( ( c = getopt( argc , argv , "a:b:p:d:N:c:h:v:m:s:t:q:r:i:m:z:" ) ) != -1 ) */
+ /* switch ( c ) { */
+ /* case 'N': */
+ /* if ( sscanf( optarg , "%d" , &N ) != 1 ) */
+ /* error( "Error parsing number of particles." ); */
+ /* if ( posix_memalign( (void *)&parts , 32 , N * sizeof(struct part) ) != 0 ) */
+ /* error( "Call to posix_memalign failed." ); */
+ /* for ( k = 0 ; k < N ; k++ ) { */
+ /* parts[k].x[0] = ((double)rand()) / RAND_MAX * dim[0]; */
+ /* parts[k].x[1] = ((double)rand()) / RAND_MAX * dim[1]; */
+ /* parts[k].x[2] = ((double)rand()) / RAND_MAX * dim[2]; */
+ /* parts[k].id = k; */
+ /* parts[k].h = r_min + ((r_max - r_min)*rand())/RAND_MAX; */
+ /* parts[k].mass = 1.0f; */
+ /* parts[k].u = 1.0f; */
+ /* } */
+ /* printf( "main: allocated memory for %i parts.\n" , N ); fflush(stdout); */
+ /* break; */
+ /* case 'a': */
+ /* if ( sscanf( optarg , "%lf" , &scaling ) != 1 ) */
+ /* error( "Error parsing cutoff scaling." ); */
+ /* printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout); */
+ /* for ( k = 0 ; k < N ; k++ ) */
+ /* parts[k].h *= scaling; */
+ /* break; */
+ /* case 'b': */
+ /* if ( sscanf( optarg , "%lf %lf %lf" , &dim[0] , &dim[1] , &dim[2] ) != 3 ) */
+ /* error( "Error parsing box dimensions." ); */
+ /* break; */
+ /* case 'c': */
+ /* printf( "main: reading parts from %s...\n" , optarg ); fflush(stdout); */
+ /* if ( parts == NULL && posix_memalign( (void *)&parts , 16 , N * sizeof(struct part) ) != 0 ) */
+ /* error( "Call to calloc failed." ); */
+ /* read_coords( optarg , parts , N ); */
+ /* break; */
+ /* case 'd': */
+ /* printf( "main: reading dt from %s...\n" , optarg ); fflush(stdout); */
+ /* read_dt( optarg , parts , N ); */
+ /* break; */
+ /* case 'h': */
+ /* printf( "main: reading cutoffs from %s...\n" , optarg ); fflush(stdout); */
+ /* read_cutoffs( optarg , parts , N ); */
+ /* break; */
+ /* case 'i': */
+ /* printf( "main: reading ids from %s...\n" , optarg ); fflush(stdout); */
+ /* read_id( optarg , parts , N ); */
+ /* break; */
+ /* case 'm': */
+ /* if ( sscanf( optarg , "%lf" , &h_max ) != 1 ) */
+ /* error( "Error parsing h_max." ); */
+ /* printf( "main: maximum h set to %e.\n" , h_max ); */
+ /* break; */
+ /* case 'p': */
+ /* if ( sscanf( optarg , "%d" , &periodic ) != 1 ) */
+ /* error( "Error parsing periodicity." ); */
+ /* printf( "main: periodicity switched %s.\n" , periodic ? "on" : "off" ); */
+ /* break; */
+ /* case 'q': */
+ /* if ( sscanf( optarg , "%d" , &nr_queues ) != 1 ) */
+ /* error( "Error parsing number of queues." ); */
+ /* break; */
+ /* case 'r': */
+ /* if ( sscanf( optarg , "%d" , &runs ) != 1 ) */
+ /* error( "Error parsing number of runs." ); */
+ /* break; */
+ /* case 's': */
+ /* if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 ) */
+ /* error( "Error parsing shift." ); */
+ /* for ( k = 0 ; k < N ; k++ ) { */
+ /* parts[k].x[0] += shift[0]; */
+ /* parts[k].x[1] += shift[1]; */
+ /* parts[k].x[2] += shift[2]; */
+ /* } */
+ /* printf( "main: shifted parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] ); */
+ /* break; */
+ /* case 't': */
+ /* if ( sscanf( optarg , "%d" , &nr_threads ) != 1 ) */
+ /* error( "Error parsing number of threads." ); */
+ /* omp_set_num_threads( nr_threads ); */
+ /* break; */
+ /* case 'z': */
+ /* if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 ) */
+ /* error( "Error parsing split size." ); */
+ /* printf( "main: split size set to %i.\n" , space_splitsize ); */
+ /* break; */
+ /* case '?': */
+ /* error( "Unknown option." ); */
+ /* break; */
+ /* } */
+
/* How large are the parts? */
- printf( "main: sizeof(struct part) is %li bytes.\n" , (long int)sizeof( struct part ) );
+ printf( "main: sizeof(struct part) is %li bytes.\n" , (long int)sizeof( struct part ));
+
+ /* Read particles and space information from (GADGET) IC */
+ tic = getticks();
+ read_ic(ICfileName, dim, &parts, &N, &periodic);
+ printf( "main: reading particle properties took %.3f ms.\n" , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
+
+ /* Apply h scaling */
+ if(scaling != 1.0)
+ for ( k = 0 ; k < N ; k++ )
+ parts[k].h *= scaling;
+
+ /* Apply shift */
+ if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )
+ for ( k = 0 ; k < N ; k++ ) {
+ parts[k].x[0] += shift[0];
+ parts[k].x[1] += shift[1];
+ parts[k].x[2] += shift[2];
+ }
+
/* Dump the kernel to make sure its ok. */
// kernel_dump( 100 );
diff --git a/src/Makefile.am b/src/Makefile.am
index addcde017e28d45a15efcc94ab77576b7a90755e..33907eb666c96b038a9cd7a1e89abfc7df847167 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -22,7 +22,7 @@ AUTOMAKE_OPTIONS=gnu
AM_CFLAGS = -g -O3 -Wall -Werror -ffast-math -fstrict-aliasing $(SIMD_FLAGS) $(CFLAGS) $(OPENMP_CFLAGS) -DTIMER -DCOUNTER -DCPU_TPS=2.67e9
# Assign a "safe" version number
-AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) -version-info 0:0:0
+AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) $(HDF5_LDFLAGS) -version-info 0:0:0
# Build the libgadgetsmp library
lib_LTLIBRARIES = libgadgetsmp.la
@@ -32,3 +32,10 @@ libgadgetsmp_la_SOURCES = space.c runner.c queue.c task.c cell.c engine.c ic.c
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h gadgetsmp.h ic.h
+# Set-up the library
+#bin_PROGRAMS = test
+
+# Sources for test
+#test_SOURCES = main.c
+#test_CFLAGS = -DCOUNTER -DTIMER $(AM_CFLAGS)
+#test_LDADD = ../src/.libs/libgadgetsmp.a
diff --git a/src/gadgetsmp.h b/src/gadgetsmp.h
index 3193cc8200cab4906a027faa775659e128cf3d62..67a01137a8b0f61b7bdfad57465a26e1890eaba5 100644
--- a/src/gadgetsmp.h
+++ b/src/gadgetsmp.h
@@ -32,4 +32,4 @@
#include "runner_iact.h"
#include "engine.h"
#include "runner.h"
-
+#include "ic.h"
diff --git a/src/ic.c b/src/ic.c
index 315395fa26760474c767e6f509b725b2e08bb01b..11025277e2443d0cfdf0d289df098b45f0f55b2b 100644
--- a/src/ic.c
+++ b/src/ic.c
@@ -22,31 +22,287 @@
#ifdef HAVE_HDF5
+
/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <stddef.h>
#include <hdf5.h>
+
#include "task.h"
#include "lock.h"
#include "part.h"
#include "space.h"
+/* Error macro. */
+#define error(s) { printf( "%s:%s():%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
+
+
+/**
+ * @brief The different types of data used in the GADGET IC files.
+ *
+ * (This is admittedly a poor substitute to C++ templates...)
+ */
+enum DATA_TYPE{INT, LONG, LONGLONG, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE};
+
+/**
+ * @brief Converts a C data type to the HDF5 equivalent.
+ *
+ * This function is a trivial wrapper around the HDF5 types but allows
+ * to change the exact storage types matching the code types in a transparent way.
+ */
+hid_t hdf5Type(enum DATA_TYPE type)
+{
+ switch(type)
+ {
+ case INT: return H5T_NATIVE_INT;
+ case UINT: return H5T_NATIVE_UINT;
+ case LONG: return H5T_NATIVE_LONG;
+ case ULONG: return H5T_NATIVE_ULONG;
+ case LONGLONG: return H5T_NATIVE_LLONG;
+ case ULONGLONG: return H5T_NATIVE_ULLONG;
+ case FLOAT: return H5T_NATIVE_FLOAT;
+ case DOUBLE: return H5T_NATIVE_DOUBLE;
+ default: error("Unknown type");
+ }
+}
+
+
+/**
+ * @brief Returns the memory size of the data type
+ */
+size_t sizeOfType(enum DATA_TYPE type)
+{
+ switch(type)
+ {
+ case INT: return sizeof(int);
+ case UINT: return sizeof(unsigned int);
+ case LONG: return sizeof(long);
+ case ULONG: return sizeof(unsigned long);
+ case LONGLONG: return sizeof(long long);
+ case ULONGLONG: return sizeof(unsigned long long);
+ case FLOAT: return sizeof(float);
+ case DOUBLE: return sizeof(double);
+ default: error("Unknown type");
+ }
+}
+
+
+/**
+ * @brief Reads an attribute from a given HDF5 group.
+ *
+ * @param grp The group from which to read.
+ * @param name The name of the attribute to read.
+ * @param type The #DATA_TYPE of the attribute.
+ * @param data (output) The attribute read from the HDF5 group.
+ *
+ * Calls #error() if an error occurs.
+ */
+void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data)
+{
+ hid_t h_attr=0, h_err=0;
+
+ h_attr = H5Aopen(grp, name, H5P_DEFAULT);
+ if(h_attr < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while opening attribute '%s'\n", name);
+ error(buf);
+ }
+
+ h_err = H5Aread(h_attr, hdf5Type(type), data);
+ if(h_err < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while reading attribute '%s'\n", name);
+ error(buf);
+ }
+
+ H5Aclose(h_attr);
+}
+
+
+/**
+ * @brief Reads a data array from a given HDF5 group.
+ *
+ * @param grp The group from which to read.
+ * @param name The name of the array to read.
+ * @param type The #DATA_TYPE of the attribute.
+ * @param N The number of particles.
+ * @param dim The dimension of the data (1 for scalar, 3 for vector)
+ * @param part_c A (char*) pointer on the first occurence of the field of interest in the parts array
+ *
+ * @todo A better version using HDF5 hyperslabs to read the file directly into the part array
+ * will be written once the strucutres have been stabilized.
+ *
+ * Calls #error() if an error occurs.
+ */
+void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, char* part_c)
+{
+ hid_t h_data=0, h_err=0, h_type=0;
+ void* temp;
+ int i=0;
+ size_t typeSize = sizeOfType(type);
+ size_t copySize = typeSize * dim;
+ size_t partSize = sizeof(struct part);
+ char* temp_c = 0;
+
+ printf("readArray: Reading '%s' array...\n", name);
+
+ /* Open data space */
+ h_data = H5Dopen(grp, name);
+ if(h_data < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while opening data space '%s'\n", name);
+ error(buf);
+ }
+
+ /* Check data type */
+ h_type = H5Dget_type(h_data);
+ if(h_type < 0)
+ error("Unable to retrieve data type from the file");
+ if(!H5Tequal(h_type, hdf5Type(type)))
+ error("Non-matching types between the code and the file");
+
+ /* Allocate temporary buffer */
+ temp = malloc(N * dim * sizeOfType(type));
+ if(temp == NULL)
+ error("Unable to allocate memory for temporary buffer");
+
+ /* Read HDF5 dataspace in temporary buffer */
+ /* Dirty version that happens to work for vectors but should be improved */
+ /* Using HDF5 dataspaces would be better */
+ h_err = H5Dread(h_data, hdf5Type(type), H5S_ALL, H5S_ALL, H5P_DEFAULT, temp);
+ if(h_err < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while reading data array '%s'\n", name);
+ error(buf);
+ }
+
+ /* Copy temporary buffer to particle data */
+ temp_c = temp;
+ for(i=0; i<N; ++i)
+ memcpy(part_c+i*partSize, &temp_c[i*copySize], copySize);
+
+ /* Free and close everything */
+ free(temp);
+ H5Tclose(h_type);
+ H5Dclose(h_data);
+}
+
+/**
+ * @brief A helper macro to call the readArrayBackEnd function more easily.
+ *
+ * @param grp The group from which to read.
+ * @param name The name of the array to read.
+ * @param type The #DATA_TYPE of the attribute.
+ * @param N The number of particles.
+ * @param dim The dimension of the data (1 for scalar, 3 for vector)
+ * @param part The array of particles to fill
+ * @param field The name of the field (C code name as defined in part.h) to fill
+ *
+ */
+#define readArray(grp, name, type, N, dim, part, field) readArrayBackEnd(grp, name, type, N, dim, (char*)(&(part[0]).field))
/**
- * @brief Reads an HDF5 initial condition file
+ * @brief Reads an HDF5 initial condition file (GADGET-3 type)
+ *
+ * @param fileName The file to read.
+ * @param dim (output) The dimension of the volume read from the file.
+ * @param parts (output) The array of #part read from the file.
+ * @param N (output) The number of particles read from the file.
+ * @param periodic (output) 1 if the volume is periodic, 0 if not.
*
* Opens the HDF5 file fileName and reads the particles contained
* in the parts array. N is the returned number of particles found
* in the file.
*
+ * @warning Can not read snapshot distributed over more than 1 file !!!
+ * @todo Read snaphsots distributed in more than one file.
+ *
+ * Calls #error() if an error occurs.
+ *
*/
-void read_ic ( char* fileName, struct part *parts, int* N)
+void read_ic ( char* fileName, double dim[3], struct part **parts, int* N, int* periodic)
{
+ hid_t h_file=0, h_grp=0;
+ double boxSize=0.; /* GADGET has only cubic boxes (in cosmological mode) */
+ int numParticles[6]={0}; /* GADGET has 6 particle types. We only keep the type 0*/
+
+ /* Open file */
+ printf("read_ic: Opening file '%s' as IC.\n", fileName);
+ h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if(h_file < 0)
+ {
+ char buf[200];
+ sprintf(buf, "Error while opening file '%s'", fileName);
+ error(buf);
+ }
+
+ /* Open header to read simulation properties */
+ printf("read_ic: Reading runtime parameters...\n");
+ h_grp = H5Gopen(h_file, "/RuntimePars");
+ if(h_grp < 0)
+ error("Error while opening runtime parameters\n");
+
+ /* Read the relevant information */
+ readAttribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
+
+ /* Close runtime parameters */
+ H5Gclose(h_grp);
+
+ /* Open header to read simulation properties */
+ printf("read_ic: Reading file header...\n");
+ h_grp = H5Gopen(h_file, "/Header");
+ if(h_grp < 0)
+ error("Error while opening file header\n");
+
+ /* Read the relevant information and print status */
+ readAttribute(h_grp, "BoxSize", DOUBLE, &boxSize);
+ readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
+ *N = numParticles[0];
+ dim[0] = dim[1] = dim[2] = boxSize;
+ printf("read_ic: Found %d particles in a %speriodic box of size [%f %f %f]\n",
+ *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]);
+
+ /* Close header */
+ H5Gclose(h_grp);
+
+ /* Allocate memory to store particles */
+ if(posix_memalign( (void*)parts , 32 , *N * sizeof(struct part)) != 0)
+ error("Error while allocating memory for particles");
+
+ printf("read_ic: Allocated %8.2f MB for particles.\n", *N * sizeof(struct part) / (1024.*1024.));
+
+ /* Open SPH particles group */
+ printf("read_ic: Reading particle arrays...\n");
+ h_grp = H5Gopen(h_file, "/PartType0");
+ if(h_grp < 0)
+ error( "Error while opening particle group.\n");
+
+ /* Read arrays */
+ readArray(h_grp, "Coordinates", DOUBLE, *N, 3, *parts, x);
+ readArray(h_grp, "Velocity", FLOAT, *N, 3, *parts, v);
+ readArray(h_grp, "Mass", FLOAT, *N, 1, *parts, mass);
+ readArray(h_grp, "SmoothingLength", FLOAT, *N, 1, *parts, h);
+ readArray(h_grp, "InternalEnergy", FLOAT, *N, 1, *parts, u);
+ readArray(h_grp, "ParticleIDs", ULONG, *N, 1, *parts, id);
+
+ /* Close particle group */
+ H5Gclose(h_grp);
+
+ printf("read_ic: Done Reading particles...\n");
+
+ /* Close file */
+ H5Fclose(h_file);
}
#endif /* HAVE_HDF5 */
+
diff --git a/src/ic.h b/src/ic.h
index c604b9a66ba7de2a2a4558468cca55ce9d140750..10216bfb8c410695b6aadc14792e0f072be8d5a5 100644
--- a/src/ic.h
+++ b/src/ic.h
@@ -18,5 +18,8 @@
******************************************************************************/
-void read_ic ( char* fileName, struct space *s, int* N);
-void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max );
+
+
+void read_ic ( char* fileName, double dim[3], struct part **parts, int* N, int* periodic);
+
+